On the Numerical Robustness of the Streaming Instability: Particle Concentration and Gas Dynamics in Protoplanetary Disks

We simulate a small patch of the protoplanetary disk using the local shearing box approximation, which is useful for studying dynamics on length scales smaller than the disk radius. With this approximation, the global cylindrical geometry of the disk is mapped onto local Cartesian coordinates with unit vectors $\hat{x}$, $\hat{y}$, and $\hat{z}$ in the radial, azimuthal, and vertical directions, respectively (Goldreich & Lynden-Bell 1965). The center of the computational domain—at (x, y, z) = (0, 0, 0)—is located at a fiducial disk radius, in the midplane. At this location, the (pressure-corrected; see below) Keplerian frequency is Ω₀. The reference frame orbits, and thus rotates, at this frequency. The local Keplerian velocity is v_K.

We solve the equations of continuity and motion for gas, as well as the equations of motion for each solid particle labeled by the subscript i:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{g}}}{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{(\partial {\rho }_{{\rm{g}}}{\boldsymbol{u}})}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{g}}}{\boldsymbol{u}}{\boldsymbol{u}}+P{\boldsymbol{I}})\\ &&\quad =\,{\rho }_{{\rm{g}}}[2{\boldsymbol{u}}\times {{\boldsymbol{\Omega }}}_{0}+3{{\rm{\Omega }}}_{0}^{2}{\boldsymbol{x}}-{{\rm{\Omega }}}_{0}^{2}{\boldsymbol{z}}]+{\rho }_{{\rm{p}}}\displaystyle \frac{\bar{{\boldsymbol{v}}}-{\boldsymbol{u}}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{v}}}_{i}}{{dt}}=2{{\boldsymbol{v}}}_{i}\times {{\boldsymbol{\Omega }}}_{0}+3{{\rm{\Omega }}}_{0}^{2}{{\boldsymbol{x}}}_{i}-{{\rm{\Omega }}}_{0}^{2}{{\boldsymbol{z}}}_{i}-\displaystyle \frac{{{\boldsymbol{v}}}_{i}-{\boldsymbol{u}}}{{t}_{\mathrm{stop}}}-2\eta {v}_{{\rm{K}}}{{\rm{\Omega }}}_{0}\hat{x},\end{eqnarray} \tag{ 3 }$$

where ρ_g, ${\boldsymbol{u}}$, and P are density, velocity, and pressure of gas; ${\boldsymbol{I}}$ is the identity matrix; and ${{\boldsymbol{\Omega }}}_{0}={{\rm{\Omega }}}_{0}\hat{z}$. For the solids, ρ_p and $\bar{{\boldsymbol{v}}}$ give the average density and velocity of the particles in a hydrodynamic grid cell; see Youdin & Johansen (2007) for details. For individual solid particles, the velocities, ${{\boldsymbol{v}}}_{i}$, and radial and vertical positions, ${{\boldsymbol{x}}}_{i}$ and ${{\boldsymbol{z}}}_{i}$, are indexed.

The stopping time, t_stop, measures how quickly drag forces decelerate a particle's motion relative to the gas. Our models consider a uniform value of the stopping time for all particles. This choice corresponds to a monodisperse distribution of particle size and mass and to the Stokes or Epstein drag accelerations, which are linear in relative velocity.

The forces on the gas in Equation (2) are the Coriolis, radial, and vertical tidal gravity (all in square brackets), with the final term giving the back reaction of drag forces on the gas, a crucial ingredient for collective drag effects such as the SI. The particle equation of motion, Equation (3), includes the corresponding acceleration terms for the solids, with the penultimate term the drag acceleration on a particle. The last term in Equation (3), $-2\eta {v}_{{\rm{K}}}{{\rm{\Omega }}}_{0}\hat{x}$, is a constant inward radial acceleration of particles. This term replaces the (equal but opposite) outward acceleration of the gas by global radial pressure gradients, which must be included by hand in a local model. Switching this acceleration to the particles means that Ω₀ is the pressure-supported, slightly sub-Keplerian orbital frequency of the gas.

We adopt an isothermal equation of state, $P={c}_{{\rm{s}}}^{2}{\rho }_{{\rm{g}}}$, where the sound speed, c_s, is constant. The units of the code are the natural units of the shearing box. The time unit is ${{\rm{\Omega }}}_{0}^{-1}$, the orbital timescale. The length unit is H = c_s/Ω₀, the vertical scale height of the gas. The mass unit is set by ρ_g,0, the midplane gas density in hydrostatic balance.

The radial and azimuthal boundary conditions are the standard shear periodic conditions for the shearing sheet, as implemented for ATHENA in Stone & Gardiner (2010). As described in Section 1, one of our goals is to study the effect of different vBCs.

