Online material
Appendix A: A measure of particle clumping
Particle clumping is difficult to measure objectively. Previous authors have measured the peak particle density, but this is a very local measure that contains limited information about the particle behavior. In particular when small particles (τ_{f} ≤ 0.003) form filaments that are clearly visible in a spacetime diagram (Fig. 4), no single grid cell achieves a very high density. At the same time, large particles (τ_{f}> 3) have may have a few cells with very high density, but the density peaks are unstable and short lived.
Figure A.1 illustrates our approach. For each time step, we take the particle surface density Σ_{solid} (top row). Particle overdensities are clearly visible on this plot. We then average Σ_{solid} over 25 orbits (⟨ Σ_{solid} ⟩ _{t}, second row). This serves to to smooth out overdensities that are unstable, shortlived, or have too much radial drift. We chose 25 orbits because it is half of the 50orbit timescale for Z to increase. We considered using the radial speed to separate clumps that are drifting into the star. The appeal of radial speed is that it provides an instantaneous measure. However, radial speed does not capture the longevity of a particle clump. For example, as in a traffic jam, the pattern speed of the clump may be lower than the particle drift speed. At the same time, particles with low radial speed may produce shortlived clumps. Therefore, it is necessary to compare the location of a clump across time. Averaging Σ_{solid} over time is a simple way to achieve this goal.
Finally, we sort the grid cells of ⟨ Σ_{solid} ⟩ _{t} in order from highest to lowest density, and compute the cumulative distribution F of the sorted grid cells. If the particle distribution is uniform, then F forms a straight diagonal line G from (0, 0) to (1, 1). We use the difference F−G as a largescale measure of clumping. The KolmogorovSmirnov (Smirnov 1948; Wall & Jenkins 2003) test gives a standard way to compare F and G. It computes a pvalue that measures probability that the data set that produced F follows the distribution Gwhere D is the maximum distance between F and G, and n is the number of independent observations. We caution the reader against interpreting these pvalues as strict probabilities; we only use p as a convenient metric to measure clumping. With that in mind, we divide the pvalues into four broad categories:

Clumping is unlikely:
p ≥ 0.45

Clumping is somewhat likely:
0.25 ≤ p< 0.45

Clumping is likely:
0.10 ≤ p< 0.25

Clumping is very likely:
p< 0.10.
Figure 8 shows the result of this analysis across our entire range of simulations. The figure is broadly consistent with the spacetime diagrams in Figs. 4 and 5.
Appendix A.1: Alternate measures of particle clumping
To confirm our results, we devised an alternate method to measure particle clumping, based as much as possible on very different principles. Our second method is to perform a chisquared test on ⟨ ρ_{p} ⟩ _{zt}, without sorting the grid cells. We use the reduced chisquared statistic, (A.3)where k = n − 1 is the degrees of freedom, O_{j} is the jth observation (value of ⟨ ρ_{p} ⟩ _{zt} on grid cell j), E_{j} = ⟨ ρ_{p} ⟩ is the jth expected density for a uniform distribution, and σ^{2} is the expected variance between different measurements. If , that indicates that the difference between observations and the model is greater than can be explained by the variance σ^{2}. This can be because the model is incorrect, or because the true variance between measurements is greater than σ^{2}. We divide into four intervals:

Clumping is unlikely:

Clumping is somewhat likely:


