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, short-lived, or have too much radial drift. We chose 25 orbits because it is half of the 50-orbit 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 short-lived 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 large-scale measure of clumping. The Kolmogorov-Smirnov (Smirnov 1948; Wall & Jenkins 2003) test gives a standard way to compare F and G. It computes a p-value 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 p-values as strict probabilities; we only use p as a convenient metric to measure clumping. With that in mind, we divide the p-values 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 chi-squared test on ⟨ ρ_p ⟩ _zt, without sorting the grid cells. We use the reduced chi-squared 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 σ² 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 σ². This can be because the model is incorrect, or because the true variance between measurements is greater than σ². We divide into four intervals:

• Clumping is unlikely:

  

• Clumping is somewhat likely:

  

• 

• Clumping is very likely:

  

That leaves σ² as the only unspecified parameter. When the simulations start, σ² is almost zero, giving unrealistically high values. As sedimentation proceeds, the value of σ² increases. We opted for a σ² that is typical for the last 25 orbits of the sedimentation phase. Because is sensitive to σ², there is necessarily some amount of subjectivity in σ². In the end, we chose σ² 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 25-orbit 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 KS-derived 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: Three-dimensional 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 (x-z 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 Kelvin-Helmholtz Instability (KHI). Whether the KHI is present is dictated by the Richardson number (Chandrasekhar 1961), where u is the gas speed, g = Ω²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 long-term 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 higher-resolution runs produce results consistent with the 128 × 128 runs.

Open with DEXTER

© ESO, 2015