ON THE FEEDING ZONE OF PLANETESIMAL FORMATION BY THE STREAMING INSTABILITY

For simplicity, we only consider a non-magnetized gas disk such that the magneto-rotational instability (e.g., Balbus & Hawley 1991) is not operating and thus the basic state of the gas is laminar. We also assume an isothermal equation of state, which remains a good approximation given that the flow is strongly subsonic, and hence that any small temperature increase is radiated away efficiently.

The equations of motion for the gas then become

$$\begin{eqnarray} \frac{\partial \rho _g}{\partial t} + u_{0,y}\frac{\partial \rho _g}{\partial y} + {\boldsymbol {\nabla }}\cdot (\rho _g{\boldsymbol {u}}) &=& 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} && \frac{\partial {\boldsymbol {u}}}{\partial t} + u_{0,y} \frac{\partial {\boldsymbol {u}}}{\partial y} + {\boldsymbol {u}}\cdot {\boldsymbol {\nabla }}{\boldsymbol {u}} = -c_s^2{\boldsymbol {\nabla }}\ln \rho _g \nonumber\\ &&\quad + \left(2\Omega _K u_y\hat{\boldsymbol {x}} - \frac{1}{2}\Omega _K u_x\hat{\boldsymbol {y}} - \Omega _K^2 z\hat{\boldsymbol {z}}\right) + 2\Omega _K\Delta v{\boldsymbol {\hat{x}}} + \frac{\rho _p}{\rho _g} \frac{{\boldsymbol {v}} - {\boldsymbol {u}}}{t_s}. \nonumber\\ \end{eqnarray} \tag{ 2 }$$

The dependent variables we solve for are the gas density ρ_g and the gas velocity ${\boldsymbol {u}}$ relative to the background shear flow ${\boldsymbol {u}}_0 = -(3/2)\Omega _K x {\boldsymbol {\hat{y}}}$, with Ω_K being the local Keplerian angular frequency. The parameter c_s is the speed of sound, being constant. The first two terms in the parenthesis of Equation (2) result from the combination of the stellar radial gravity and the centrifugal and the Coriolis forces, while the last term accounts for the linearized stellar vertical gravity. The next-to-last term on the right-hand side of Equation (2) resembles the centrifugal support due to the large-scale radial pressure gradient in the gas disk such that the azimuthal speed of the gas is reduced by approximately

$$\begin{equation} \Delta v = \frac{1}{2}c_s\left(\frac{c_s}{v_K}\right) \left(-\frac{\partial \ln {p}}{\partial \ln R}\right), \end{equation} \tag{ 3 }$$

in which v_K is the Keplerian velocity, p is the gas pressure, and R is the radial distance to the central star.¹ In this work, we set Δv/c_s = 0.05, which is a typical value in the inner regions of a minimum-mass-solar-nebular disk model (Hayashi 1981; Bai & Stone 2010b). The gas experiences the frictional drag from the solid particles through the last term in Equation (2), where ρ_p is the volume density of solids, ${\boldsymbol {v}}$ is the local velocity of particles relative to the background shear, and t_s is the stopping time for the particles (see Section 2.2). The factor ρ_p/ρ_g stems from the conservation of linear momentum in the friction between the gas and the particles.

The primary objective of this work is to simulate the streaming instability in relatively large boxes. In this regard, the vertical density stratification of the gas becomes significant, and preserving this stratification numerically is to our advantage here. We define

$$\begin{equation} \rho _g(x,y,z,t) \equiv \rho _{g,0}(z) [1 + \xi (x,y,z,t)], \end{equation} \tag{ 4 }$$

where

$$\begin{equation} \rho _{g,0}(z) = \rho _0 \exp \left(-\frac{z^2}{2H^2}\right) \end{equation} \tag{ 5 }$$

is a constant background density stratification, which is set by the balance between stellar vertical gravity and gas pressure, and the gas scale height is thus H = c_s/Ω_K. The arbitrary constant ρ₀ is the mid-plane density of this equilibrium stratification, which depends on the location of the shearing box in the protoplanetary disk. Note that positive densities imply ξ > −1. The equations of motion formulated in ξ now read