Clumping is very likely:
That leaves σ^{2} as the only unspecified parameter. When the simulations start, σ^{2} is almost zero, giving unrealistically high values. As sedimentation proceeds, the value of σ^{2} increases. We opted for a σ^{2} that is typical for the last 25 orbits of the sedimentation phase. Because is sensitive to σ^{2}, there is necessarily some amount of subjectivity in σ^{2}. In the end, we chose σ^{2} so that the τ_{f} = 0.003 simulations wold give the same result as in Fig. 8. Therefore, the true test is whether the method remains consistent with the KS method across all the other particle sizes.
Figure A.2 shows the final results for the method. The similarity between this figure and Fig. 8 is striking. Considering how different the two techniques are, the agreement between Figs. 8 and A.2 indicates that our core results are robust.
Fig. A.1
Particle density distribution for three particle sizes. Particle size is measured in the stopping time τ_{f}, where τ_{f} = 0.001, 0.003 and 3 correspond to the suspension regime, streaming regime and radial drift regime respectively. Row a) is the Surface density of solids Σ_{solid} as a function of the radial coordinate x, evaluated at two points in time 25 orbits apart (first blue, then green). Row b) averages Σ_{solid} over the 25orbit interval ⟨ Σ_{solid} ⟩ _{t}. Row c) shows ⟨ Σ_{solid} ⟩ _{t} with the grid cells sorted from highest to lowest density. A steep curve indicates stable particle clumps. Row d) shows the cumulative distribution of ⟨ Σ_{solid} ⟩ _{t} in row c) (solid red curve) along with the cumulative distribution of a perfectly uniform particle distribution (dashed blue line). The distance between these two curves serves as a measure of particle clumping. In these plots the particle concentrations Z = ⟨ Σ_{solid} ⟩ / ⟨ Σ_{total} ⟩ are (left to right) 0.095, 0.095 and 0.05. 

Open with DEXTER 
Fig. A.2
As in Fig. 8, but using the method to estimate the likelihood of clumping instead of the KSderived method. The figure marks the region of the particle size vs. concentration phase space where the streaming instability is active. Particle size is measured in the stopping time τ_{f} and the particle concentration is Z = Σ_{solid}/ Σ_{total}. The colors mark the probability that particle clumps can form, where red is “unlikely”, magenta is “somewhat likely”, blue is “likely”, and green is “very likely”. When different simulations give different results, the two extreme values are shown. The symbols indicate the number of simulations available. The circle, cross and triangle indicate, respectively, one, two and three simulations. Regions without symbols are extrapolations. To facilitate the comparison with Fig. 8, the boundary marked by the black lines has been copied from Fig. 8 without modification. 

Open with DEXTER 
Appendix B: Convergence tests
Appendix B.1: Threedimensional box
Fig. B.1
Spacetime diagram of two 3D runs with particle size τ_{f} = 0.03 and Δ = 0.05. The color indicates the column density Σ_{solid}/ ⟨ Σ_{solid} ⟩ using the same color bar as Figs. 4 and 5. Both runs begin with Z = 0.01 and have Z increased only for a short interval. In one run (left), Z grows to 0.02 and in the other (right) Z grows to 0.03. Clumps are visible for Z ≥ 0.02, consistent with the 2D runs. 

Open with DEXTER 
Fig. B.2
A snapshot of run 3D.Z2 taken at t = 150 orbits (Z = 0.02). Image a) is the view from the side (x − z plane). Image b) is the view from above (x − y plane). The color scale marks the column solid density normalized to the initial gas density (Σ_{solid}/ Σ_{gas,0}). The two images have a different color scale, since the column density in the azimuthal direction is much higher. Note that the particle structure is nearly axisymmetric. 

Open with DEXTER 
Fig. B.3
A sequence of vertical snapshots (xz plane) showing the formation of distinct particle clumps in one of our 2D simulations (left) and our two 3D simulations (middle, right). For each run, the images are taken 10 orbits apart. The time and Z values are shown in each plot. The times were chosen to show the formation of the first distinct clumps, and to give a similar peak density on the final image (ρ_{p,max} = 1.8,1.5,2.1left to right). The color indicates the column density using the same color scale as in Fig. B.2a. The other 2D runs give similar results to the one shown. The 2D run is consistent with the 3D runs. In particular, the 3D runs do not show any evidence of additional turbulence stirring compared to the 2D runs. 

Open with DEXTER 
Fig. B.4
Maximum particle density (ρ_{p,max}) for one of the 2D simulations and the two 3D simulations. The top plot compares run 3D.Z2 (blue) and a 2D run (red). The bottom plot compares run 3D.Z3 (blue) and the same 2D run (red). In each plot the purple line shows ρ_{p,max} for the 3D run after averaging along the azimuthal direction. The vertical dashed lines mark the places where the 3D runs transition from Z = 0.01 to Z = 0.019 (for 3D.Z2) or Z = 0.029 (for 3D.Z3). Similarly, the vertical solid lines mark the progressive increase of Z during the course of the 2D runs. 