The classic or standard BCs for SI simulations in a stratified shearing box are the shearing periodic boundaries in the radial direction, purely periodic in the azimuthal direction and periodic or reflecting in the vertical direction. While applying the periodic or reflecting vBCs, simulations are assuming there are an infinite amount of disks stacking up together with even spacing. Thus, gas wave motion owing to the interplay with the particle layer in the midplane will propagate up and down to affect other disks, which numerically is the disk itself. This configuration is not the favorable situation. Physically speaking, all truncated disks are artificial. However, due to the nature of numerical experiments, it is only computationally feasible to study SI within a local shearing box for now, which means there exist vertical limits based on box size. Hence, we propose to adopt a new vBC to improve gas behaviors.

Following the Appendix of Simon et al. (2011), we implement a modified outflow boundary condition to handle gas flow at the boundaries. It extrapolates gas density outside the upper and lower boundaries via an exponential function into ghost zones. Taking the upper ghost zone as an example, the gas density at Z_end + δZ is

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}({Z}_{\mathrm{end}}+\delta Z)={\rho }_{{\rm{g}}}({Z}_{\mathrm{end}})\exp \left(-\displaystyle \frac{{(\delta Z)}^{2}}{2{H}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where Z_end denotes the $\hat{z}$-coordinates of the physical cells right at the upper box boundary (let us call them parent cells) and δZ indicates the $\hat{z}$-coordinate excesses of the ghost cells outside the upper boundary. This extrapolation naturally creates the hydrostatic balance since the vertical gravity takes effect everywhere. Then for each cell in the ghost zone, its horizontal velocity is directly copied from the corresponding parent cell. But, its vertical velocity, u_z, will only be copied when the parent cell has a u_Z that does not induce inflow. Otherwise, the u_Z of the ghost cell will be set to zero to prevent introducing garbage information. In such a manner, the total amount of the materials in the numerical box will slowly monotonotically decrease, but the net outward flow is very tiny in each time step (see Table 1). Therefore, we renormalize the total mass in the entire domain after each time step to make it constant (Ogilvie 2012). This design robustly maintains the hydrostatic equilibrium of gas (see the Appendix).

Table 1.  Simulation Parameters

Runa	Domain Size	Resolution					Nparb	Vertical Fluxc	vBCs
	(LX × LY × LZ)H3	NX	×	NY	×	NZ		ρg,0cs
re22	0.2 × 0.2 × 0.2	64	×	64	×	64	219	⋯	Reflecting
re24	0.2 × 0.2 × 0.4	64	×	64	×	128	219	⋯	Reflecting
re28	0.2 × 0.2 × 0.8	64	×	64	×	256	219	⋯	Reflecting
re42	0.4 × 0.4 × 0.2	128	×	128	×	64	221	⋯	Reflecting
re44	0.4 × 0.4 × 0.4	128	×	128	×	128	221	⋯	Reflecting
re48	0.4 × 0.4 × 0.8	128	×	128	×	256	221	⋯	Reflecting
re82	0.8 × 0.8 × 0.2	256	×	256	×	64	223	⋯	Reflecting
re84	0.8 × 0.8 × 0.4	256	×	256	×	128	223	⋯	Reflecting
re88	0.8 × 0.8 × 0.8	256	×	256	×	256	223	⋯	Reflecting
pe22	0.2 × 0.2 × 0.2	64	×	64	×	64	219	6.42e–3	Periodic
pe24	0.2 × 0.2 × 0.4	64	×	64	×	128	219	4.54e–3	Periodic
pe28	0.2 × 0.2 × 0.8	64	×	64	×	256	219	2.46e–3	Periodic
pe42	0.4 × 0.4 × 0.2	128	×	128	×	64	221	5.52e–3	Periodic
pe44	0.4 × 0.4 × 0.4	128	×	128	×	128	221	5.03e–3	Periodic
pe48	0.4 × 0.4 × 0.8	128	×	128	×	256	221	1.88e–3	Periodic
pe82	0.8 × 0.8 × 0.2	256	×	256	×	64	223	3.78e–3	Periodic
pe84	0.8 × 0.8 × 0.4	256	×	256	×	128	223	5.12e–3	Periodic
pe88	0.8 × 0.8 × 0.8	256	×	256	×	256	223	1.36e–3	Periodic
ou22	0.2 × 0.2 × 0.2	64	×	64	×	64	219	1.78e–3	Outflow
ou24	0.2 × 0.2 × 0.4	64	×	64	×	128	219	2.02e–3	Outflow
ou28	0.2 × 0.2 × 0.8	64	×	64	×	256	219	1.70e–3	Outflow
ou42	0.4 × 0.4 × 0.2	128	×	128	×	64	221	1.65e–3	Outflow
ou44	0.4 × 0.4 × 0.4	128	×	128	×	128	221	2.19e–3	Outflow
ou48	0.4 × 0.4 × 0.8	128	×	128	×	256	221	1.47e–3	Outflow
ou82	0.8 × 0.8 × 0.2	256	×	256	×	64	223	1.55e–3	Outflow
ou84	0.8 × 0.8 × 0.4	256	×	256	×	128	223	1.69e–3	Outflow
ou88	0.8 × 0.8 × 0.8	256	×	256	×	256	223	1.19e–3	Outflow

Notes. For all runs: the solid-to-gas ratio Z = 0.02, the dimensionless particle stopping time τ_s = 0.314, the headwind speed ηv_K = 0.05c_s, the (cubic) grid cells are H/320 wide, and the run time is $314{{\rm{\Omega }}}_{0}^{-1}$ (=50P).

^aRun names consist of a two-letter abbreviation for the vertical BC followed by the first significant digit of the horizontal and then vertical size of the computational box. ^bThe number of particles. For reference, 2¹⁹ ≈ 5.2 × 10⁵, 2²¹ ≈ 2.1 × 10⁶, and 2²³ ≈ 8.4 × 10⁶. ^cThe time-averaged gas mass flux through vertical boundaries over the initial strong clumping phase. More specifically, for each snapshot, the horizontal average of the absolute value of mass flux $\langle | {\rho }_{{\rm{g}}}{c}_{{\rm{s}}}| \rangle$ over both the upper and lower boundaries was obtained.

Download table as:  ASCIITypeset image

The physical behavior of our simulations is controlled by three dimensionless parameters, which we hold fixed in order to focus on numerical issues. First, the dimensionless particle stopping time is τ_s ≡ Ω₀t_stop = 0.314. Physically, this choice corresponds to compact solid particles with sizes from 10 cm (at ∼1 au) to 1 mm (at ∼100 au; Youdin et al. 2010). Second, we fix the total solid-to-gas ratio at Z = 0.02, which sets the strength of drag feedback. The total gas mass used in Z includes the total vertical column density of gas, $\sqrt{2\pi }{\rho }_{{\rm{g}},0}H$, which extends beyond the computational domain. Particles remain close to the midplane and require no such correction. Third, the headwind parameter ηv_K/c_s = 0.05 sets the strength of global radial pressure gradients—which drive radial drift and streaming—to a typical value at several au.

Table 1 summarizes the numerical parameters for our simulations. For all of the runs, the initial vertical density profiles of gas and particles are Gaussians with scale heights of H and H_p = 0.025H, respectively, where H_p is defined as the rms of particle z_i coordinates. Both gas and particles are initialized with the Nakagawa–Sekiya–Hayashi (NSH) equilibrium drift velocities in Nakagawa et al. (1986).

Following the standard local shearing box approach, we apply the periodic BCs in both radial and azimuthal directions with azimuthal shear in the radial boundaries. For the vertical direction, in addition to previously used periodic and reflecting vBCs, we implement and explore the outflow vBCs. Under each of these three vBCs, we perform a series of nine runs with different numerical domain sizes, serving as a convergence test. Their horizontal sizes L_X = L_Y span from 0.2H to 0.8H systematically, while their vertical length L_Z spans the same range in the same way, but independently. Their names in Table 1 are composed of the abbreviation of the vBCs, the first significant digit of the box horizontal length, and the first significant digit of the box vertical length.

In order to minimize the influence of numerical resolution, we keep it fixed at 320 grid cells per H (cell size δl = 0.003125H). Bai & Stone (2010b) claimed that an equal number of particles and grid cells is sufficient for numerical convergence in the development of SI. Considering the disk stratification, the number of the particles should not depend on the height of the numerical box but the typical initial thickness of particle layer, 2H_p. Therefore, to reduce the number of free parameters in simulations, the number of particles is set by

$$\begin{eqnarray}&&{N}_{\mathrm{par}}={n}_{{\rm{p}}}\cdot {N}_{\mathrm{cell}}^{* }={n}_{{\rm{p}}}\times \displaystyle \frac{{L}_{X}{L}_{Y}\cdot 2{H}_{{\rm{p}}}}{{({\rm{\Delta }}x)}^{3}}=8388608\times \displaystyle \frac{{L}_{X}{L}_{Y}}{{(0.8H)}^{2}},\end{eqnarray} \tag{ 5 }$$

where n_p = 8 is the number of particles per cell and ${N}_{\mathrm{cell}}^{* }$ is the number of effective cells around the midplane between 2H_p. Thus, in the biggest run, there are over 8 million particles in a numerical box consisting of ∼16 million grid cells. Superparticle approximation is applied in our calculations. To be more specific with the ATHENA code, we use the semi-implicit particle integrator and a triangular shaped cloud (TSC) scheme to interpolate the particle properties to the grid cell for the need of computing interactions between gas and particles (Bai & Stone 2010b). Our simulations are computationally expensive due to high resolution and the huge amount of particles. Therefore, the ending point of each simulation is set to $314{{\rm{\Omega }}}_{0}^{-1}$, which is equivalent to 50P, where P is the orbital period and $P=2\pi {{\rm{\Omega }}}_{0}^{-1}$.

The SI is primarily of interest because it can concentrate particles to high densities and facilitate planetesimal formation. Figure 3 shows the maximum particle density, max(ρ_p)⁵ of all grid cells, as a function of time for all the runs.

Zoom In Zoom Out Reset image size

The evolution of max(ρ_p) is similar for most simulations, and can be divided into four stages. First, in the sedimentation phase, max(ρ_p) increases rapidly and steadily to ∼20ρ_g,0 at t ∼ 2P. In the second stage, the SI causes increasingly strong particle clustering, with max(ρ_p) reaching $\sim {10}^{3}{\rho }_{{\rm{g}},0}$ in most simulations by t ∼ 15P.

We define the third stage as lasting from t ∼ 15P to t ∼ 20P when peak particle densities saturate in all our simulations. The fourth stage, after t ∼ 20P, has similarly saturated values of max(ρ_p), with some fluctuations. We later show (in Section 3.3) that particle filaments undergo mergers during stage 4 (see also Yang & Johansen 2014). However, the particle densities in stage 3 are already high enough to trigger gravitational collapse in a self-gravitating simulation. We thus consider the clustering properties during stage 3 to be the most relevant, especially in terms of initial conditions for gravitational collapse calculations. We focus on the early strong clumping of stage 3 in most of our analysis.

Our most anomalous simulation is the case of the smallest box ((0.2H)³) with the reflecting vBCs (Run re22; see the blue curve in the upper left frame). This case shows significantly weaker clumping than all of the others, especially in stage 3. Aside from this anomaly, all of the simulations exhibit very strong particle clumping. One apparent trend is that the horizontally largest boxes (with L_X = L_Y = 0.8H) show more consistently strong clumping during stage 3. In the following subsection, we analyze clumping statistics in more detail.

Figure 4 compares the evolution of particle scale heights, H_p. In all simulations, H_p drops sharply to 0.001H in t ∼ 3P, as particles sediment. Then, H_p rises sharply to ∼0.005H, as the SI rapidly develops, with simultaneous particle clustering underway.

Zoom In Zoom Out Reset image size

After the initial sedimentation and rebound, the evolution of H_p varies significantly with box size and especially vBCs. With the outflow vBCs, the values of H_p remain below 0.0075H as indicated by the black dashed lines in Figure 4, whereas with the other two vBCs, H_p tends to increase and exceed 0.0075H. The reflecting vBCs cases evolve to the largest H_p values, especially in the shorter L_Z = 0.2H boxes. However, H_p is similarly thin during stage 3 in all tall boxes (L_Z = 0.8H), largely independent of vBCs or horizontal box size.

To better visualize what causes H_p variations in short boxes, Figure 5 presents the vertical extent of wave-shaped particle layers. The corrugated structures of these particle layers are well illustrated thanks to the axisymmetric particle clumping. From top to bottom, a larger H_p value results from a thicker particle layer as well as the larger amplitudes of the radial waves on it.

Zoom In Zoom Out Reset image size

Figure 5 also shows that the densest particle clumps (filaments viewed on-axis) lie very close to the midplane. The parts of the corrugated particle layer that are centered away off the midplane have lower particle densities.

The max(ρ_p) gives a useful but limited view of the overall particle clustering. While the peak particle density is useful for estimating whether any fragmentation should occur, it ignores the statistics of particle clustering on different length scales. Moreover, the results for max(ρ_p) depend on the cell size and thus do not numerically converge. For max(ρ_p), only the densest grid cell is considered, which disregards the rest of the domain and furthermore depends on numerical resolution. Thus, Figure 6 plots the maximum particle density on a range of length scales, time-averaged over t = 15–20P. This scale-dependent $\langle \max ({\rho }_{{\rm{p}}})\rangle ({\ell })$ gives the highest particle density within a sphere of radius ℓ located anywhere in the domain. Essentially all choices of vBC and box size give the same particle overdensity of ∼2ρ_g0 at ℓ = ηr, the characteristic radial length scale of disk pressure gradients and also of the linear (unstratified) SI for τ_s ≃ 1 and ρ_p ≲ ρ_g (see Figure 2 of Youdin & Johansen 2007). The only discrepant case is, again, the shortest box with reflecting vBCs. Particle overdensities increase toward smaller scales, down to the resolution limit of the simulations (as shown in Johansen et al. 2012).

Zoom In Zoom Out Reset image size

The upper panel of Figure 6 compares different box sizes for each of the vBCs. In all cases, larger boxes produce stronger particle overdensities, especially on small scales. The larger mass budget in boxes with a larger horizontal area offers the opportunity for larger density fluctuations, and the SI takes advantage. This horizontal area effect is also seen in Figure 3 and indicates that large boxes are needed to study the maximum extent of small-scale clustering.

The lower panel of Figure 6 presents the effect of vBCs on clumping by simply regrouping the contents in the upper panel. The clumping of particles across all scales is remarkably independent of the vBCs, again except for the short box reflecting case. This result is reassuring for the robustness of particle clustering by the SI since all imposed vBCs are artificial in some way.

Johansen et al. (2012) found that a power law of $\langle \max ({\rho }_{{\rm{p}}}){\rangle }_{t}\propto {{\ell }}^{-2}$ (indicated by the diagonal reference line) is a good approximation to the scale dependence of clumping. Such a scaling is expected for a linearly elongated structure, whose mass ∝ℓ averaged over a volume ∝ℓ³ gives a density ∝ℓ⁻². At small scales, this argument breaks down due to the finite width of elongated structures and the presence of overdense knots, as seen in Figure 1. Respectively, these two effects can give either a flatter or a steeper slope for $\langle \max ({\rho }_{{\rm{p}}})\rangle ({\ell })$. Thus, the fact that the slope stays as steep as −2, and even gets steeper in the larger boxes, is a sign of strong small-scale clumping. The existence of strong clumping on length scales ≪ηr is consistent with the rough expectations of linear SI theory. As μ ≡ ρ_p/ρ_g increases above unity, the radial length scale of the linear SI decreases sharply, while the growth rate increases (Youdin & Goodman 2005; Youdin & Johansen 2007; Squire & Hopkins 2018). These linear properties seem consistent with a cascade of stronger concentration to smaller length scales. Nevertheless, it is difficult to use the linear theory (which is also axisymmetric and unstratified) to predict this nonlinear clumping. Moreover, the box-size dependence of small-scale clustering shows that the small-scale clustering is not fully numerically converged, apparently because it is not a local process.

Figure 7 investigates whether the vertical extent of the computation domain affects particle clustering by comparing runs with the same horizontal area. Overall, the effect of L_Z on particle clustering is small, especially in the boxes with the largest horizontal area (L_X = L_Y = 0.8H). For all boundary conditions, shorter boxes do have somewhat lower particle overdensities on the smallest scales. (We again see that the smallest reflecting box (Run re22) is unique among our runs for its weak particle clustering.) Our finding that particle clustering depends only weakly on box height is reassuring since many self-gravitating simulations use short boxes (L_Z = 0.2H) to reduce computational costs.

Zoom In Zoom Out Reset image size

Figure 8 plots the probability distribution functions (PDFs) of the particle densities at the grid scale,⁶ again time-averaged over t = 15–20P. The left panel compares different box sizes for the outflow vBCs. Larger boxes produce higher peak densities, with PDFs that extend above 10³ρ_g0. The right panel shows that vBCs have an impressively negligible effect on density PDFs in the largest simulation box. The convergence properties of particle density PDFs on the grid scale thus agrees with the analysis of maximum density as a function of scale.

Zoom In Zoom Out Reset image size

To study the largest structures in our simulations, we analyze the azimuthally averaged surface density of both particles and gas. We first focus on the particles to better understand the spacing of azimuthal filaments. The spacing of these filaments affects the mass reservoir available to form planetesimals in the gravitational fragmentation phase (not modeled here). We then analyze axisymmetric gas structures to highlight two important points. First, the SI produces zonal flow structures, which, to our knowledge, have not been previously reported, and which depend strongly on horizontal box size and vBCs. Second, the pressure bumps associated with these zonal flows are too weak to trap particles, and thus are not directly related to the particle clumping mechanism of the SI.

3.3.1. Particle Ring Spacing

The left panels of Figures 9–11 show spacetime plots of the azimuthally averaged particle surface density, $\langle {{\rm{\Sigma }}}_{{\rm{p}}}{\rangle }_{y}(x)$, relative to Σ_g,0, the initial gas surface density. These three figures present outflow, periodic, and reflecting vBCs, respectively. The evolution is similar in all cases. Radial structures begin to develop at t ∼ 3P, first in many closely spaced filaments. Over time these filaments merge, and the simulations end with a small number of dense filaments, with only a single filament in some of the radially narrow boxes. The spacing of filaments clearly depends on time and on box size, as addressed below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, we address the radial drift. The oscillations in the radial location of the filaments occur on a timescale of one orbital period, which reflect epicyclic oscillations. Drift speeds vary between filaments, with denser filaments drifting slower. When a filament is disrupted, partially or totally, the escaping particles drift at higher speeds, approaching the value of a test particle. In the NSH equilibrium state, the radial drift speed of particles is

$$\begin{eqnarray}&&{v}_{x}=-\displaystyle \frac{2{\tau }_{{\rm{s}}}}{{(1+{\mu }^{2})}^{2}+{\tau }_{{\rm{s}}}^{2}}\eta {v}_{{\rm{K}}},\end{eqnarray} \tag{ 6 }$$

where μ ≡ ρ_p/ρ_g. For the case of a test particle (μ = 0),

$$\begin{eqnarray}&&{v}_{x,\mathrm{test}}=-\displaystyle \frac{2}{1+{\tau }_{{\rm{s}}}^{2}}\eta {v}_{{\rm{K}}}{\tau }_{{\rm{s}}}.\end{eqnarray} \tag{ 7 }$$

The black arrows in those three figures illustrate v_x,test, which are almost parallel with the trajectories of escaping particles. Equation (6) also provides a comparison between the radial drift speed of a filament and the equilibrium radial drift speed (represented by a purple arrow) with a certain μ value (labeled by purple text). Since there are so many substructures inside a filament, the averaged dust-to-gas density ratio, μ, is only of order 10¹, which corresponds to a width of ℓ ≈ 2 × 10 ⁻²H. The oscillations in the radial location of the filaments occur on a timescale of one orbital period, which reflect epicyclic oscillations.

At late times in our simulations, the spacing between (and the number of) dense particle filaments varies with vBC and box size, as seen by comparing Figures 9–11. It is dynamically interesting that vBCs can affect the growth of large-scale particle structures. However, the late-time behavior of our simulations is less relevant for planetesimal formation, as argued above (because if self-gravity were included then fragmentation would occur earlier). Thus, we analyze the clump spacing during the early strong clumping phase of t = 15–20P.

To show the radial spacing of particle filament, Figure 12 plots the Fourier power spectrum of $\langle {{\rm{\Sigma }}}_{{\rm{p}}}{\rangle }_{y}$. For boxes with larger L_X, the power spectrum extends to lower radial wavenumbers, k_X. In the narrowest boxes with L_X = 0.2H, the power is largest at the smallest k_X, consistent with the existence of a single dominant clump in these simulations (see Figures 9–11). Thus, L_X = 0.2H boxes are too narrow to determine the physical spacing between particle filaments.

Zoom In Zoom Out Reset image size

For wider boxes with L_X = 0.4H and 0.8H, the power has a maximum at intermediate wavenumbers of k_XH/(2π) = 7 ± 2 for L_X = 0.4H cases and k_X H/(2π) = 8 ± 2 for L_X = 0.8H cases. The corresponding wavelength, i.e., radial distance between peaks, is ∼3ηr, somewhat larger than the length scale, ηr, that is set by pressure gradients.

We thus conclude that boxes with L_X ≥ 6ηr should adequately capture the formation of multiple (i.e., at least two) particle filaments for any vBC. For our adopted value of ηr/H = 0.05, this width corresponds to L_X ≳ 0.3H. We emphasize that the spacing of filaments will also depend on τ_s, i.e., particle sizes, which we do not explore here.

3.3.2. SI-induced Zonal Flows

The right columns of Figures 9–11 plot perturbations to the azimuthally averaged gas surface density, ${\rm{\Delta }}\langle {{\rm{\Sigma }}}_{{\rm{g}}}{\rangle }_{y}$, for the three vBCs. The amplitude of the perturbations varies with box size and with vertical BC. The strongest gas density fluctuations are for the outflow conditions in the widest box, Run ou88. For this run, the amplitude of gas overdensities reaches 10⁻³ at late times (Figure 9). The periodic and reflecting cases (in Figures 10 and 11) only reach an amplitude of 1.5 × 10⁻⁴. Though we have previously argued that the late-time behavior of our simulations is not the most relevant, here we consider the late behavior for a couple reasons. First, the ability of the SI to generate these gas perturbations and (as described below) the corresponding zonal flows is not well known. Second, we want to emphasize that the largest amplitude gas bumps in our simulations are still not responsible for (and indeed not capable of) trapping particles by reversing the global pressure gradient.

The gas fluctuations typically contain only one radial wavelength per box, indicative of an inverse cascade to the largest scales. The inverse cascade is strongest and most persistent in the wide box with outflow vBCs (ou88). The wide periodic and reflecting cases (Runs pe88 and re88) also show an inverse cascade, but with more small-scale features and time-varying structures. The inverse cascade seen for the gas density is not shared by all aspects of the gas flow or particle dynamics. For instance, particle filaments and the gas vertical momentum contain significant small-scale structure (including axisymmetric structure), as seen in Figures 1, 2, and 5.

In a zonal flow, large-scale radial fluctuations in gas surface density give rise to pressure gradients that balance the Coriolis force of a radially varying azimuthal flow, i.e.,

$$\begin{eqnarray}&&2{\rho }_{{\rm{g}}}{{\rm{\Omega }}}_{0}\delta {u}_{y}\approx \displaystyle \frac{\partial P}{\partial x}={c}_{{\rm{s}}}^{2}\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial x},\end{eqnarray} \tag{ 8 }$$

where δu_y = u_y + 3/2Ω₀x is the deviation of the azimuthal flow from Keplerian. Geostrophically balanced zonal flows arise in simulations of the magnetorotational instability (MRI; Johansen et al. 2009a; Bai & Stone 2014). They tend to grow to very large radial scales, with a maximum scale of ∼6H seen in large box simulations with L_X up to 20 H (Simon et al. 2012; Dittrich et al. 2013). We do not measure the maximum radial extent of SI-induced zonal flows as this would require wider simulation domains.

Now we show the large-scale gas density perturbations induced by the SI are in fact zonal flows. Consider axisymmetric gas density fluctuations with a radial wavelength equal to the box width, with i.e., 2π/k₀ = L_X, as

$$\begin{eqnarray}&&\delta {\rho }_{{\rm{g}}}=\langle {\rho }_{{\rm{g}}}{\rangle }_{y}-\langle {\rho }_{{\rm{g}}}{\rangle }_{x,y}={\mathfrak{R}}({\hat{\rho }}_{{\rm{g}}}^{{\prime} }{e}^{{{ik}}_{0}x}),\end{eqnarray} \tag{ 9 }$$

where $\langle \cdot {\rangle }_{y}$ means an azimuthal average, $\langle \cdot {\rangle }_{x,y}$ denotes a horizontal average, and ${\hat{\rho }}_{{\rm{g}}}^{{\prime} }$ is a Fourier amplitude. Then Equation (8) gives a zonal flow speed

$$\begin{eqnarray}&&\delta {u}_{y}=\displaystyle \frac{{k}_{0}{c}_{{\rm{s}}}^{2}}{2{{\rm{\Omega }}}_{0}\langle {\rho }_{{\rm{g}}}{\rangle }_{{xy}}}{\mathfrak{R}}(i{\hat{\rho }}_{{\rm{g}}}^{{\prime} }{e}^{{{ik}}_{0}x}).\end{eqnarray} \tag{ 10 }$$

Figure 13 shows that gas density and azimuthal flow perturbations are indeed in geostrophic balance.

Zoom In Zoom Out Reset image size

Figure 13 also shows that the zonal flow is too weak to trap particles. The zonal flow would have to exceed the headwind speed of ηv_K = 0.05c_s to halt the inward drift of test particles. The zonal flow amplitude is only 0.01c_s in Figure 13, which shows the strongest zonal flow in all our simulations. Thus none of our simulations are close to creating a pressure trap for particles.

To further emphasize that the zonal flows and their associated pressure bumps appear unconnected with particle clustering. Note that (i) the particle clustering is essentially saturated when zonal flows are still growing, (ii) the radial spacing of particle filaments is initially much smaller than zonal flow widths, and (iii) zonal flows show a dramatic variation with box size and vBC, quite unlike particle clustering.

The fact that particle clumping appears unconnected from pressure traps is important for a better understanding of the mechanism of the SI. Some attempts to intuitively explain the SI appeal to the creation of zonal flows and pressure traps (e.g., Jacquet et al. 2011; Johansen et al. 2014). Our analysis shows that the full explanation is more complicated and that the dynamical gas motions are an essential part of the SI trapping mechanism.

Other works have started to explain this complexity. For instance, the (very closely related to the SI) secular drag layer instability of Goodman & Pindor (2000) is explained by a toy model of overstable oscillations (see also Chiang & Youdin 2010). Lin & Youdin (2017) show that many linear drag instabilities in disks (including the SI) are powered by "PdV" work caused by out-of-phase gas pressure and particle density perturbations. Thus, while pressure perturbations are essential to the SI (as they are for all hydrodynamics), any connection between quasi-static particle trapping in pressure bumps and the SI is not obvious.

Even though particle–gas interactions are confined to a thin midplane layer, these interactions trigger hydrodynamic activity that extends through the computational domain (see Figures 1 and 2). The vBCs have a direct and very strong affect on gas flow through the vertical boundaries. To quantify this behavior, Table 1 lists the time-averaged gas mass flux through vertical boundaries. Specifically, the table shows the horizontal average of the absolute value of mass flux ($\langle | {\rho }_{{\rm{g}}}{u}_{z}| \rangle$) over both the upper and lower boundaries. Thus, for the outflow vBCs, we get the outflow rate (through both boundaries) as ${\dot{{\rm{\Sigma }}}}_{{\rm{g}}}=2\langle | {\rho }_{{\rm{g}}}{u}_{z}| \rangle$. For the periodic vBCs, $\langle | {\rho }_{{\rm{g}}}{u}_{z}| \rangle$ measures the gas mass flux through the vertical boundary in either direction, from top to the bottom or vice versa. The direction of this periodic flow fluctuates in time. For the reflecting vBCs, there is no data since $\langle | {\rho }_{{\rm{g}}}{u}_{z}| \rangle =0$.

As Table 1 shows, the vertical fluxes of the outflow runs are all ∼10 ⁻³ρ_g,0c_s and vary little with box size. Each time step, a small fraction of gas mass (∼10⁻⁶) is lost before being replenished. Without this replenishment, the mass-loss timescale of ${10}^{3}{{\rm{\Omega }}}_{0}^{-1}$ is not only very short compared to observationally known disk lifetimes, but also physically implausible from the perspective of the energy budget. In order to overcome the gravitational potential of the star to escape globally, those gas flows need an energy flux of ${v}_{{\rm{K}}}^{2}{\dot{{\rm{\Sigma }}}}_{{\rm{g}}}$. However, the energy powering the SI only provides an energy flux of

$$\begin{eqnarray}&&\left[\displaystyle \frac{4{\tau }_{{\rm{s}}}}{{(1+\epsilon )}^{2}+{\tau }_{{\rm{s}}}^{2}}\right]{(\eta {v}_{{\rm{K}}})}^{2}{{\rm{\Sigma }}}_{{\rm{p}}}{{\rm{\Omega }}}_{0}\approx {(\eta {v}_{{\rm{K}}})}^{2}{{\rm{\Sigma }}}_{{\rm{p}}}{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 11 }$$

where [ · ] is order unity or smaller (see Equation (22) in Youdin & Johansen 2007). The SI is nowhere near sufficient to power outflows on a global scale (by a factor of ∼10⁻⁵).⁷

These vigorous outflow rates are thus unrealistic. In our solutions, the outflow is always subsonic, which explains why the rates are not physically realistic. Outflow rates are only reliable if all critical points (including MHD wave speeds when magnetic fields are included) are inside the computational domain so that the outflow rate is causally disconnected from the vBCs (Fromang et al. 2013). Unfortunately, an increased vertical extent does not allow us to find a physical outflow solution with a sonic point.⁸ The fundamental limitation appears to be the local geometry of the shearing box. Local simulations of MRI outflows find that the critical point for fast magnetosonic waves is not inside the computational domain (Bai & Stone 2013a, 2013b).

For the periodic vBCs, the flow of mass through the connected top and bottom boundaries does not attempt to represent physical outflow or inflow. This artificial connection between the vertical boundaries can introduce numerical artifacts. As can be seen from Table 1, the vertical flux in periodic cases are even stronger than that in the outflow solutions. Moreover, these periodic flows generally decrease with increasing L_Z, indicating that the vertical boundaries have fewer effects when placed farther away.

While reflecting vBCs have no flow through the vertical boundaries, reflection at these boundaries can still introduce numerical artifacts. Such effects could be much stronger in shorter boxes, resulting in the most anomalous case (Run re22) in the smallest box in our simulations.

We conclude that all vBCs introduce numerical artifacts that can affect gas motions throughout the simulation domain. We are fortunate that the clumping of larger (τ_s ≳ 0.1) particles in the midplane is not strongly affected by this imperfection. However, hydrodynamics away from the midplane is important for understanding the lofting of small dust. We thus caution that we cannot present physically robust results for dust lofting until either (i) different vBCs give convergent behavior or (ii) a physically robust outflow solution—with all sonic points—is found. Resolving this issue is an important topic for future work.