$$\begin{eqnarray} \frac{\partial \xi }{\partial t} + u_{0,y}\frac{\partial \xi }{\partial y} + {\boldsymbol {u}}\cdot {\boldsymbol {\nabla }}\xi + (1 + \xi){\boldsymbol {\nabla }}\cdot {\boldsymbol {u}} &=& \frac{z u_z}{H^2}(1 + \xi), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} && \frac{\partial {\boldsymbol {u}}}{\partial t} + u_{0,y}\frac{\partial {\boldsymbol {u}}}{\partial y} + {\boldsymbol {u}}\cdot {\boldsymbol {\nabla }}{\boldsymbol {u}} = -c_s^2{\boldsymbol {\nabla }}\ln (1 + \xi) \nonumber\\ &&\quad + \left(2\Omega _K u_y\hat{\boldsymbol {x}} - \frac{1}{2}\Omega _K u_x\hat{\boldsymbol {y}}\right) + 2\Omega _K\Delta v{\boldsymbol {\hat{x}}} + \frac{\rho _p}{\rho _g} \frac{{\boldsymbol {v}} - {\boldsymbol {u}}}{t_s}. \quad \end{eqnarray} \tag{ 7 }$$

The stellar vertical gravity exactly cancels the pressure gradient of the equilibrium density stratification in the momentum Equation (7), while a source term appears in the continuity Equation (6) to account for the imbalance in vertical mass flux due to the stratification. With this formulation, we solve for the dimensionless variable ξ instead of ρ_g.

We use the usual sheared periodic boundary conditions (e.g., Hawley et al. 1995), where f(x, y, z) = f(x + L_x, y − (3/2)Ω_KL_xt, z) in which f is any field in question and L_x is the x-dimension of the computational domain. We adopt periodic boundary conditions for the vertical direction. Note that in the formulation of Equations (6) and (7), a periodic ξ in the vertical direction does not introduce discontinuity in either density or pressure gradient (see, Davis et al. 2010). We set ξ = 0 and ${\boldsymbol {u}} = 0$ at t = 0 as our initial conditions.

Instead of treating particles as pressureless fluid, we consider the motion of each individual particle according to

$$\begin{eqnarray} \frac{d{\boldsymbol {x}}_p}{dt} &=& -\frac{3}{2}\Omega _K x_p {\boldsymbol {\hat{y}}} + {\boldsymbol {v}}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \frac{d{\boldsymbol {v}}}{dt} &=& \left(2\Omega _K v_y\hat{\boldsymbol {x}} - \frac{1}{2}\Omega _K v_x\hat{\boldsymbol {y}} - \Omega _K^2 z_p {\boldsymbol {\hat{z}}}\right) + \frac{{\boldsymbol {u}} - {\boldsymbol {v}}}{t_s}, \end{eqnarray} \tag{ 9 }$$

where ${\boldsymbol {x}}_p = (x_p, y_p, z_p)$ is the position of the particle relative to the center of the box and ${\boldsymbol {v}} = (v_x, v_y, v_z)$ is the relative velocity of the particle with respect to the Keplerian shear as defined above. Equation (8) is the total velocity of the particle while Equation (9) is the acceleration of the particle with the contributions parallel to those for the gas except the pressure gradient (see Equation (2)).

The stopping time t_s in Equations (2), (7), and (9) is the damping time for the relative speed between gas and each solid particle due to their mutual viscous drag. It is often expressed in terms of the dimensionless parameter τ_s ≡ Ω_Kt_s, which is a measure of the Stokes number for the particles. In this work, we set τ_s = π/10 ≃ 0.314, for which the radius of the particles is about 0.7 m at 1 AU or about 4 mm at 30 AU in the minimum-mass solar nebula (Johansen et al. 2007, 2014; Bai & Stone 2010b).

It is impractical to simulate all millimeter-to-meter-sized solid particles even with a computational box as small as one of 0.2H on an edge. Instead, we consider super-particles, each of which represents a swarm of real, identical particles. The mass of each super-particle is the total mass of the constituent particles, while the damping time for the super-particle remains the same as that of the individual members. It has been demonstrated that this approach is numerically convergent when on average more than one super-particles per grid cell exist in the sedimented mid-plane layer (Youdin & Johansen 2007; Bai & Stone 2010c).