Open with DEXTER 
Our simulations are axisymmetric, they neglect the effect of the KelvinHelmholtz Instability (KHI). Whether the KHI is present is dictated by the Richardson number (Chandrasekhar 1961), where u is the gas speed, g = Ω^{2}z is the vertical gravitational acceleration, and ρ is the effective fluid density (treating gas and solids as a single fluid). Bai & Stone (2010a) have shown that the streaming instability is able to keep the Richardson number above the critical value needed for the KHI, so that the KHI is absent in this problem. To confirm this, we performed two 3D simulations with τ_{f} = 0.03 and Δ = 0.05. The experiment performed in our 3D simulation is slightly different from that of our 2D simulations. Instead of increasing Z gradually, we only test the behavior of the system at three discrete values: Z = 0.01, Z = 0.02 and Z = 0.03. The reason for this is that 3D simulations are very computationally expensive, especially for small particles, so that a longterm simulation is prohibitive. Our 3D runs have three phases:

Phase I: each run starts with Z = 0.01, and Z is held constant for 30 orbits.

Phase II: the value of Z is increased quickly. In run “3D.Z2” we increase Z to Z = 0.02, and in run “3D.Z3” we increase Z to Z = 0.03.

Phase III: the value of Z is once again held constant for the remainder of the run.
Figure B.1 shows the spacetime diagrams for both runs. The figure shows that both runs produce visible particle clumps for Z ≥ 0.02 and no clumps for Z = 0.01. It is also clear that clumps form more readily with Z = 0.03. These results are consistent with our 2D simulations (Fig. 4). For particles smaller than τ_{f} = 0.03, the small integration steps make 3D simulations prohibitive. However, as explained in Sect. 2.6, the long response time of very small particles (τ_{f} ≤ 0.0064) means that our simulations will overestimate Z; so our results can be considered robust upper bounds for this size range.
Figure B.2 gives a side view (x − z axis) and a top view (x − y axis) of the 3D.Z2 run at t = 150 orbits. The clumps and their filamentary structure are very prominent. An important feature of the 3D runs is that they form visible clumps more easily than the 2D runs. Figure B.3 shows the formation of clumps in the two 3D runs and one of the 2D runs. This figure shows that the 3D runs are consistent with the 2D run. In particular, the clumps are readily visible at Z = 0.02 and the 3D runs show no evidence of additional stirring from the KHI.
Figure B.4 shows the maximum particle density (ρ_{p,max}) for the two 3D simulations and one of the 2D simulations. One salient feature of the figure is that in the early phase of the simulations, before clumping occurs, the 3D runs consistently show a peak density 2−3 times greater than the 2D runs. This occurs because particles can also accumulate along the azimuthal direction. The figure also shows ρ_{p,max} for the 3D runs after averaging along the azimuthal direction (purple lines). This averaging removes the apparent discrepancy between the 2D and 3D runs. In the later stages the 2D runs reach higher peak densities because they have very high Z values.
Appendix B.2: Resolution
We performed two simulations at 256 × 256 grid resolution with τ_{f} = 0.003 particles. Figure B.5 shows the spacetime diagrams for both sets of simulations. Since the behavior of the 256 × 256 runs is almost identical to the corresponding 128 runs, we conclude that, in terms of resolution, our simulations are fully converged.
Fig. B.5
Spacetime diagrams showing the solid surface density Σ_{solid} as a function of the radial coordinate x and simulation time τ_{f} = 0.003 and τ_{f} = 0.010 particles. The surface density is shown by color, using the same color scale as in Figs. 4 and 5. The figure shows both 128 × 128 simulations and the corresponding 256 × 256 simulations. In all simulations Δ = 0.05. The higherresolution runs produce results consistent with the 128 × 128 runs. 

Open with DEXTER 
© ESO, 2015