The mass of each super-particle m_p is determined by the initial solid-to-gas ratio, Z ≡ Σ_{p, 0}/Σ_{g, 0}, of the medium, where Σ_{p, 0} and Σ_{g, 0} are the initial (uniform) column densities of the particles and the gas, respectively, integrated over the full vertical extent of the disk. From Equation (5), $\Sigma _g = \sqrt{2\pi }H\rho _0$. Since virtually all particles we consider settle within a distance much less than H of the mid-plane of the disk, Σ_p is well approximated by m_pN_p/(L_xL_y), where N_p is the total number of particles used in a simulation, and L_x and L_y are the sizes of the computational domain in the x and the y directions, respectively. In this work, we set Z = 0.02, which is just above the critical solid-to-gas ratio required to trigger strong clumping of particles by the streaming instability (Johansen et al. 2009). In addition, we use as many particles as the total number of grid points in each simulation, i.e., N_p = N_xN_yN_z, where N_x, N_y, and N_z are the number of grid points in the x, y, and z directions. Therefore, the average number of particles per cell near the mid-plane is roughly N_z/N_mid after sedimentation, where N_mid is the number of vertical grid cells resolving the particle scale height.

As our initial conditions, we use a uniform distribution to randomly place the particles throughout the computational domain while setting ${\boldsymbol {v}} = 0$. The noise inherent in the initial positions of the particles serve as the seed for the ensuing growth of the streaming instability. The boundary conditions for the particles are such that when a particle crosses a boundary plane, it reemerges in the opposite plane with sheared periodic positional mapping as the gas while preserving its relative velocity ${\boldsymbol {v}}$.

To solve the system of Equations (6)–(9), we use the Pencil Code,² a cache-efficient, parallelized magnetohydrodynamical code founded by Brandenburg & Dobler (2002). For the (magneto-)hydrodynamics, the code employs sixth-order finite differences in space while integrating the system of equations in time by third-order Runge–Kutta steps. For the particle dynamics, the position and velocity of each individual particle is evolved simultaneously with the Runge–Kutta time steps for the fluid. The interactions between the fluid and the particles, i.e., the frictional drag, are computed via the standard particle-mesh method of triangular-shaped clouds, which ensures conservation of total momentum (Youdin & Johansen 2007; Johansen et al. 2007). The particle-block-decomposition algorithm implemented by Johansen et al. (2011) is used in order to achieve better load balance among processors.

We have implemented in the Pencil Code the new formalism for balanced stratification of gas density, Equations (6) and (7). We find this formalism is in general more numerically stable than evolving the system of Equations (1) and (2) in the sense that much less artificial diffusion is required to maintain hydrostatic equilibrium against perturbations. This is especially true when the vertical extent of the computational domain is greater than a few disk scale heights. Furthermore, this formalism relieves the necessity of implementing special boundary conditions in order to incorporate the background density stratification. Finally, we find that the total mass in the system remains well conserved with this new formalism, despite the continuity equation (Equation (6)) not being written in conservative form.

In this work, we also increase the time-step constraint limited by the drag force calculation from one-fifth of the decay time constant adopted in previous works to one time constant. This amounts to a possible over-damping in relative velocity between gas and particles by a maximum relative error of about 9%. We find that this relaxation does not noticeably alter the results while allowing us a five-time speed up for otherwise the same simulation.

We systematically adjust the dimensions of the simulation box and investigate the difference in the results. The horizontal sizes L_x = L_y span from 0.2H up to 1.6H, while the vertical size L_z is up to 1.6H with L_z ⩽ L_x. The maximum number of grid cells in each dimension we have explored is 256, which translates to a resolution of 160 H⁻¹ for our largest 1.6H × 1.6H × 1.6H box.

The top row of Figure 1 shows the particle scale height as a function of time for boxes with various resolutions and horizontal sizes but a constant vertical size of L_z = 0.2H. Since the particles are initially uniformly distributed, the initial scale height is virtually infinite. Nevertheless, the particles quickly settle down toward the mid-plane within t = 2P due to their vertical motion and gas drag, where P = 2π/Ω_K is the orbital period. At this point, the particles have concentrated near the mid-plane to the extent that the gas experiences enough perturbation from the particles to become turbulent. The random motion of the gas-dragged particles becomes dominant and puffs the particle layer back up. Within another time interval of ∼1P, the particle scale height reaches its final, roughly constant value, which lasts till the end of the simulations at t = 100P.

Zoom In Zoom Out Reset image size

Little variation in the equilibrium particle scale height is observed between boxes of different resolutions and/or different horizontal sizes. Only the boxes with a resolution of 40H⁻¹, the lowest we have investigated in this work, show a slightly larger particle scale height. The equilibrium scale height is on the order of ∼10⁻²H, which is resolved beyond a resolution of 320H⁻¹. Nevertheless, it is not clear if resolving the particle layer is critical in the saturated state of the streaming instability except perhaps in predicting the correct peak local particle density (see below).

The top row of Figure 2 also shows the particle scale height as a function of time, but for boxes of various horizontal and vertical dimensions at a fixed resolution of 160H⁻¹. Similarly to the case of L_z = 0.2H discussed above, the horizontal size has little effect on the particle scale height for other vertical dimensions. On the other hand, we see a factor of close to two increase in the particle scale height from a box of L_z = 0.2H to that of L_z = 0.4H. Boxes of L_z ≳ 0.4H have a consistent equilibrium value. However, it appears that the larger the vertical size of the box, the longer timescale is required to reach the equilibrium, which might be an effect of the more stochastic particle clumping for boxes with larger L_z (see Section 3.3).

Zoom In Zoom Out Reset image size

The bottom rows of Figures 1 and 2 show the corresponding maximum particle density ρ_{p, max} as a function of time for various box dimensions and resolutions. Due to the vertical sedimentation of the particles, the particle density coherently increases within t = 2P to above the initial gas density ρ₀ in the mid-plane. As can be seen in Figure 1, the particle density for the boxes with a resolution of 40H⁻¹ then remains at a constant low level of a few ρ₀ and no significant particle concentration occurs. On the other hand, appreciable particle concentration driven by the streaming instability starts to appear and drives ρ_{p, max} further up for all boxes with a resolution of ≳80H⁻¹, before it reaches a roughly constant state at t ∼ 20P.

In general, the higher the resolution, the larger the maximum particle density ρ_{p, max} results, which has been reported in previous studies (Johansen et al. 2007, 2012; Bai & Stone 2010c). Here, we also find the higher the resolution, the smaller the increase in the final level of ρ_{p, max}. This indicates the numerical convergence with resolution in the saturated stage of the streaming instability. Also, Figure 1 hints that a resolution of ∼160–320H⁻¹ might already give a converged result in ρ_{p, max}, at least in the case of this work. A level of ρ_{p, max} close to 10³ρ₀ has been reached.

As shown in Figure 2, we find little variation in the maximum particle density ρ_{p, max} in the saturated stage of the streaming instability with respect to the vertical box size L_z. On the other hand, there exists a slight increase of a factor of a few in ρ_{p, max} from horizontal box size L_x = L_y = 0.2H to L_x = L_y = 0.4H, while the results are consistent for all boxes with L_x = L_y ≳ 0.4H.

Combining with the similar behavior of the particle scale height discussed in Section 3.1, this suggests that simulation boxes with either horizontal extent L_x = L_y = 0.2H or vertical extent L_z = 0.2H are insufficient to capture all necessary scales perturbed by the streaming instability in its nonlinear saturation stage. More evidence on this is presented in Section 3.3.

The streaming instability predominantly concentrates sedimented particles radially into filamentary structures, extended in the azimuthal direction (Johansen & Youdin 2007; Bai & Stone 2010b; Kowalik et al. 2013). Therefore, the column density of particles, while averaged over the azimuthal dimension of the simulation box, well describes the time evolution of the particle layer driven by the streaming instability. Particular interest here is to investigate the dependence of these particle radial concentrations on the box dimensions as well as the resolution.

Figure 3 shows the averaged column density of particles 〈Σ_p〉 as a function of time and radial position for the simulation boxes with the same resolution (160H⁻¹) and vertical extent (0.2H) but varying horizontal dimensions (from 0.2H to 1.6H). The 0.2H × 0.2H × 0.2H box shows only one major filament of solids, which is consistent with previous works (e.g., Johansen et al. 2009). The 0.4H × 0.4H × 0.2H box, however, generates two filaments initially, but one of them disperses and later merges into the other, forming one dense filament. More interestingly, the 0.8H × 0.8H × 0.2H box and the 1.6H × 1.6H × 0.2H box demonstrate much richer dynamics driven by the streaming instability, with multiple filaments of solids forming, dispersing, splitting, and merging, in drastic contrast to one single dominant filament for smaller simulation boxes. For the 1.6H × 1.6H × 0.2H box, about five or six dense filaments coexist at any given time. Noticeable is that their separation remains quite regular for a long period of time, due to roughly the same radial drift speed of the filaments. Since the radial drift speed is determined by the particle density (Nakagawa et al. 1986; Youdin & Goodman 2005), this in turn implies that the filaments have roughly the same column density, as is also seen in the figure.

Zoom In Zoom Out Reset image size

Figure 4 compares 〈Σ_p〉 for the simulation boxes of various vertical sizes (from L_z = 0.2H to 1.6H) but the same resolution (160H⁻¹) and horizontal dimensions (L_x = L_y = 1.6H). Extending the vertical size of the box evidently introduces much more complexity in the evolution of the particle layer. The densities of the particle filaments for the tall boxes have significantly larger variance than those for their short counterpart (the 1.6H × 1.6H × 0.2H box) such that their radial drift speeds noticeably differ. This in turn makes the merging and splitting events of the filaments occur relatively more frequently. Nevertheless, multiple filaments still remain at any given time for these tall boxes.

Zoom In Zoom Out Reset image size

In order to make a more quantitative statement on their characteristics, we devise a simple algorithm to capture the particle filaments. At any given time, we use a stencil of fixed physical length w to scan through the radial position x, appending enough ghost cells near the two ends of the computational domain. If the azimuthally averaged particle column density 〈Σ_p〉 at the center of a stencil is the maximum for all points in the stencil and is larger than a certain threshold, we define that a concentrated particle filament occurs at this location with its peak density the same as the maximum 〈Σ_p〉. We choose w = 0.05H and use for the threshold 5σ Poisson noise in a uniform particle distribution, which is represented by N_z particles for each column of cells in our simulations, i.e.,

$$\begin{equation} \max _{x_i - w/2 \le x_j \le x_i + w/2}\langle \Sigma _p\rangle _{x = x_j} = \langle \Sigma _p\rangle _{x = x_i} > \Sigma _{p,0}\left(1 + \frac{5}{\sqrt{N_z}}\right). \end{equation} \tag{ 10 }$$

Note that the minimum separation between adjacent filaments that can be resolved is then w/2 = 0.025H. In our simulations, especially those with high resolutions (≳160H⁻¹), we do find that dense filaments undergo splitting and merging events with structures of length scale less than 0.025H. However, the majority of these events are intermittent and the filaments recover their original states on a relatively short timescale. The remaining ones do significantly change the properties of the filaments and are detectable by this simple algorithm.

Figure 5 shows the mean separation D and the mean peak averaged column density $\overline{\langle \Sigma _p\rangle _{\max }}$ of the particle filaments as a function of time for simulation boxes of various resolutions and horizontal dimensions but fixed vertical size L_z = 0.2H. For boxes of the same dimensions, increasing the resolution generally results in more filaments being identified and thus smaller mean separation. As can be seen in the figure, we find that the mean separation tends to converge with resolution toward D ∼ 0.2H, which is well above the radius of the stencil we use. For smaller boxes with L_x = L_y = 0.2H and 0.4H, the mean peak density $\overline{\langle \Sigma _p\rangle _{\max }}$ covers a wide range of values (between 0.05Σ_{g, 0} and 0.35Σ_{g, 0}) with no obvious trend of convergence. This is due to small-number statistics (only one or two filaments form in these cases) and indicates the stochastic nature in the formation and evolution of these filaments. On the other hand, the larger boxes with L_x = L_y = 0.8H and 1.6H does show relatively consistent results between different resolutions and horizontal box sizes, where $\overline{\langle \Sigma _p\rangle _{\max }} \sim 0.08\Sigma _{g,0}$.

Zoom In Zoom Out Reset image size

Figure 6 further compares the mean separation D and the mean peak averaged column density $\overline{\langle \Sigma _p\rangle _{\max }}$ for boxes of various vertical size L_z. As discussed above, increasing L_z makes the formation and evolution of the particle filaments even more stochastic. This behavior is also manifest in D and $\overline{\langle \Sigma _p\rangle _{\max }}$, where more time variation with more significant amplitude occurs for boxes of larger L_z. This is yet another view of the more frequent formation, dispersal, splitting, and merging of the particle filaments along with larger variance in their densities for taller boxes. Nevertheless, D oscillates around ∼0.2H for our largest boxes, indicating that the number of the major, persistent filaments is about the same over time. These results demonstrate the importance of the vertical coverage of the gas disk to fully capture the dynamics driven by the streaming instability, even though the particle layer is thin compared to the gas scale height.

Zoom In Zoom Out Reset image size

From Figures 5 and 6, we also notice that there exists a trend of increasing mean separation D and the mean peak averaged column density $\overline{\langle \Sigma _p\rangle _{\max }}$ at late times for various cases. This indicates the tendency for the particle filaments to merge and keep accreting surrounding materials. However, we note that this phenomenon might not be relevant in reality, since from Figures 1 and 2, the local particle density, being the decisive factor to drive the ultimate gravitational collapse of the particles to form planetesimals, would have reached its maximum level already in the early saturated stage of the streaming instability.

The results presented in Section 3.3 indicate that even though the particle layer is thin compared with the gas disk, the gas dynamics over at least one gas scale height cannot be ignored in the nonlinear stage of the streaming instability. It appears that a significant fraction of the column of the gas mass is still required to better describe the interaction between the gas and the solids. In this regard, we need to study the correlation between the gas and the solids and its dependence on the vertical dimension of the simulation box.

We plot in Figure 7 the correlation coefficient between the yz averages of the gas density deviation ξ and of the solid density ρ_p for various vertical box sizes. The former is a proxy for the azimuthal average of the gas column density 〈Σ_g〉 and the latter is directly proportional to that of the solid column density 〈Σ_p〉, both of which are functions of the radial position x and time t. Since most of the solid particles are concentrated in azimuthal filamentary structures, the correlation coefficient represent the spatial correlation between the gas mass and the particle filaments. As shown in the figure, although there exists significant time variation in the correlation coefficient, the gas distribution and the particle distribution in general anti-correlate, with an average of about −0.2. This implies that the gas tends to be entrapped in between the particle filaments and slightly enhance the pressure there. We also see that this anti-correlation decreases noticeably with increasing vertical box size, indicating the lessening of the pressure buildup with increased vertical dynamical range.

Zoom In Zoom Out Reset image size

Given these findings, we speculate that even though the perturbation in the gas due to the streaming instability is only at about 0.1% level, it is enough for the gas to participate in regulating the dynamics of the particle filaments. The slightly enhanced gas pressure between the filaments may help inhibit the filaments from approaching, leading to the exceptionally regular spacing and similar migration speeds as we see in Figure 3 for short simulation boxes. With increased vertical dynamical range, on the other hand, the gas in the mid-plane gains freedom to escape vertically when the particle filaments tend to merge or split, relieving otherwise the pressure enhancement of the gas. This may explain why there exists a noticeable decrease in the anti-correlation between the gas and the solid column densities with increasing vertical box size, as seen in Figure 7. In any case, this further demonstrates the importance of the vertical dimension in the nonlinear evolution driven by the streaming instability.