© ESO 2020

1 Introduction

Turbulence crucially influences various stages of the formation of planets: such as (1) the vertical settling of dust grains to a midplane layer whose thickness is determined by the equilibrium between sedimentation and turbulent diffusion (Dubrulle et al. 1995; Johansen & Klahr 2005; Fromang & Papaloizou 2006; Youdin & Lithwick 2007); (2) collisional grain growth (Ormel & Cuzzi 2007; Birnstiel et al. 2010); (3) planetesimal formation owing to passive concentration of grains (Barge & Sommeria 1995; Johansen et al. 2007; Cuzzi et al. 2008); and (4) planetary migration (Nelson & Papaloizou 2004; Oishi et al. 2007; Yang et al. 2009, 2012; Baruteau et al. 2011).

However, it is challenging to observationally constrain the strength of the turbulence in the gas and the dust in protoplanetary disks, whose motions are coupled via drag. This is particularly because these turbulent motions are weaker than both the thermal gas motions and the orbital motions of gas and dust (Flaherty et al. 2018).

Recently, high-resolution ALMA observations of protoplanetary disks permitted a number of authors to assess the turbulent strength in these disks: Flaherty et al. (2015, 2017, 2018) derive upper limits of the strength of the turbulent vertical gas motions in the disks surrounding HD 163296 and TW Hya from the nonthermal broadening of molecular emission lines. These upper limits correspond to Mach numbers of the order of 0.01¹.

Pinte et al. (2016) employ a model of micron- to millimeter-sized grains in the disk around HL Tau to estimate the strength of their vertical turbulent diffusion from the observed dust scale height, which is equal to ~10% of the gas scale height. These latter authors also find a Mach number of ~0.01. Similarly, Ohashi & Kataoka (2019) constrain the dust grain sizes and dust scale height in the disk around HD 163296 using polarization measurements of the dust emission. From a grain size of ~100 μm and a scale height of less than one-third of the gas scale height inside the ring that is located at an orbital radius of 70 au in this disk, these latter authors infer an upper limit of the Mach number of the vertical gas velocity of ~0.01. On the other hand, from the dust-to-gas scale height ratio of two-thirds outside of the ring the authors estimate a Mach number of the order of 0.1.

Furthermore, Dullemond et al. (2018) obtain a lower limit on the strength of the radial turbulent motions from the width of the dust rings that characterize the majority of the observed disks. This limit also amounts to a Mach number of about 0.01 for grains of 0.2 mm in size, but is proportional to the grain size.

In addition to the uncertainties in the observational determination of the turbulent strength, the physical processes that are thedominant sources of protoplanetary disk turbulence remain to be theoretically established. In large fractions of disks, particularly in the midplane, magnetohydrodynamic (MHD) turbulence is suppressed because of nonideal MHD effects (Gressel et al. 2015; Bai 2017). Turbulence in these regions must therefore be driven either from the MHD turbulent disk surface (Oishi & Mac Low 2009; Bai 2015) or by purely hydrodynamic instabilities.

One of the most promising among the hydrodynamical instabilities is the vertical shear instability (Arlt & Urpin 2004; Nelson et al. 2013), which is similar to the Goldreich-Schubert-Fricke instability in differentially rotating stars (Goldreich & Schubert 1967; Fricke 1968). It arises when the disk rotation profile depends on the height. This can, for instance, be due to baroclinity, which is a misalignment between the density and the pressure gradient, resulting from a radial temperature profile. The source of energy of the instability is the free energy associated with the vertical shear (Barker & Latter 2015).

The vertical shear instability can overcome both the radial angular momentum gradient, as its modes are characterized by a large radial-to-vertical wavenumber ratio, and the vertical buoyancy if the gas cooling timescale is sufficiently short (Nelson et al. 2013; Lin & Youdin 2015). Analytical and numerical analyses of the linear growth of the instability have found two classes of modes: short-wavelength surface modes with higher growth rates and vertically global body modes with lower growth rates. The former appear at the artificial vertical simulation domain boundary, where the vertical shear is strongest – although natural transitions in the density can also give rise to these modes. Their growth rate increases with the vertical shear at theboundary, and thus with vertical domain size (Nelson et al. 2013; Barker & Latter 2015; Lin & Youdin 2015).

In nonlinear simulations, the vertical shear instability grows over at least approximately 30 orbital periods until it attains a saturated state (Stoll & Kley 2014; Flock et al. 2017). The Mach numbers of the turbulent vertical motions in this state are of the order of 10⁻² to 10⁻¹ (Flock et al. 2017). This turbulent strength is higher if the radial temperature gradient is steeper (Nelson et al. 2013; Lin 2019). Perturbations associated with the surface modes appear first and grow towards the disk midplane. The later emerging body modesare characterized by perturbations that evolve from an odd symmetry with respect to the midplane into an even symmetry, and therefore the instability saturates last in this plane (Nelson et al. 2013; Stoll & Kley 2014).

The turbulence that is induced by the vertical shear instability entails angular momentum being transported radially outwards and vertically away from the midplane. Since the latter eliminates the vertical shear, external heating (of locally nonisothermal disks) is necessary to sustain the instability (Stoll & Kley 2014). In numerical models including dust and the drag exerted on it by the gas, the turbulence further gives rise to vertical dust motions, and radial motions leading to accumulation in overdensities of up to five times the initial dust density (Stoll & Kley 2016; Flock et al. 2017).

However, Lin & Youdin (2017) and Lin (2019) show that, if the drag back-reaction of the dust onto the gas is taken into account as well, the dust, which sediments to the midplane, introduces an effective vertical buoyancy – it “weighs down” the gas – that can quench the vertical shear instability. If they are tightly coupled, gas and dust can be described as a single fluid, with a density equal to the sum of the gas and the dust density, but a pressure that is solely due to the gas. We assume that the pure gas is locally isothermal, that is, its cooling timescale is infinitely short. Under this assumption, the equation of state of the mixture is given by $P=c_{\text{s}}^2{\rho }_{\text{g}}=c_{\text{s}}^2(1-{\rho }_{\text{d}}/{\rho }_{\text{tot}}){\rho }_{\text{tot}}$, where P is the gas pressure, c_s the sound speed, and ρ_g, ρ_d, and ρ_tot are the gas, dust, and total density, respectively. In other words, the gas-dust mixture is not locally isothermal, and its cooling timescale is finite. As noted above, the instability is suppressed by vertical buoyancy if the cooling timescale is too long.

The streaming instability (Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007) results from the inwards radial drift of dust, which is caused by the difference in orbital velocity between the gas and the dust as well as their mutual coupling via drag. In contrast to the dust, the gas orbits with a sub-Keplerian velocity because it is supported against the radial stellar gravity by a pressure gradient. This pressure gradient constitutes a source of free energy that is tapped by the streaming instability (Youdin & Johansen 2007). The linear growth rate of the instability is highest if the dust drift is fastest, that is when the dust stopping time is comparable to the inverse of the orbital frequency (Weidenschilling 1977), and if the dust-to-gas density ratio is slightly greater than one (Youdin & Goodman 2005). Physical interpretations of the instability have been devised by Lin & Youdin (2017) and Squire & Hopkins (2018).

The turbulence that the streaming instability gives rise to in its nonlinear regime can result in dust concentration in axisymmetric filaments. In these filaments, the dust accumulates in clumps that are sufficiently dense to collapse under their self-gravity and form planetesimals (Johansen et al. 2007; Bai & Stone 2010b; Yang & Johansen 2014; Simon et al. 2016; Schäfer et al. 2017). However, whether the streaming instability causes dust accumulation that is strong enough for planetesimal formation depends on the dust-to-gas surface density ratio and the size of the dust (Johansen et al. 2009; Bai & Stone 2010b; Carrera et al. 2015; Yang et al. 2017, 2018), as well as the strength of the radial pressure gradient (Bai & Stone 2010c). The required surface density ratio is higher than 1% – the canonical value in the interstellar medium – for all dust sizes. It can be enhanced sufficiently to reach the critical value globally by photoevaporation (Carrera et al. 2017; Ercolano et al. 2017), and locally in radial dust pile-ups (Drążkowska et al. 2016)or at ice lines (Ida & Guillot 2016; Schoonenberg & Ormel 2017; Schoonenberg et al. 2018; Drążkowska & Alibert 2017).

Observations of binaries among the Trans-Neptunian objects provide evidence for planetesimal formation owing to the streaming instability: their observed orbital inclinations relative to their heliocentric orbit as well as the ones that are found in simulations of binary formation by gravitational collapse are predominantly prograde (Grundy et al. 2019; Nesvorný et al. 2010, 2019). In contrast, dynamical binary capture leads to either mostly retrograde inclinations or a similar number of prograde and retrograde ones, depending on the capture mechanism (Schlichting & Sari 2008).

While itis considered to be and is studied as one of the most promising mechanisms to induce planetesimal formation, the streaming instability has received comparably little attention as a process to drive turbulence, although this is its fundamental effect. Moreover, the streaming instability has been studied numerically almost exclusively in local shearing box simulations. Only Kowalik et al. (2013) present global two- and three-dimensional simulations, which reproduce the dust accumulation in dense axisymmetric filaments.

In this paper, we study the streaming instability as a source of turbulence, employing axisymmetric global simulations with considerablylarger domains than the ones that were simulated by Kowalik et al. (2013). In contrast to these authors, we model the dust as Lagrangian particles rather than as a pressureless fluid, and take into account the vertical stellar gravity, which leads to the dust sedimenting to the midplane. We further apply adaptive mesh refinement to enhance the resolution of the dust midplane layer.

The paper is structured as follows: in Sect. 2, the simulations, their initial conditions, and the parameters that govern their evolution are described. In Sects. 3–5, respectively, we present our study of the turbulence in models of the vertical shear instability only, both the vertical shear and the streaming instability, and the streaming instability only. Implications and limitations of the study are discussed in Sect. 6. We summarize our results in Sect. 7.

2 Numerical model

We perform two-dimensional numerical simulations of the gas and the dust components of protoplanetary disks, including the mutual drag between the two components as well as the radial and the vertical stellar gravity. Magnetic fields, the self-gravity of the gas and the dust, and (nonnumerical) viscosity are neglected. We employ version 4.5 of the adaptive-mesh refinement (AMR) finite volume code FLASH Code² (Fryxell et al. 2000).

2.1 Domains, boundary conditions, and resolutions

The cylindrical, axisymmetric simulation domains extend from either 1 to 10 au or from 10 to 100 au in the radial dimension and 1 or 2 gas scale heights above and below the midplane in the vertical dimension. Since in our model the gas scale height increases nonlinearly with the radius (see Eq. (9)), the domains are shaped like isosceles trapezoids with curved rather than straight legs. To model this shape as accurately as possible with respect to the initial resolution, we initially create rectangular domains with a vertical size of two or four gas scale heights at the outer radial boundary. From these domains, we then exclude all blocks of 10 × 10 grid cells whose distance to the midplane is greater than one or two local gas scale heights. For this purpose, we apply the obstacle block implementation that is included in the FLASH Code³.

Diode conditions are applied at both the radial and the vertical boundaries, that is, the boundaries are permeable to gas and dust moving out of the domain, but reflect gas and dust that would move into it. At the vertical boundaries, the pressure in the guard cells is quadratically interpolated to maintain hydrostatic equilibrium. The temperature in the guard cells at all boundaries is reset to the initial value because we find this to be conducive to the stability of our simulations. In addition, the orbital velocity is corrected to account for the difference in temperature and stellar gravity between the guard and the nonguard cells³.

The FLASH Code employs the PARAMESH package (MacNeice et al. 2000) to perform block-structured AMR. In other words, the domain is subdivided into blocks of 10 × 10 grid cells, which as a whole are refined or derefined if a refinement or derefinement criterion is fulfilled in any cell inside them. To resolve the gas and dust dynamics that are induced by the streaming instability in the disk midplane, we apply a (de-) refinement criterion which is based on the spatial dust distribution: the resolution is doubled if more than ten particles, which we use to model the dust, are located in one cell. On the other hand, the resolution is halved if no particles remain in a cell³. In addition, we initially increase the resolution of the midplane density at the inner radial domain boundary by a factor of two or of four where the ratio of the gas density to the midplane density at the inner radial domain boundary exceeds 1 or 10%, respectively.

In the domains with a radial size of 9 au, the fiducial initial and maximum resolution amount to 160 and 5120 cells per astronomical unit, respectively, while in the domains with an extent of 90 au they are equal to 10 and 320 cells per astronomical unit, respectively. At the maximum resolution, this corresponds to more than 200 cells per gasscale height at all radii, which has been shown to be sufficient to resolve the formation of gravitationally unstable dust clumps owing to the streaming instability in local shearing box simulations (Yang & Johansen 2014). To investigate whether or not our findings are dependent on resolution, we also conduct simulations with a doubled initial and maximum resolution.

The simulation names and parameters are compiled in Table 1. Sets of simulations which are analyzed in different sections are separated by a double horizontal line. All names are composed of (1) the gas equation of state (isothermal or adiabatic); (2) the dust-to-gas surface density ratio (this is omitted if dust is not included in a simulation); and (3) the radial simulation domain size. Where applicable, the names further indicate that (4) the vertical domain extent amounts to four gas scale heights; (5) the resolution is twice the fiducial resolution; (6) the dust particles are introduced after 50 kyr rather than at the beginning of the simulation; (7) the initial dust-to-gas scale height ratio is set to 1 or 100%, with the fiducial value being 10%; and (8) the dust particle size deviates from the fiducial size of a = 3 cm.

2.2 Gas

The equations of motion of the gas can be expressed as $$\frac{\partial {\rho }_{\text{g}}}{\partial t}+\nabla \cdot \left({\rho }_{\text{g}}{v}_{\text{g}}\right)=0\text{ and}$$(1) $$\frac{\partial v_{\text{g}}}{\partial t}+\left(v_{\text{g}}\cdot \nabla \right)v_{\text{g}}=-\frac{1}{{\rho }_{\text{g}}}\nabla P-\nabla {\Phi }_{\text{S}}+\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}\frac{v_{\text{d}}-v_{\text{g}}}{t_{\text{stop}}},$$(2)

where v is the velocity, and the subscripts g and d refer to the gas and the dust, respectively. The stellar gravitational potential is given by ${\Phi }_{\text{S}}=-G{M}_{\text{S}}/\sqrt{r^2+z^2}$, where G is the gravitational constant, M_S = 1 M_⊙ is the stellar mass, and r and z, respectively,are the radial distance to the star and the height above or below the disk midplane. The last term on theright-hand side of Eq. (2) results from the drag exerted by the dust on the gas, with t_stop being the stopping time (see Sect. 2.3).

To closethe system of equations, we consider either an isothermal or an adiabatic equation of state. In the former case, the pressure is calculated as $$P=\frac{RT}{\mu }{\rho }_{\text{g}},$$(3)

where R is the ideal gas constant, T the temperature, and μ = 2.33 the mean molecular weight³. In this case, the adiabatic index γ = 1. To model the latter case, we employ a polytropic equation of state that is given by $$P=K{\rho }_{\text{g}}^{\gamma },$$(4)

where $K=RT{\rho }_{\text{g}}^{1-\gamma }/\mu$ is the polytropic constant and the adiabatic index γ = 5∕3³. While it is a local constant in time, the polytropic constant depends on the global temperature and density distributions (see Eqs. (5)–(7)). If the gas is locally adiabatic, the vertical shear instability is stabilized by vertical buoyancy.

The initial temperature is adopted from the minimum mass solar nebula (MMSN) model (Hayashi 1981), $$T=280 \text{K} {\left(\frac{r}{1 \text{au}}\right)}^{-1/2}.$$(5)

The steepness of this profile is in agreement with that of power-law fits to observed temperature distributions (Andrews & Williams 2005). The radial temperature gradient gives rise to a variation of the orbital speed with height. This vertical shear in turn is the source of energy of the vertical shear instability.

The gas initially orbits with a sub-Keplerian velocity, which is determined by the balance between stellar gravity, centrifugal force, and pressure gradient. It is furthermore in vertical hydrostatic equilibrium. As it is vertically isothermal, the vertical density profile is thus given by $${\rho }_{\text{g}}={\rho }_{\text{g}}(z=0)\exp \left[-\frac{\gamma G{M}_{\text{S}}}{c_{\text{s}}^2}\left(\frac{1}{r}-\frac{1}{\sqrt{r^2+z^2}}\right)\right],$$(6)

where $c_{\text{s}}={(\gamma RT/\mu )}^{1/2}\propto r^{-1/4}$ (see Eq. (5)) is the sound speed. We set the initial midplane density to $${\rho }_{\text{g}}(z=0)={10}^{-9} \text{g}\hspace{0.17em}{\text{cm}}^{-3} {\left(\frac{r}{1 \text{au}}\right)}^{-9/4}.$$(7)

Numerically integrating over Eq. (6) yields the surface density $${\Sigma }_{\text{g}}={\displaystyle {\int }_{-1 H_{\text{g}}}^{1 H_{\text{g}}}{\rho }_{\text{g}}} \text{d}z={10}^3 \text{g}\hspace{0.17em}{\text{cm}}^{-2} {\left(\frac{r}{1 \text{au}}\right)}^{-1},$$(8)

where $$H_{\text{g}}=\sqrt{\frac{c_{\text{s}}^2{r}^3(2\gamma G{M}_{\text{S}}-c_{\text{s}}^{2}r)}{{(c_{\text{s}}^{2}r-\gamma G{M}_{\text{S}})}^2}},$$(9)

is the gas scale height. This surface density profile is shallower than the one in the MMSN model (Hayashi 1981), but is consistent with that of observed young, massive protoplanetary disks (Andrews et al. 2009, 2010).

We notethat Eq. (8) gives the surface density as the density integrated over two gas scale heights, not integrated from −∞ to ∞ as in the commonly used definition. Consequently, the total mass M_g,tot in each of our domains depends not only on the radial, but also on the vertical domain extent: $$M_{\text{g},\text{tot}}=\{\begin{array}{ll}7.1\times {10}^{-3} M_{\odot } \left(\frac{L_r}{9 \text{au}}\right)\hfill & L_z\approx 2 H_{\text{g}}\text{ or}\hfill \\ 8.0\times {10}^{-3} M_{\odot } \left(\frac{L_r}{9 \text{au}}\right)\hfill & L_z\approx 4 H_{\text{g}},\hfill \end{array}$$(10)

where L_r and L_z are the radial and vertical domain size, respectively.

The orbital speed can be expressed as v_g,ϕ = v_K − Πc_s, where $$v_{\text{K}}=\sqrt{\frac{G{M}_{\text{S}}}{{\left(r^2+z^2\right)}^{3/2}}{r}^2}$$(11)

is the Keplerian speed. We adopt the dimensionless parameter Π, which is introduced by Bai & Stone (2010b) to indicate the strength of the radial pressure gradient $$\Pi =-\frac{1}{2c_{\text{s}}{\rho }_{\text{g}}{\Omega }_{\text{K}}}\frac{\text{d}P}{\text{d}r},$$(12)

where Ω_K = v_K∕r is the Keplerian orbital frequency. In the midplane, the parameter is equal to $$\Pi (z=0)=4.6\times {10}^{-2}\frac{1}{\sqrt{\gamma }}{\left(\frac{r}{1 \text{au}}\right)}^{1/4};$$(13)

(see Eqs. (5) and (7)).

2.3 Dust

We model the dust as Lagrangian particles using the massive active particle implementation that is included in the FLASH Code. We adopt an approach that is commonly used in local shearing box studies of the streaming instability (Youdin & Johansen 2007; Bai & Stone 2010a): the mass and momentum of every simulated particle are equal to the total mass and momentum of a large number of dust pebbles, while the drag coupling to the gas is the same as that of a single pebble.

The dust is initially uniformly distributed in the radial dimension. The mass of the dust particles is determined by their total number N_d = 10⁶, the dust-to-gas surface density ratio Z, and the gas surface density. It can be expressed as $$m_{\text{d}}=\frac{1}{N_{\text{d}}}{\displaystyle {\int }_{L_r}2}\pi r{\Sigma }_{\text{d}} \text{d}r=\frac{Z}{N_{\text{d}}}{\displaystyle {\int }_{L_r}2}\pi r{\Sigma }_{\text{g}} \text{d}r.$$(14)

Because the gas surface density is inversely proportional to the radius (see Eq. (8)), the mass of all particles in a simulation is given by $$\begin{array}{ll}{m}_{\text{d}}\hfill & =1.27\times {10}^{23} \text{g} \left(\frac{Z}{0.01}\right)\left(\frac{L_r}{9 \text{au}}\right)\hfill \\ \hfill & =0.14 M_{\text{Ceres}} \left(\frac{Z}{0.01}\right)\left(\frac{L_r}{9 \text{au}}\right),\hfill \end{array}$$(15)

where M_Ceres = 9.3 × 10²³ g is the mass of Ceres. The simulated dust-to-gas surface density ratios range from Z = 0.01 to Z = 0.1, with Z = 0.02 being the fiducial value. We verified that our results are converged with respect to the total number of particles by comparing simulations with N_d = 5 × 10⁵ and N_d = 10⁶.

The initial vertical positions are randomly sampled from a Gaussian distribution with a scale height of 10% of the gas scale height. This scale height is in agreement with the thickness of the dust midplane layer observed by Pinte et al. (2016). To investigate the dependence of our findings on the initial scale height, we additionally perform simulations with a dust-to-gas scale height ratio of 0.01 and of 1. The noise inthe vertical distribution serves as a seed for the streaming instability. A comparison of simulations with two different vertical distributions has shown that our results do not noticeably depend on the random seed.

We simulate dust with a fixed size of a = 3 cm or of a = 3 mm, or with a fixed dimensionless stopping time, which is equivalent to the Stokes number, of τ_stop = 0.1. Even if a = 3 cm, the dust is smaller than the gas mean free path length at all densities in our model. Under this condition, that is to say in the Epstein regime, the dimensionless stopping time in the midplane is given by $$\begin{array}{lll}{\tau }_{\text{stop}}(z=0)\hfill & =\hfill & t_{\text{stop}}(z=0){\Omega }_{\text{K}}(z=0)=\frac{a{\rho }_{\text{s}}}{c_{\text{s}}{\rho }_{\text{g}}(z=0)}{\Omega }_{\text{K}}(z=0)\hfill \\ \hfill & =\hfill & 6\times {10}^{-4} \frac{1}{\sqrt{\gamma }}\left(\frac{a}{3 \text{mm}}\right)\left(\frac{r}{1 \text{au}}\right),\hfill \end{array}$$(16)

where ρ_s = 1 g cm⁻³ is the dust material density (see also Eq. (7)).

We note that the collisional growth to sizes greater than millimeters is prevented by dust grains bouncing or fragmenting under mutual collisions (Güttler et al. 2010; Zsom et al. 2010; Birnstiel et al. 2011), except for in the innermost regions of protoplanetary disks (Birnstiel et al. 2012) and at ice lines (Ros & Johansen 2013; Ros et al. 2019). However, simulating centimeter-sized grains allows us to probe dimensionless stopping times of the order of τ_stop = 0.1. These are pertinent to a study of the turbulence that is driven by the streaming instability since the linear growth rate of the instability is highest if the dimensionless stopping time is close to one (Youdin & Goodman 2005).

The dust particles initially orbit with the Keplerian speed (see Eq. (11)). They are initialized either at the beginning of the simulations or after 50 kyr. The latter isto give the vertical shear instability time to attain a saturated state before the introduction of the dust. To advance the particlesin time, we use the Leapfrog algorithm, which was adapted for cylindrical geometries by Boris (1970), as detailed in Appendix A³.

The implementation of the mutual drag between the gas and the dust is based on the cloud-in-cell mapping between the grid and the particles that is included in the FLASH Code. The algorithm can be described as follows: firstly, the gas properties are mapped to the particles. Secondly, for each particle, the stopping time (see Eq. (16)) and the change in velocity dueto the drag exerted by the gas is computed. The corresponding change in particle momentum is then mapped to the grid. Finally, the change of the gas velocity in every grid cell is calculated from the particle momentum change. A more detailed description of the implementation can be found in Appendix B³.

3 Vertical shear instability

In this section, we verify that our model can reproduce the findings of previous studies of the vertical shear instability, in particular the turbulent strength in its saturated state, despite the comparably small vertical domain sizes we consider. For this purpose, we analyze our model of locally isothermal gas in which dust is not included.

In Fig. 1, the vertical gas motions in the model are depicted. The characteristic perturbations that the vertical shear instability gives rise to are reproduced. The perturbations are bent outwards, their radial wavelength is much less than their vertical wavelength, and they are symmetric with respect to the midplane in the saturated state of the instability (compare with, e.g., Figs. 2 and 3 of Nelson et al. 2013 and Fig. 2 of Stoll & Kley 2014).

We find that the strength of the turbulence induced by the vertical shear instability depends on the extent of the vertical domain. This can be seen from Fig. 2, in which we depict the vertical gas velocity in our models with a vertical domain size of two and of four gas scale heights. Throughout the simulations, the turbulence is considerably weaker in the vertically smaller domain than in the vertically larger one. The reason for this is likely that the vertical shear at the vertical domain boundaries is less in smaller domains, which entails a decline of the linear growth rate of the surface modes of the instability (Lin & Youdin 2015).

In the vertically larger domain, the vertical shear instability saturates at a Mach number of the vertical gas velocity of ${\mathcal{M}}_{\text{g},z}\approx 0.1$. This value does not depend significantly on the resolution and is consistent with the Mach number which Flock et al. (2017) find. This is regardless of the fact that these authors simulate a domain with an aspect ratio of z∕r = 0.35, while our domain extends to between z∕r = 0.17 at r = 10 au and z∕r = 0.3 at r = 100 au. Furthermore, Flock et al. (2017) employ a radiative transfer model rather than assuming the gas to be locally isothermal as we do. Therefore, in the following we investigate the vertical shear instability and its coexistence with the streaming instability using our model with a vertical domain size of two gas scale heights above and below the midplane.

Fig. 1 Mach number of the vertical gas velocity ${\mathcal{M}}_{\text{g},z}$ as a functionof radius r and height z in a simulation including a locally isothermal gas, but no particles (iso_Lr=90au_Lz=4Hg). The simulation domain spans four gas scale heights in the vertical dimension. Black lines mark one gas scale height above and below the midplane. After 30 kyr, the vertical shear instability has saturated at all radii and heights; it gives rise to perturbations which are characterized by a radial-to-vertical wavelength ratio much smaller than unity, by being bent outwards, and by a mirror symmetry relative to the midplane.

Fig. 2 Root mean square of ${\mathcal{M}}_{\text{g},z}$ as a functionof time t. The root mean square is computed using the mass-weighted average over the vertical domain size and over 1 au extending from r = 50  to 51 au. In the domain with a vertical dimension of four gas scale heights (blue and orange lines), the Mach number saturates at ${\mathcal{M}}_{\text{g},z,\text{rms}}\approx {10}^{-1}$. This value is consistent with the one which Flock et al. (2017) find in their numerical study of the vertical shear instability. The turbulent strength is similar in the simulation with the fiducial resolution (iso_Lr=90au_Lz=4Hg; blue line) and the one with the doubled resolution (iso_Lr=90au_Lz=4Hg_dou_res; orange line); however, it remains lower until the simulations end in the domain with a vertical size of two gas scale heights (simulation iso_Lr=90au; green line).

4 Coexistence of vertical shear instability and streaming instability

Using our model with dust and an isothermal gas equation of state, we investigate how the vertical shear instability and the streaming instability interact. We consider two scenarios: in the scenario that we refer to as SIwhileVSI,the streaming instability grows simultaneously with the vertical shear instability. In the scenario referred to here as SIafterVSI, on the other hand, the streaming instability is not active before the vertical shear instability has saturated. We model the latter scenario by introducing the dust into the simulations after 50 kyr. At this point, the vertical shear instability has reached a saturated state at all radii (see Fig. 1).

4.1 Source of turbulence

Figure 3 shows the Mach number of the vertical gas motions in SIwhileVSI (left panel) and SIafterVSI (right panel). In both scenarios, away from the midplane the vertical shear instability is the primary source of turbulence and induces the characteristic large-scale perturbations (compare with Fig. 1). However, in the midplane, small-scale perturbations can be seen at all radii in the former scenario. Similar perturbations exist in the latter scenario, although only at certain radii and not as limited in vertical extent. These small-scale perturbations are not present in our model of the vertical shear instability only. They resemble the perturbations that Li et al. (2018) observe in the midplane of their local shearing box simulations of the streaming instability (see their Fig. 2).

In Fig. 4, we compare the time-dependence of the Mach number in the two scenarios. The mass-weighted vertical average of the Mach number is depicted, that is, the turbulent strength in the midplane is weighted more heavily than the strength away from it. In SIafterVSI, we find that the initialization of the dust does not cause the Mach number to deviate significantly from the value of ${\mathcal{M}}_{\text{g},z}\approx 0.1$ induced by the vertical shear instability after it has saturated (compare with Fig. 2).

In contrast, in the SIwhileVSI scenario the Mach number saturates at a lower value of ${\mathcal{M}}_{\text{g},z}\approx 0.01$. In the local shearing box simulations without vertical stellar gravity presented by Johansen & Youdin (2007), the streaming instability drives turbulence with a similar strength. After saturation, the Mach number does not vary significantly until, owing to its radial drift, locally no dust is left. Subsequently, it increases to approximately the value the vertical shear instability in its saturated state gives rise to in SIafterVSI.

Figure 5 shows the ratio of the radial to the vertical Mach number in both scenarios. The turbulent velocity in the radial direction is significantly less than the one in the vertical direction in SIafterVSI. This is consistent with Stoll & Kley (2016) and Stoll et al. (2017) showing that the vertical shear instability drives anisotropic turbulence. On the other hand, in SIwhileVSI the radial and the vertical Mach number are comparable near and in the midplane. Johansen & Youdin (2007) find the streaming instability to cause isotropic turbulence.

From the fact that the turbulence is comparably weak and isotropic as well as from the presence of small-scale perturbations that are notobservable in our model of the vertical shear instability only, we conclude that the streaming instability is the primary source of turbulence in the dust midplane layer if it starts to operate at the same time as the vertical shear instability. One possible explanation for this is that the streaming instability grows faster in turbulent strength than the vertical shear instability. This is evident when comparing the time taken for the vertical shear instability (see Fig. 2) and the streaming instability (see the right panel of Fig. 4) to saturate.

However, the turbulence is mainly driven by the vertical shear instability if it has attained a saturated state before the streaming instability can begin to grow. In this state, we find it to give rise to turbulence with a higher Mach number than the streaming instability in SIwhileVSI. Nevertheless, we note that the small-scale perturbations, which we find to be induced by the streaming instability at all radii in SIwhileVSI, are also observable locally in the SIafterVSI scenario (see the right panel of Fig. 3).

The vertical shear instability would likely be the main source of turbulence in SIwhileVSI as well if the dust grains were too small to settle and trigger the streaming instability before the vertical shear instability had saturated. Equating a saturation timescale of 30 orbital periods (Stoll & Kley 2014) with the settling timescale as given by Eq. (10) of Chiang & Youdin (2010) yields a critical dimensionless dust stopping time of ~0.005. However, since simulating such small dust grains as particles is computationally very expensive, we do not explore this scenario.

Fig. 3 Mach number ${\mathcal{M}}_{\text{g},z}$ as a functionof r and z in two simulations with dust and an isothermal gas equation of state. As in the dust-free equivalent (compare with Fig. 1), the vertical shear instability induces large-scale perturbations away from the midplane. In the simulation iso_Z=0.02_Lr=90au_Lz=4Hg (left panel), in which we initialize the dust at the beginning, small-scale perturbations can be seen at all radii in the midplane. These perturbations are also observable in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr (right panel), in which the dust is introduced after 50 kyr, although not at all radii and extending to greater heights. The inlays with a vertical extent of 0.6 gas scale heights show the enlarged perturbations. However, such small-scale perturbations do not exist in the corresponding simulation without dust.

Fig. 4 Root mean square of ${\mathcal{M}}_{\text{g},z}$, averaged over the vertical domain extent and the radii given in the legend, as a function of t after the initialization of the dust at td,init. The mean is weighted by mass, i.e., more weight is assigned to the Mach number in the midplane than to the one away from it. In the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr (left panel), the vertical shear instability has saturated before the dust is introduced. The Mach number remains approximately constant at a value ${\mathcal{M}}_{\text{g},z,\text{rms}}\approx {10}^{-1}$ after the dust initialization (compare with Fig. 2). That is, in this simulation the vertical shear instability is the primary source of turbulence in the midplane. In contrast, in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg (right panel), in which the vertical shear instability and the streaming instability begin to grow at the same time, the Mach number saturates at a lower value of ${\mathcal{M}}_{\text{g},z,\text{rms}}\approx {10}^{-2}$. After no dust remains at a radius because of its radial drift (this point is marked with a dashed line), the Mach number increases until it is approximately equal to the value at the same radius in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr. We find that, as long as dust is present, the streaming instability drives the turbulence in the midplane in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg.

Fig. 5 Ratio of the root mean square Mach number of the radial gas velocity ${\mathcal{M}}_{\text{g},r}$ to the root mean square of ${\mathcal{M}}_{\text{g},z}$ as a functionof z in units of gas scale heights. The abscissa is scaled linearly between − 0.1 and 0.1 gas scale heights and logarithmically otherwise. To compute the root mean square, we take the mass-weighted average over r = 50  to 51 au and the average over 1 kyr through t = 55 kyr in thecase of the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr (blue line) and through t = 10 kyr in thecase of the simulation iso_Z=0.02_Lr=90au_Lz=4Hg (orange line). Away from the mid-plane, the radial-to-vertical Mach number ratio is less than unity in both simulations. While this anisotropy extends over all heights if the dust is introduced after 50 kyr, the Mach number ratio is close to one in the midplane if the dust is initialized at the start. This is indicative of the vertical shear instability being the primary source of turbulence in the midplane in the former case, but the streaming instability in the latter case.

4.2 Dependence of turbulent strength on dust density

As explained in Sect. 1, settling dust introduces an effective vertical gas buoyancy that can suppress the vertical shear instability. This is a second possible reason – besides the more rapid growth in turbulent strength of the streaming instability – why the streaming instability is the source of turbulence in the dust layer in the scenario SIwhileVSI.

Lin (2019) performs simulations of the vertical shear instability which include dust from the beginning, as in the SIwhileVSI scenario. This latter author shows that the dust-induced buoyancy leads to a decrease in the Mach number in the midplane from ${\mathcal{M}}_{\text{g},z}\approx {10}^{-1}$ to ${\mathcal{M}}_{\text{g},z}\approx {10}^{-2}$ if the dust-to-gas surface density or the stopping time of the dust exceed a threshold value. These threshold values are correlated: for a surface density ratio of 1%, the critical dimensionless stopping time of the dust is equal to 0.005 (see their Fig. 5). For a stopping time of 0.001, the critical surface density ratio amounts to 3% (see their Fig. 9). This is consistent with the Mach number of ${\mathcal{M}}_{\text{g},z}\approx {10}^{-2}$ we find in SIwhileVSI, in which the surface density ratio is equal to 2%⁴ and the stopping time of the dust ranges from 0.046 to 0.46.

However, while Lin (2019) associates the turbulence in the midplane of their simulations with the vertical shear instability, we find it tobe driven by the streaming instability if both instabilities begin to operate at the same time. Lin (2019) notes that they apply a diffusive numerical scheme which probably suppresses the growth of the streaming instability. In addition, they model tightly coupled dust and gas employing a one-fluid approach.

To investigate whether the vertical shear instability is quenched by the dust-induced buoyancy in the SIafterVSI scenario, we analyze howthe turbulent strength depends on the dust density. We note that this analysis does not involve information about the spatial distribution of the dust. In other words, the strength at a certain dust density is the sum of the bulk motions of regions with this density and the internal turbulence in these regions.

In Fig. 6, the Mach number is shown as a function of the dust-to-gas density ratio in both scenarios. Overall, the Mach number is higher in SIafterVSI (orange, green, and purple lines) than in SIwhileVSI (blue line). This is consistent with stronger turbulence in the dust layer being induced by the vertical shear instability in the former scenario, and weaker turbulence by the streaming instability in the latter scenario.

In the SIafterVSI scenario, at low to intermediate dust-to-gas volume density ratios the Mach number decreases with increasing dust-to-gas surface density ratio, although only slightly. This can be explained by the vertical shear instability being suppressed if the surface density ratio, that is, the total dust mass, is greater, which leads to the vertical bulk motions of regions with these volume density ratios being weaker.

Furthermore, the vertical shear instability is gradually quenched if the dust density increases to values greater than the gas density. The Mach number is aslow as ${\mathcal{M}}_{\text{g},z}\approx 0.01$ for the highest density ratios. This turbulent strength is similar to that caused by the streaming instability throughout the dust layer in the SIwhileVSI scenario. However, we note that even in the regions with the highest density ratios in SIafterVSI, the Mach number is greater than in the regions with the same density ratios in SIwhileVSI. This is because the vertical shear instability gives rise to stronger bulk motions of these regions in the former scenario than the streaming instability in the latter scenario.

The turbulent strength that the streaming instability induces in the scenario SIwhileVSI is also lower if the density ratio exceeds a few. This is in agreement with Johansen et al. (2009) who find that the collision speeds of dust grains in the filaments forming in their simulations of the streaming instability are smaller if the dust density is higher.

Fig. 6 Root mean square of ${\mathcal{M}}_{\text{g},z}$ as a functionof the dust-to-gas density ratio ρd∕ρg. The orange, green, and purple lines represent the turbulent strength in simulations in which we initialize the dust after the vertical shear instability has saturated. The strength at low to intermediate values of the volume density ratio in the three simulations is marginally smaller if the dust-to-gas surface density ratio is greater. Nevertheless, the Mach number is higher in all of these simulations than in the simulation in which we introduce the dust at the start (blue line). We find that the stronger turbulence in the former simulations is driven by the vertical shear instability, while the weaker turbulence in the latter one is caused by the streaming instability. Independent of the source of the turbulence, its strength decreases at dust densities larger than the gas density in all four simulations.

Fig. 7 Ratio of the dust scale height to the gas scale height Hg (left panel) and maximum of the dust-to-gas density ratio ρd∕ρg (right panel)as functions of t after the dust is introduced at td,init. We compute the dust scale height as the root mean square of the vertical particle positions zd, and average the scale height ratio over 1 au spanning from r = 20  to 21 au. For a dust-to-gas surface density ratio of Z = 2%, the dust scale height amounts to ~1% of the gas scale height if the streaming instability induces the vertical diffusion of the dust (blue line). On the other hand, it is equal to ~10% for the same surface density ratio if the diffusion is caused by the vertical shear instability (orange line). The scale height which is induced by the vertical shear instability decreases with increasing surface density ratio (green and purple line). However, it is higher than the value that the streaming instability gives rise to for all surface density ratios that we consider. In contrast, the maximum dust volume density is significantly smaller if the streaming instability drives the turbulence in the dust layer than if the turbulence is caused by the vertical shear instability. In the former case, the dust-to-gas volume density ratio amounts to a few hundred for a surface density ratio of 2%. In the latter case, on the other hand, it exceeds 103 for the same surface density ratio, and 104 for a surface density ratio of 10%. In all cases, the maximum dust density is greater than the Roche density ρR in the midplane at the inner radial boundary of the simulation domains, i.e., the maximum Roche density in the simulations. The ratio of this Roche density to the gas density is marked as a black line.

4.3 Vertical and radial dust concentration

The dust-to-gas scale height ratio in SIwhileVSI (blue line) and SIafterVSI (orange, green, and purple lines) is shown in the left panel of Fig. 7. In the former scenario, the streaming instability is the source of turbulence in the dust layer. The dust sedimentation and the vertical diffusion of the dust which it causes reach a balance at a dust-to-gas scale height ratio of ~1% (see also Sect. 5).

Since the vertical shear instability drives stronger turbulence than the streaming instability, we find the equilibrium dust scale height to be greater in the SIafterVSI scenario. On the other hand, it decreases with increasing dust-to-gas surface density ratio in this latter scenario: for a ratio of 2%, the dust scale height amounts to ~10% of the gas scale height. In comparison, if the ratio is equal to 10% the scale height is close to the value in SIwhileVSI. We show in Sect. 5 that the scale height which is induced by the streaming instability is largely independent of the surface density ratio.

The dust settling to smaller scale heights for higher surface density ratios is most probably a consequence of the vertical shear instability being more suppressed by the dust-induced buoyancy. Lin (2019) finds the instability to diffuse dust with a dimensionless stopping time of 0.001 to a scale height that is similar to the gas scale height if the surface density ratio amounts to 1%, but that the dust scale height is an order of magnitude smaller if the ratio is equal to 5% (see their Fig. 9). These scale heights are significantly greater than the ones in the SIafterVSI scenario. This is likely because we simulate dust with stopping times of between 0.046 and 0.46, which is more weakly coupled to the gas and is therefore less elevated by the vertical gas motions.

To investigate whether or not in the SIafterVSI scenario the vertical shear instability can be quenched if the dust is initialized with a smaller scale height and thus a higher midplane density, we conducted a simulation with an initial surface density ratio of 2% and an initial dust scale height of 1% of the gas scale height rather than the fiducial value of 10%. Despite thisscale height being comparable to the one which is induced by the streaming instability in the SIwhileVSI scenario and the initial midplane dust-to-gas density ratio being of order unity, we find that the vertical shear instability is not noticeably affected. The dust is elevated to a dust-to-scale height ratio of ~10% in less than an orbital period in this simulation.

In the right panel of Fig. 7, we depict the maximum dust-to-gas density ratio in the two scenarios. We note that this maximum is stochastic and dependent on the resolution (Johansen & Youdin 2007; Bai & Stone 2010c). Nonetheless, we find the maximum dust density to exceed the maximum Roche density,

Fig. 8 Dust-to-gas volume density ratio ρd∕ρg as a functionof r and z (upper panels) as well as dust-to-gas surface density ratio Σd∕Σg as functions of r and t (lower panels). The upper panels show the spatial dust distribution 5 kyr after the dust initialization. We compare a simulation in which the turbulence in the dust layer is driven by the streaming instability (left panels) and two simulations in which the vertical shear instability is the main source of turbulence (middle and right panels). For the same dust-to-gas surface density ratio of Z =0.02, the dust scale height is smaller, but the radial dust concentration is weaker in the former case. (The values of Z are specified in the titles.) In the latter case, more dust is accumulated in overdensities if the surface density ratio is higher. Some of the accumulations are sufficiently dense for their radial drift to cease almost entirely.

$${\rho }_{\text{R},\text{max}}=\frac{9{\Omega }_{\text{K}}(r=10 \text{au}, z=0)}{4\pi G},$$(17)

in both scenarios. In other words, if the dust self-gravity were included in our model, local dust overdensities could undergo gravitational collapse and form planetesimals. This is in line with expectations for the SIwhileVSI scenario since the surface density ratio of 2% in this scenario exceeds the critical value for the dust concentration by the streaming instability to be strong enough to lead to planetesimal formation (Carrera et al. 2015; Yang et al. 2017).

If the radial concentration of the dust were comparable in both scenarios, the maximum dust density would by tendency be greater in SIwhileVSI because the dust scale height is smaller in this scenario. On the contrary, we find the maximum dust density to be considerably higher in SIafterVSI. The maximum of the density ratio amounts to a few hundred in SIwhileVSI, but to a few thousand for the same surface density ratio in SIafterVSI, and to more than 10⁴ for higher surface density ratios. We investigate the radial dust concentration in the following.

In the upper panels of Fig. 8, we show the spatial dust distribution 5 kyr after the dust is introduced in SIwhileVSI (left panel) and in SIafterVSI (middle and right panels). In agreement with what can be seen from Fig. 7, the oscillating dust midplane layer extends to greater heights in the latter scenario than in the former. In addition, the dust is more uniformly distributed in the radial dimension in the former scenario. Interestingly, the region between r = 15  and 30 au is depleted of dust in both simulations of the SIafterVSI scenario that are presented in Fig. 8. The reason for this depletion is unclear.

The time- and radius-dependence of the dust-to-gas surface density ratio in SIwhileVSI is depicted in the lower-left panel of Fig. 8. The streaming instability causes the dust to radially accumulate in overdensities. Comparable, azimuthally elongated filaments are found in three-dimensional simulations of the instability (Johansen et al. 2007; Bai & Stone 2010b; Kowalik et al. 2013). As noted above, the dust concentration in this scenario is sufficiently strong for the streaming instability to induce planetesimal formation.

From the lower-middle and lower-right panels, it can be seen that similar, but significantly denser dust concentrations form in SIafterVSI. If the surface density ratio is higher, that is to say the total dust mass is greater, more dust is contained inoverdensities. The radial drift of overdensities is reduced in both scenarios. However, it is only in the SIafterVSI scenario that some of the accumulations are dense enough for their radial drift to be nearly completely halted.

Stoll & Kley (2016) perform simulations of the vertical shear instability in which dust and the drag exerted on it by the gas are included. These latter authors find that pressure fluctuations which are induced by the vertical shear instability lead to a radial concentration of the dust. This concentration is strongest if the dimensionless stopping time of the dust is close to unity. Nevertheless, because they do not take the drag back-reaction of the dust on the gas into account, the radial drift speed does not depend on the dust density in their simulations.

In the vicinity of the dust overdensities in the SIafterVSI scenario, we find the streaming instability to contribute to the driving of turbulence. At the end of the simulation of the SIafterVSI scenario with a surface density ratio of 2% – 55 kyr after the beginning of the simulation and 5 kyr after the initialization of the dust – several overdensities can be found between r = 30 au and 40 au, and around 60 au (see the lower-middle panel of Fig. 8). From the right panel of Fig. 3, it is evident that at this time and these radii small-scale perturbations are present. We associate these perturbations with the streaming instability because they are not observable in our model of the vertical shear instability only. Therefore, it is most probably a combination of the vertical shear instability and the streaming instability that induces the dust accumulation in this scenario. We speculate that the vertical shear instability first concentrates the dust in weak overdensities, as shown by Stoll & Kley (2016). Subsequently, the streaming instability causes the growth of these seeds to the dense accumulations that we find in our simulations.

The fact that the dust is radially more concentrated in SIafterVSI despite its vertical diffusion being stronger in this scenario is consistent with the results of the numerical study by Yang et al. (2018). These latter authors show that the radial dust concentration owing to the streaming instability is not significantly affected by the vertical dust diffusion that is induced by nonideal MHD turbulence.

Fig. 9 Mach number of the vertical gas velocity as a function of r and z in simulations of dust and a locally adiabatic gas. The Mach number at radii of less than ~ 3 au is similarlyhigh in the simulation without (adi_Lr=9au; left panel) or with dust (adi_Z=0.02_Lr=9au; right panel). This is in spite of the streaming instability operating in the latter simulation, but not in the former one. At larger radii in the simulation that includes dust, the streaming instability induces large-scale perturbations with ${\mathcal{M}}_{\text{g},z}\approx {10}^{-3}$ away from the midplane. These perturbations are similar in shape to the ones caused by the vertical shear instability, but they are weaker, are bent inwards rather than outwards, and are not symmetric with respect to the midplane (compare with Fig. 1). In the inlay, which extends to 0.1 gas scale heights Hg above and below the midplane, small-scale perturbations with ${\mathcal{M}}_{\text{g},z}\approx {10}^{-2}$ in the midplane can be seen.

5 Global simulations of the streaming instability

To study the vertical and radial gas motions induced by the streaming instability, as well as the dust scale height, we employ our model including dust and locally adiabatic gas. Because the gas is locally adiabatic, the vertical shear instability is quenched by vertical buoyancy.

The Mach number of the vertical motions in this model as well as its dust-free equivalent are shown in Fig. 9. Since in the latter model (left panel) the streaming instability is not active, the turbulence in this model is most probably a numerical artifact caused by vertical shear that results from imperfect boundary conditions. At radii greater than ~3 au, this artificial turbulence is negligibly weak compared to the turbulence driven by the streaming instability in the model with dust (right panel). However, at smaller radii, we cannot distinguish between artificial turbulence and turbulence caused by the streaming instability. We therefore exclude these radii from the following analysis.

In the model including dust, the streaming instability induces small-scale perturbations with a turbulent strength of ${\mathcal{M}}_{\text{g},z}\approx 0.01$ in the dust midplane layer (see the inlay in the right panel) and somewhat weaker large-scale perturbations away from it. The small-scale perturbations are similar to the ones caused by streaming instability in the model in which both this latter instability and the vertical shear instability operate (see Fig. 3). We find both kinds of perturbations to be largely isotropic (see also Fig. 5).

The large-scale perturbations resemble the perturbations caused by the vertical shear instability in that their radial-to-vertical wavelength ratio is much less than one (see Fig. 1). In contrast to these, their symmetry is odd with respect to the midplane and they are bent inwards. This bending of the perturbations can be explained by gas moving both vertically away from the midplane and radially outwards, the latter being caused by the inward radial drift of the dust and the conservation of angular momentum.

In the local shearing box simulations presented by Li et al. (2018), the streaming instability gives rise to the same two kinds of perturbations with a similar turbulent strength (compare with their Fig. 2). However, the large-scale perturbations observed by these latter authors are bent outwards rather than inwards. Nevertheless, these perturbations are likely suppressed by the strongerones caused by the vertical shear instability under more realistic conditions, that is, in simulations including heating and cooling by radiation rather than an adiabatic equation of state (Stoll & Kley 2014; Flock et al. 2017).

5.1 Vertical gas velocity and dust scale height

In Fig. 10, we show the dust-to-gas scale height ratio as well as the Mach number of the vertical gas velocity in the midplane after dust settling and vertical diffusion have reached an equilibrium. We find both the dust scale height and the turbulent strength to be similar in the model with a dust-to-gas surface density ratio of 2% and the model with a ratio of 1%. The equilibrium scale height and the Mach number do not depend significantly on the initial dust scale height either, regardless of the increase in time taken for the dust to sediment to the midplane if the scale height is greater. Furthermore, our results are converged with respect to resolution.

Bai & Stone (2010b) study local shearing box simulations of the streaming instability with surface density ratios of 1, 2, or 3%. These latter authors derive the strength of the vertical diffusion as well as the scale height of the dust from fitting a Gaussian distribution to the vertical dust density profile. In contrast to us, they find the vertical diffusion to be stronger and the scale height to increase with the dust-to-gas surface density ratio if the ratio is less than a threshold value. Above this threshold value, which depends on the particle size, the turbulent strength decreases again, and the dust settles to smaller scale heights.

The dust scale height is equal to ~1% of the gas scale height at all radii in our model. The streaming instability gives rise to a similar dust scale height in our model in which both this latter instability and the vertical shear instability are active and begin to grow at the same time (see the left panel of Fig. 7). Comparable dust scale heights are further found in local shearing box simulations, for example, the ones that were conducted by Yang & Johansen (2014) and Carrera et al. (2015). For these scale heights and the dust-to-gas surface density ratios of the order of 1% which we simulate, the ratio of the dust to the gas density in the midplane is of order unity, and the linear growth rate of the instability is largest (Youdin & Goodman 2005).

The thickness of the dust layer is set in a self-regulatory manner if the turbulent diffusion of the dust is caused by the streaming instability (Bai & Stone 2010b): if the dust scale height is greater than the equilibrium value, the instability induces weaker turbulence, which does not balance the dust settling towards the midplane. On the other hand, if the scale height is less than the equilibrium value, the instability drives overly strong turbulence, which lifts the dust away from the midplane.

At radii r ≳10 au, the Mach number is consistent with the observed values of the order of 10⁻² (Flaherty et al. 2015, 2017, 2018; Pinte et al. 2016; Ohashi & Kataoka 2019). The turbulent strength is comparable in our model in which both the vertical shear instability and the streaming instability operate, but the latter drives the turbulence in the dust layer. Pinte et al. (2016) and Ohashi & Kataoka (2019) derive the turbulent strength from the dust scale height, which is regulated by the streaming instability if it is the main source of turbulence in the dust layer of protoplanetary disks. We note that these latter authors observe a dust-to-gas scale height ratio of 10% and one-third, respectively, which is an order of magnitude larger than the ratio we find. This can be explained by the fact that these latter authors consider micron- and millimeter-sized dust grains, which are elevated to greater heights than the centimeter-sized ones in the simulations that we present in Fig. 10.

Fig. 10 Ratio of the dust scale height, as the root mean square of zd, to the gas scale height Hg (left panel) as well as mass-weighted root mean square of ${\mathcal{M}}_{\text{g},z}$ in the midplane (right panel). Both quantities are depicted as functions of r. To compute the root mean square, we average over 50  or 500 yr after the dust scale height has reached an equilibrium value in the simulation domains spanning 1 au ≤ r ≤ 10 au or 10 au ≤ r ≤ 100 au, respectively. The dashed line marks the boundary between these domains. The model with an initial dust-to-gas scale height ratio of ${(H_{\text{d}}/H_{\text{g}})}_{\text{init}}=0.1$, a dust-to-gas surface density ratio of Z = 2%, and the fiducial resolution (blue line) is shown together with the models that deviate from this fiducial one in that ${(H_{\text{d}}/H_{\text{g}})}_{\text{init}}=1$, that Z = 1% (green line),or that the initial and maximum resolution are doubled (purple line). While both the dust scale height and the Mach number depend on the radius, they are largely independent of the initial dust scale height, the surface density ratio, and the resolution.

5.2 Dependence of turbulent strength on dust stopping time and radial gas pressure gradient

We find that the strength of the turbulence that is driven by the streaming instability, like the linear growth rate of the instability (Youdin & Goodman 2005), increases with the speed of the radial dust drift. The drift is faster if the dimensionless stopping time of the dust is closer to one and the radial gas pressure gradient is stronger. Both quantities increase with the radius in our model of dust with a fixed size (see Eqs. (13) and (16)). From Fig. 10, it is evident that both the Mach number and the dust scale height are also consequently greater at larger radii.

To analyze the extent to which the turbulent strength depends on the stopping time and the magnitude of the pressure gradient individually, we compare three simulations with different dimensionless stopping times: in one, it is fixed at 0.1, while in the other two it ranges from 0.0046 to 0.046 and from 0.046 to 0.46, respectively. The Mach number in each of these simulations is depicted in Fig. 11 as a function of the radius and the pressuregradient strength. We find an equivalent radius-dependence of the dust scale height in the three simulations.

The turbulent strength increases with both the dimensionless pressure gradient parameter and the dimensionless stopping time until it saturates for respective values of Π ≈ 0.1 and of τ_stop ≈ 0.05. This is consistent with Bai & Stone (2010a,b) finding the scale height of the dust and its vertical diffusion to increase with the strength of the pressure gradient and the stopping time. The dependence on the pressure gradient strength is evident when considering only the simulation with a fixed stopping time (blue line in the figure).

At all radii in the simulation with a dimensionless stopping time of 0.1 and the simulation with stopping times greater than 0.046 (green line),the turbulent strength is about the same, and therefore it is largely independent of the stopping time if τ_stop ≳ 0.05. In contrast, in the simulation with stopping times of less than 0.046 (orange line) the turbulence is overall weaker, but its strength increases more strongly with the radius than in the other two simulations. That is, if τ_stop ≲ 0.05 the strength depends not only on the pressure gradient magnitude, but also on the stopping time.

Fig. 11 Rootmean square of ${\mathcal{M}}_{\text{g},z}$ in the midplane as a function of r (lower abscissa) and the dimensionless gas pressure gradient parameter Π in the midplane (upper abscissa). The root mean square is calculated using the mass-weighted mean over 500 yr after an equilibrium value of the dust scale height is reached. In the simulation adi_Z=0.02_Lr=90au_taustop=0.1 (blue line), in which the dimensionless dust stopping time τstop is fixed at 0.1, the Mach number increases with the pressure gradient parameter if Π ≲ 0.09, but is approximately constant if the pressure gradient is stronger. The Mach number in this one and the simulation adi_Z=0.02_Lr=90au (green line), in which the stopping time increases from 0.046 to 0.46, is comparable at all radii. That is, the turbulent strength does not depend significantly on the stopping time in this stopping time regime. The Mach number is lower overall, but depends more strongly on the radius in the simulation adi_Z=0.02_Lr=90au_a=3mm (orange line) in which the stopping time ranges from 0.0046 to 0.046. In this regime, the turbulent strength increases with both the strength of the pressure gradient and the stopping time.

5.3 Radial gas velocity

In Fig. 12, we show the radial velocity of the gas. The gas moves outward on average as a consequence of the inward radial drift of the dust and, as the total angular momentum is conserved, the transfer of angular momentum from the dust to the gas. Yet, the turbulent motions that are induced by the streaming instability lead to the variance in the gas velocity being greater thanthe speed of this mean outward motion. That is, the streaming instability does not cause a sustained inward transport of gas mass (and outward transport of angular momentum) that contributes to the observed stellar accretion (Alcalá et al. 2017).

6 Discussion

6.1 Implications for turbulence and planetesimal formation in protoplanetary disks

We have studied the interaction of two instabilities, the vertical shear instability and the streaming instability. Both instabilities appear to be robust: operation of the vertical shear instability requires a vertical rotation profile and a sufficiently short gas cooling timescale. Lin & Youdin (2015) show that in the MMSN model, the instability can grow between r ≈ 5  and 50 au. On the other hand, only dust and a radial gas pressure gradient are necessary for the streaming instability to be active.

We considered two scenarios: in the first scenario, the vertical shear instability and the streaming instability start to grow at the same time. In this scenario, the turbulence in the dust layer is driven by the streaming instability. If the vertical shear instability has already saturated beforehand, as is the case in the second scenario, it remains the main source of turbulence even while the streaming instability is active.

It is unclear which of these scenarios is more realistic: Monte-Carlo simulations of dust coagulation show that (sub-)micron-sized grains grow to millimeter- or centimeter-sized aggregates within 10³ to 10⁴ orbital periods (Zsom et al. 2010; Lorek et al. 2018). However, the point at which dust growth commences during the formation of a protoplanetary disk is unclear, as is the point at which the conditions in a disk become conducive to the growth of the vertical shear instability and the streaming instability.

The vertical shear instability gives rise to stronger turbulence, and thus a stronger vertical dust diffusion, than the streaming instability. Nonetheless, we find this instability, in combination with the streaming instability, to cause the dust to be substantially more concentrated in the radial dimension than the streaming instability alone. We note that previous studies show the dust accumulation owing to the streaming instability to be sufficient to lead to planetesimal formation for the dust-to-gas surface density ratios we consider.

In other words, for a given dust size, the critical surface density ratio that is required for planetesimal formation may be significantly lower if not only the streaming instability, but both the vertical shear instability and the streaming instability together induce the dust concentration. The critical surface density ratio exceeds 1% for all dust sizes in the former case (Carrera et al. 2015; Yang et al. 2017), but might be lower than this canonical interstellar medium value in the latter case.

Fig. 12 Radial gas velocity (solid line) as a function of r in the simulations adi_Z=0.02_Lr=9au and adi_Z=0.02_Lr=90au. (The dashed lines marks the boundary between the domains of the two simulations.) The velocity is computed as the mass-weighted average over the vertical domain extent and a time-span of 50  and 500 yr, respectively,after an equilibrium dust scale height has been attained in each simulation. The average radial velocity at r ≥ 3 au is plotted as a dotted line and amounts to 0.035 m s−1. This outward motion is caused by the dust drifting radially inwards and its angular momentum being transferred to the gas. However, it is evident that the mean of the velocity is less than its standard deviation, which is equal to 0.16 m s−1. This is a result of the turbulent motions caused by the streaming instability.

6.2 Limitations of our numerical study

Both the linear vertical shear instability (Nelson et al. 2013; Barker & Latter 2015) and the linear streaming instability (Youdin & Goodman 2005) are axisymmetric in nature. Nonetheless, slight deviations from this symmetry are found in simulations of the nonlinear regime of both the former instability (Nelson et al. 2013; Stoll & Kley 2014) and the latter one (Kowalik et al. 2013). These deviations are not captured by our two-dimensional simulations.

Umurhan et al. (2019) and Chen & Lin (2020) analytically examined the linear streaming instability in connection with the α-model for protoplanetary disk turbulence (Shakura & Sunyaev 1973). These latter authors find that turbulence reduces the growth rate of the instability compared to the purely laminar case. However, this analysis is not applicable if the streaming instability itself is the dominant source of turbulence.

Similarly, Gole et al. (2020) study planetesimal formation owing to the streaming instability in local shearing box simulations with driven Kolmogorov-like turbulence that is not affected by the presence of the dust. In the simulations of these latter authors, planetesimals formation is hampered by turbulence and quenched if the turbulent Mach number is greater than ~10⁻².

We studied only models with a single dust size (or a single dimensionless dust stopping time). Krapp et al. (2019) show that the linear growth rate of the streaming instability can decrease if multiple dust species with a distribution of sizes rather than a single species are simulated. Nevertheless, simulations of the nonlinear streaming instability including dust size distributions do not seem to show indications of this effect (Bai & Stone 2010b; Schaffer et al. 2018).

Furthermore, we did not take into account other hydrodynamic instabilities like the convective overstability (Klahr & Hubbard 2014; Lyra 2014), the subcritical baroclinic instability (Klahr & Bodenheimer 2003; Klahr 2004; Lyra 2014), or the zombie vortex instability (Marcus et al. 2015, 2016; Lesur & Latter 2016). These instabilities do not operate in our model because the former two require finite gas cooling timescales – we simulate either a locally adiabatic or a locally isothermal gas – while the latter is three-dimensional.

In addition, our simulations do not include magnetic fields. Cui & Bai (2020) study the vertical shear instability in nonideal MHD simulations. These latter authors find that magnetized disk winds and the vertical shear instability can coexist. Under conditions that are typical of protoplanetary disks, the instability induces a comparable turbulent strength in their simulations and in purely hydrodynamical ones like ours. Nonetheless, if the magnetization is enhanced or the gas is more strongly coupled to the magnetic field, the instability causes weaker turbulence.

7 Summary

We present two-dimensional axisymmetric global numerical simulations of protoplanetary disks spanning orbital radii between 1  and 100 au. The simulations include Lagrangian particles to model the dust, the mutual drag between dust and gas, and the radial and vertical stellar gravity. We used the FLASH Code to conduct these simulations, which allowed us to apply adaptive mesh refinement to increase the resolution locally in and close to the dust layer in the midplane of the disks.

Employing these simulations, we investigated the turbulence driven by the vertical shear instability and that driven by the streaming instability individually, as well as the interaction of the two instabilities. The results of our study can be summarized as follows:

• We conducted simulations of the vertical shear instability only, with vertical domain extents of one or two gas scale heights above and below the midplane. Only the latter vertical size is sufficient to reproduce the turbulent strength found by previous numerical studies of vertically larger domains. The Mach number of the vertical gas motions is of the order of 10⁻¹ in the saturated state of the instability (Flock et al. 2017).

• If both the vertical shear instability and the streaming instability start to grow simultaneously, we find the turbulence in the dust midplane layer to be mainly driven by the streaming instability. This is most likely the result of a combination of two effects: in the midplane, the streaming instability grows faster in turbulent strength than the vertical shear instability. Furthermore, the weight of the dust induces an effective buoyancy in the gas that quenches the vertical shear instability (Lin & Youdin 2017; Lin 2019).

• The vertical dust settling and the turbulent diffusion that is induced by the streaming instability attain an equilibrium if the dust-to-gas scale height ratio is equal to ~1%. The dust scale height is set in a self-regulatory way if the streaming instability gives rise to the diffusion of the dust (Bai & Stone 2010b): if the scale height is less than the equilibrium value, then the turbulent strength is greater than the equilibrium strength, and the dust is lifted away from the midplane.

• We show that the streaming instability drives isotropic turbulence with a Mach number of up to ~10⁻². This is in agreement with observed values in protoplanetary disks (Flaherty et al. 2015, 2017, 2018; Pinte et al. 2016; Ohashi & Kataoka 2019). In particular, Pinte et al. (2016) and Ohashi & Kataoka (2019) obtain this Mach number from the scale height of the dust disks surrounding HL Tau and HD 163296, respectively; that is, they would probe the turbulent strength induced by the streaming instability if this instability were the primary source of turbulence in the dust layer of protoplanetary disks.

• Both the equilibrium dust scale height and the Mach number that are induced by the streaming instability are largely independent of the dust-to-gas surface density ratio and the initial dust scale height. The turbulent strength, and with it the scale height, increases with the speed of the radial dust drift. In other words, the turbulence is stronger if the dust stopping time or the radial gas pressure gradient is greater. The strength saturates for dimensionless stopping times of ~0.05 and dimensionless pressure gradient parameters, as defined by Bai & Stone (2010b), of ~0.1.

• In contrast, if the vertical shear instability has attained a saturated state before we introduce the dust into our simulations, then this instability remains the primary source of turbulence in the dust layer; it gives rise to stronger turbulence than the streaming instability, which elevates the dust to greater scale heights. For a dust-to-gas surface density ratio of 2%, the instability induces a Mach number of ~10⁻¹ and a dust scale height of ~10% of the gas scale height. Nevertheless, if the surface density ratio is higher, the instability is more strongly quenched by the dust-induced buoyancy.

• We find that a combination of the vertical shear instability and the streaming instability leads to a considerably strongerradial concentration of the dust than the streaming instability only. The dust accumulations are dense enough for their radial drift to be halted almost completely. This is despite the vertical shear instability inducing stronger vertical diffusion than the streaming instability. We speculate that the vertical shear instability induces the formation of weak overdensities that seed the streaming instability. The streaming instability in turn causes strong dust concentration that would likely lead to planetesimal formation in simulations including the self-gravity of the dust.

Acknowledgements

We thank the anonymous referee and Chao-Chin Yang for their constructive feedback that helped to improve thispaper. To analyze and visualize the simulations, the Python packages yt (http://yt-project.org) (Turk et al. 2011), Matplotlib (https://matplotlib.org) (Hunter 2007), and NumPy (https://numpy.org) (Oliphant 2006) have been used. The FLASH Code has in part been developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. Computational resources employed to conduct the simulations presented in this paper were provided by the Regionales Rechenzentrum at the University of Hamburg, by the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN), and by the Swedish Infrastructure for Computing (SNIC) at LUNARC at Lund University. US thanks the University of Hamburg for granting him a scholarship to fund his doctoral studies. A.J. is thankful for research support by the European Research Council (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg Academy Fellow Grant 2017.0287), and the Swedish Research Council (Project Grant 2018-04867). US and R.B. gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG), grant BA 3706/18-1. R.B. is thankful for support by the Excellence Cluster 2121 “Quantum Universe” which is funded by the DFG.

Appendix A: Leapfrog algorithm for cylindrical geometries

A.1 Implementation

We implement a second-order-accurate, explicit Leapfrog algorithm for the time integration of particles in cylindrical geometries. The algorithm is based on the one developed by Boris (1970) for charged particles in simulations including electric and magnetic fields. It has been described as well by, for example, Delzanno & Camporeale (2013).

For the update of the vertical components of the particle velocity and position, we adopt the Leapfrog algorithm for Cartesian geometries that is part of the FLASH Code: $$v_{z,i}^{1/2}=v_{z,i}^0+\frac{1}{2}{a}_{z,i}^0\Delta t^0,$$(A.1) $$v_{z,i}^{n+1/2}=v_{z,i}^{n-1/2}+A_n{a}_{z,i}^n+B_n{a}_{z,i}^{n-1},\text{ and}$$(A.2) $$z_i^{n+1}=z_i^n+v_{z,i}^{n+1/2}\Delta t^n,$$(A.3)

where $z_i^n$, $v_{z,i}^n$, and $a_{z,i}^n$ are the vertical position, velocity, and acceleration of the ith particle at the nth time-step Δtⁿ, respectively. The acceleration is computed employing cloud-in-cell mapping between the grid and the particles; in our simulations, it is due to the stellar gravity only. We explain in Appendix B how we take account of the acceleration caused by the drag of the gas onto the particles. The coefficients A_n and B_n are given by $$A_n=\frac{1}{2}\Delta t^n+\frac{1}{3}\Delta t^{n-1}+\frac{1}{6}\frac{\Delta t^{n^2}}{\Delta t^{n-1}}\text{ and}$$(A.4) $$B_n=\frac{1}{6}\left(\Delta t^{n-1}-\frac{\Delta t^{n^2}}{\Delta t^{n-1}}\right).$$(A.5)

The radial and azimuthal velocity and position components are advanced as in Cartesian coordinates, followed by a transformation from these to cylindrical coordinates:

• 1.

  The velocity is calculated as in a Cartesian geometry, i.e., inertial forces are disregarded: $$v_{\{r,\varphi \},i}^{1/2^{\prime }}=v_{\{r,\varphi \},i}^0+\frac{1}{2}{a}_{\{r,\varphi \},i}^0\Delta t^0\text{ and}$$(A.6) $$v_{\{r,\varphi \},i}^{n+1/2^{\prime }}=v_{\{r,\varphi \},i}^{n-1/2}+A_n{a}_{\{r,\varphi \},i}^n+B_n{a}_{\{r,\varphi \},i}^{n-1},$$(A.7)

  where the subscript {r, ϕ} denotes the radial or azimuthal component.

• 2.

  The position is updated in Cartesian geometry and then transformed from Cartesian to cylindrical geometry: $$x_i^{n+1}=r_i^n+v_{r,i}^{n+1/2^{\prime }}\Delta t^n,$$(A.8) $$y_i^{n+1}=v_{\varphi ,i}^{n+1/2^{\prime }}\Delta t^n,$$(A.9) $$r_i^{n+1}=\sqrt{x_i^{n+1^2}+y_i^{n+1^2}},\text{ and}$$(A.10) $${\varphi }_i^{n+1}={\varphi }_i^n+{\alpha }_i^{n+1},$$(A.11)

  where $r_i^n$ and ${\varphi }_i^n$ are the radial and azimuthal position, respectively, and the angle ${\alpha }_i^{n+1}$ can be computed as $${\alpha }_i^{n+1}=\arccos \left(\frac{x_i^{n+1}}{r_i^{n+1}}\right)=\arcsin \left(\frac{y_i^{n+1}}{r_i^{n+1}}\right).$$(A.12)

• 3.

  The velocity is corrected to reflect the transformation from Cartesian to cylindrical geometry, which entails the consideration of inertial forces: $$v_{r,i}^{n+1/2}=\cos \left({\alpha }_i^{n+1}\right)v_{r,i}^{n+1/2^{\prime }}+\sin \left({\alpha }_i^{n+1}\right)v_{\varphi ,i}^{n+1/2^{\prime }}\text{ and}$$(A.13) $$v_{\varphi ,i}^{n+1/2}=-\sin \left({\alpha }_i^{n+1}\right)v_{r,i}^{n+1/2^{\prime }}+\cos \left({\alpha }_i^{n+1}\right)v_{\varphi ,i}^{n+1/2^{\prime }}.$$(A.14)

The velocity at a full time-step can be computed from the one at a half time step as follows: $$v_i^n=v_i^{n-1/2}+\frac{1}{2}\left(A_n{a}_i^n+B_n{a}_i^{n-1}\right).$$(A.15)

We note that this full-step velocity is second-order accurate, while the error in the half-step velocity is proportional to the time-step.

A.2 Test

We test our implementation using the Leapfrog algorithm for Cartesian geometries that is already included in the FLASH Code as a benchmark: we conduct two analogous simulations of a particle in the gravitational potential of a point mass, one with a two-dimensional cylindrical and one with a three-dimensional Cartesian geometry.

The initial position of the particle, relative to the position of the point mass, can be expressed as (r, z) = (3 au, −1 au) in cylindrical and by (x, y, z) = (3 au, 0 au, −1 au) in Cartesiancoordinates. To establish an initial balance between the centrifugal and the radial gravitational force which are exerted on the particle, we set its velocity to $$\{\begin{array}{ll}{v}_{\varphi }=\sqrt{\frac{GM}{{\left(r^2+z^2\right)}^{3/2}}{r}^2}\hfill & \text{(cylindrical geometry) or}\hfill \\ v_y=\sqrt{\frac{GM}{{\left(x^2+y^2+z^2\right)}^{3/2}}\left(x^2+y^2\right)}\hfill & \text{(Cartesian geometry),}\hfill \end{array}$$(A.16)

where M = 1 M_⊙ is the point mass. Cloud-in-cell mapping is applied to calculate the gravitational acceleration of the particle. The grid cell edge length is fixed at 0.025 au.

In Fig. A.1, we show the relative error of the cylindrical radial coordinate r_cyl with respect to the coordinate $r_{\text{Car}}=\sqrt{x^2+y^2}$ in Cartesian geometry. As is evident, the error does not exceed 2 × 10⁻³ throughout the simulations.

Fig. A.1 Error of the radial coordinate rcyl that is computed by our implementation of the Leapfrog algorithm for cylindrical geometries relative to the coordinate $r_{\text{Car}}=\sqrt{x^2+y^2}$ that is calculated by the already implemented algorithm for Cartesian geometries. We conduct two equivalent simulations of a particle in a point mass gravitational potential, one applying the former and one the latter algorithm. Both simulations end after ten Keplerian orbital periods at the initial position of the particle PK,0. As can be seen, the relative error, while varying, is at most 2 × 10−3.

Appendix B: Drag

B.1 Implementation

For our implementation of the drag exerted by the gas on the dust and vice versa, we take advantage of the first-order cloud-in-cellmapping between the grid and the particles that is part of the FLASH Code.

For each of the radial, azimuthal, and vertical velocity components, we execute the following steps:

• 1.

  The gas density, pressure, and velocity component are mapped to the particles.

• 2.

  For every particle:

  • (a)

    Thefull-step dust velocity component is calculated from the stored half-step component according to Eq. (A.15).

  • (b)

    The dust stopping time t_d,stop is computed as $$t_{\text{d},\text{stop}}=\{\begin{array}{ll}\frac{{\rho }_{\text{s}}a}{{\rho }_{\text{g}}{c}_{\text{s}}}\hfill & a\le 9/4\hspace{0.17em}{\lambda }_{\text{g},\text{mfp}}\text{ (Epstein regime) or}\hfill \\ \frac{4{\rho }_{\text{s}}{a}^2}{9{\rho }_{\text{g}}{c}_{\text{s}}{\lambda }_{\text{g},\text{mfp}}}\hfill & a>9/4\hspace{0.17em}{\lambda }_{\text{g},\text{mfp}}\text{ (Stokes regime)},\hfill \end{array}$$(B.1)

    where λ_g,mfp = 1∕(σ_gn_g) = μm_H∕(σ_gρ_g) is the gas mean free path length, n_g the gas number density, σ_g = 2 × 10⁻¹⁵ cm² the molecular collision cross section (Chapman & Cowling 1970), and c_s = γP∕ρ_g the sound speed.

  • (c)

    To the stored dust half-step velocity component, the drag source term − Δv_d,drag = −(v_d − v_g)∕t_d,stop Δt is added, where v_d and v_g are the (full-step) dust and gas velocity components, respectively, and Δt is the current time-step.

• 3.

  The change in dust momentum Δp_d,drag = m_dΔv_d,drag is mapped to the grid.

• 4.

  The dragsource term Δv_g,drag = Δp_d,drag∕(ρ_gV) is added to the gas velocity in each grid cell, where V is the cell volume.

We employ the global minimum of the particle stopping time and the gas stopping time t_g,stop = ρ_g∕ρ_d t_d,stop as an upper limit of the simulation time-step.

B.1 Test

To evaluate our algorithm, we adopt the test problem introduced by Laibe & Price (2011) and applied by Bai & Stone (2010a), Mignone et al. (2019) as well as Yang & Johansen (2016), who refer to it as DUSTYBOX, particle-gas deceleration, and uniform streaming test, respectively. This problem provides an opportunity to simultaneously test the drag that is exerted by the gas on the dust and the drag back-reaction of the dust onto the gas. In addition, it allows us to compare the numeric solution computed by our algorithm to an analytic one.

We conduct simulations with a two-dimensional cylindrical geometry, as this is the geometry of the simulations we present in the main text. Dust and gas initially move in the radial direction with v_d,r,init = c_s and v_g,r,init = −c_s, respectively. The dust and gas density as well as the gas temperature are constant. The equations of motion of the dust and the gas reduce to $$\frac{\text{d}{v}_{\text{d},r}}{\text{d}t}=-\frac{v_{\text{d},r}-v_{\text{g},r}}{t_{\text{d},\text{stop}}}\text{ and}$$(B.2) $$\frac{\text{d}{v}_{\text{g},r}}{\text{d}t}=ϵ\frac{v_{\text{d},r}-v_{\text{g},r}}{t_{\text{d},\text{stop}}},$$(B.3)

where ϵ is the solid-to-gas density ratio. The dust and the gas velocity can be solved for analytically, yielding $$\begin{array}{ll}{v}_{\text{d},r}(t)=\hfill & v_{\text{d},r,\text{init}}\exp \left[-(1+ϵ)\frac{t}{t_{\text{d},\text{stop}}}\right]\hfill \\ \hfill & +v_{\text{com},r,\text{init}}\left(1-\exp \left[-(1+ϵ)\frac{t}{t_{\text{d},\text{stop}}}\right]\right)\text{ and}\hfill \end{array}$$(B.4) $$\begin{array}{ll}{v}_{\text{g},r}(t)=\hfill & v_{\text{g},r,\text{init}}\exp \left[-(1+ϵ)\frac{t}{t_{\text{d},\text{stop}}}\right]\hfill \\ \hfill & +v_{\text{com},r,\text{init}}\left(1-\exp \left[-(1+ϵ)\frac{t}{t_{\text{d},\text{stop}}}\right]\right),\hfill \end{array}$$(B.5)

where $$v_{\text{com},r,\text{init}}=\frac{v_{\text{g},r,\text{init}}+ϵv_{\text{p},r,\text{init}}}{1+ϵ}$$(B.6)

is the initial velocity of the center of mass of dust and gas.

The displacement of every dust particle relative to its initial position is given by $$\begin{array}{ll}\Delta r_{\text{d}}(t)=\hfill & r_{\text{d}}(t)-r_{\text{d},\text{init}}\hfill \\ =\hfill & \frac{(v_{\text{d},r,\text{init}}-v_{\text{com},r,\text{init}})\hspace{0.17em}{t}_{\text{d},\text{stop}}}{1+ϵ}\left(1-\exp \left[-(1+ϵ)\frac{t}{t_{\text{d},\text{stop}}}\right]\right)\hfill \\ \hfill & +v_{\text{com},r,\text{init}}\hspace{0.17em}t.\hfill \end{array}$$(B.7)

Initially, one particle is positioned at the center of every cell.

It is natural to choose the dust stopping time t_d,stop and the sound speed c_s as the units of time and velocity, respectively. The unit of length is therefore c_s t_d,stop. The domains of our simulations span 100 c_s t_d,stop in the radial dimension. The domain boundaries are periodic, which is necessary to maintain a constant dust and gas density and to conserve the total momentum of dust and gas. The simulations end after 3 t_d,stop.

We employ two quantities to measure the accuracy of our implementation: the absolute error of the total momentum with respect to the initial value and the absolute error of the particle displacement relative to the analytic solution (see Eq. (B.7)). In upper and lower panels of Fig. B.1, respectively,we show the momentum and the displacement error at the end of simulations with varying time-steps Δt (left panels), grid cell edge lengths Δx (middle panels), and solid-to-gas density ratios ϵ (right panels).

First and foremost, we note that our drag algorithm should not be used in combination with the periodic boundary implementation that is part of the FLASH Code⁵. In the figure, we mark with red crosses the errors in simulations in which particles cross the domain boundaries. The absolute error of the total momentum, while generally of the order of 10⁻¹⁰ or 10⁻⁹ otherwise, can be larger than unity in these simulations (see the upper panels).

Fig. B.1 Absolute error of the total momentum ptot of dust and gas relative to the initial total momentum ptot, init (upper panels) and mean absolute error of the simulated particle displacement Δrsim with respect to the analytic solution Δrana (see Eq. (B.7); lower panels). In the latter case, the mean is calculated by averaging over all particles, with the standard deviations plotted as error bars. (With one exception, the standard deviations are too small for the error bars to be visible, though.) Both the momentum and the displacement errors are computed after 3 td,stop. We show the errors as functions of the time-step Δt (left panels),of the grid cell edge length Δx (middle panels), and of the dust-to-gas density ratio ϵ (right panels). In the title of every column of panels, the fixed values of the other two quantities are given. Simulations in which particles cross the domain boundaries are marked with red crosses. The error in the total momentum can exceed one in these simulations, but is in general of the order of 10−10 to 10−9 otherwise. The error in the displacement increases linearly with the time-step since our drag algorithm is first-order accurate. On the other hand, it is independent of the cell size for sizes of at least 1 cs td,stop. At these resolutions, the particles are displaced by less than half a cell within 3 td,stop and therefore do not transverse the boundaries. The error for the fiducial dust-to-gas density ratio of ϵ = 0.9 amounts to 10−2; while for a ratio of ϵ = 0.01 it is greater by a factor of a few, it is as small as ~10−6 if the ratio is equal to ϵ = 100. This is despite a number of particles crossing the boundaries if the density ratio is much less or much greater than the fiducial value.

We choose ϵ = 0.9 as the fiducial dust-to-gas density ratio and Δx = 1 c_s t_d,stop as the fiducial resolution. This is to reduce the influence of the boundary conditions on our examination of the drag implementation as much as possible. For the fiducial density ratio, the analytic displacement remains less than 0.436 c_s t_d,stop, that is, less than half a cell edge length at thefiducial resolution, within 3 t_d,stop. Indeed, no particle crosses the boundaries in our simulations with this density ratio and the fiducial or a lower resolution.

From the lower-left panel of the figure, it can be seen that our implementation is first-order accurate in time. For the fiducial time-step of Δt = 10⁻², the absolute error of the displacement amounts to 1%. If the time-step is ten or a hundred times smaller, the error is less by one order or two orders of magnitude, respectively.

Reducing the resolution by a factor of ten or a hundred with respect to the fiducial resolution does not lead to an increase in the displacement error (see the lower-middle panel). This is despite the displacement within 3 t_d,stop not being resolved even at the fiducial resolution. Nonetheless, the momentum error is greater at lower resolutions; it is of the order of 10⁻⁷ in the simulation with a resolution of a hundredth of the fiducial resolution, in which only one particle is present. On the other hand, if the resolution is higher than the fiducial one, the error in both the displacement and the total momentum is considerable. This is because a large number of particles transverse the domain boundaries.

Compared to the error for the fiducial dust-to-gas density ratio, the displacement error increases by a factor of a few for low density ratios, but decreases by orders of magnitude for high ratios. This is evident from the lower-right panel, in which the errors in simulations with the fiducial density ratio as well as ϵ = 0.01, ϵ = 1, and ϵ = 100 are shown. A number of particles cross the boundaries in both the simulation with the highest density ratio and the one with the lowest ratio. Thus, the error in the total momentum, though not the error in the displacement, is significantly greater for these ratios than for the fiducial one.

References

All Tables

All Figures

Fig. 1 Mach number of the vertical gas velocity ${\mathcal{M}}_{\text{g},z}$ as a functionof radius r and height z in a simulation including a locally isothermal gas, but no particles (iso_Lr=90au_Lz=4Hg). The simulation domain spans four gas scale heights in the vertical dimension. Black lines mark one gas scale height above and below the midplane. After 30 kyr, the vertical shear instability has saturated at all radii and heights; it gives rise to perturbations which are characterized by a radial-to-vertical wavelength ratio much smaller than unity, by being bent outwards, and by a mirror symmetry relative to the midplane.

In the text

Fig. 2 Root mean square of ${\mathcal{M}}_{\text{g},z}$ as a functionof time t. The root mean square is computed using the mass-weighted average over the vertical domain size and over 1 au extending from r = 50  to 51 au. In the domain with a vertical dimension of four gas scale heights (blue and orange lines), the Mach number saturates at ${\mathcal{M}}_{\text{g},z,\text{rms}}\approx {10}^{-1}$. This value is consistent with the one which Flock et al. (2017) find in their numerical study of the vertical shear instability. The turbulent strength is similar in the simulation with the fiducial resolution (iso_Lr=90au_Lz=4Hg; blue line) and the one with the doubled resolution (iso_Lr=90au_Lz=4Hg_dou_res; orange line); however, it remains lower until the simulations end in the domain with a vertical size of two gas scale heights (simulation iso_Lr=90au; green line).

In the text

Fig. 3 Mach number ${\mathcal{M}}_{\text{g},z}$ as a functionof r and z in two simulations with dust and an isothermal gas equation of state. As in the dust-free equivalent (compare with Fig. 1), the vertical shear instability induces large-scale perturbations away from the midplane. In the simulation iso_Z=0.02_Lr=90au_Lz=4Hg (left panel), in which we initialize the dust at the beginning, small-scale perturbations can be seen at all radii in the midplane. These perturbations are also observable in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr (right panel), in which the dust is introduced after 50 kyr, although not at all radii and extending to greater heights. The inlays with a vertical extent of 0.6 gas scale heights show the enlarged perturbations. However, such small-scale perturbations do not exist in the corresponding simulation without dust.

In the text

Fig. 4 Root mean square of ${\mathcal{M}}_{\text{g},z}$, averaged over the vertical domain extent and the radii given in the legend, as a function of t after the initialization of the dust at td,init. The mean is weighted by mass, i.e., more weight is assigned to the Mach number in the midplane than to the one away from it. In the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr (left panel), the vertical shear instability has saturated before the dust is introduced. The Mach number remains approximately constant at a value ${\mathcal{M}}_{\text{g},z,\text{rms}}\approx {10}^{-1}$ after the dust initialization (compare with Fig. 2). That is, in this simulation the vertical shear instability is the primary source of turbulence in the midplane. In contrast, in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg (right panel), in which the vertical shear instability and the streaming instability begin to grow at the same time, the Mach number saturates at a lower value of ${\mathcal{M}}_{\text{g},z,\text{rms}}\approx {10}^{-2}$. After no dust remains at a radius because of its radial drift (this point is marked with a dashed line), the Mach number increases until it is approximately equal to the value at the same radius in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr. We find that, as long as dust is present, the streaming instability drives the turbulence in the midplane in the simulation iso_Z=0.02_Lr=90au_Lz=4Hg.

In the text

Fig. 5 Ratio of the root mean square Mach number of the radial gas velocity ${\mathcal{M}}_{\text{g},r}$ to the root mean square of ${\mathcal{M}}_{\text{g},z}$ as a functionof z in units of gas scale heights. The abscissa is scaled linearly between − 0.1 and 0.1 gas scale heights and logarithmically otherwise. To compute the root mean square, we take the mass-weighted average over r = 50  to 51 au and the average over 1 kyr through t = 55 kyr in thecase of the simulation iso_Z=0.02_Lr=90au_Lz=4Hg_tdinit=50kyr (blue line) and through t = 10 kyr in thecase of the simulation iso_Z=0.02_Lr=90au_Lz=4Hg (orange line). Away from the mid-plane, the radial-to-vertical Mach number ratio is less than unity in both simulations. While this anisotropy extends over all heights if the dust is introduced after 50 kyr, the Mach number ratio is close to one in the midplane if the dust is initialized at the start. This is indicative of the vertical shear instability being the primary source of turbulence in the midplane in the former case, but the streaming instability in the latter case.

In the text

Fig. 6 Root mean square of ${\mathcal{M}}_{\text{g},z}$ as a functionof the dust-to-gas density ratio ρd∕ρg. The orange, green, and purple lines represent the turbulent strength in simulations in which we initialize the dust after the vertical shear instability has saturated. The strength at low to intermediate values of the volume density ratio in the three simulations is marginally smaller if the dust-to-gas surface density ratio is greater. Nevertheless, the Mach number is higher in all of these simulations than in the simulation in which we introduce the dust at the start (blue line). We find that the stronger turbulence in the former simulations is driven by the vertical shear instability, while the weaker turbulence in the latter one is caused by the streaming instability. Independent of the source of the turbulence, its strength decreases at dust densities larger than the gas density in all four simulations.

In the text

Fig. 7 Ratio of the dust scale height to the gas scale height Hg (left panel) and maximum of the dust-to-gas density ratio ρd∕ρg (right panel)as functions of t after the dust is introduced at td,init. We compute the dust scale height as the root mean square of the vertical particle positions zd, and average the scale height ratio over 1 au spanning from r = 20  to 21 au. For a dust-to-gas surface density ratio of Z = 2%, the dust scale height amounts to ~1% of the gas scale height if the streaming instability induces the vertical diffusion of the dust (blue line). On the other hand, it is equal to ~10% for the same surface density ratio if the diffusion is caused by the vertical shear instability (orange line). The scale height which is induced by the vertical shear instability decreases with increasing surface density ratio (green and purple line). However, it is higher than the value that the streaming instability gives rise to for all surface density ratios that we consider. In contrast, the maximum dust volume density is significantly smaller if the streaming instability drives the turbulence in the dust layer than if the turbulence is caused by the vertical shear instability. In the former case, the dust-to-gas volume density ratio amounts to a few hundred for a surface density ratio of 2%. In the latter case, on the other hand, it exceeds 103 for the same surface density ratio, and 104 for a surface density ratio of 10%. In all cases, the maximum dust density is greater than the Roche density ρR in the midplane at the inner radial boundary of the simulation domains, i.e., the maximum Roche density in the simulations. The ratio of this Roche density to the gas density is marked as a black line.

In the text

Fig. 8 Dust-to-gas volume density ratio ρd∕ρg as a functionof r and z (upper panels) as well as dust-to-gas surface density ratio Σd∕Σg as functions of r and t (lower panels). The upper panels show the spatial dust distribution 5 kyr after the dust initialization. We compare a simulation in which the turbulence in the dust layer is driven by the streaming instability (left panels) and two simulations in which the vertical shear instability is the main source of turbulence (middle and right panels). For the same dust-to-gas surface density ratio of Z =0.02, the dust scale height is smaller, but the radial dust concentration is weaker in the former case. (The values of Z are specified in the titles.) In the latter case, more dust is accumulated in overdensities if the surface density ratio is higher. Some of the accumulations are sufficiently dense for their radial drift to cease almost entirely.

In the text

Fig. 9 Mach number of the vertical gas velocity as a function of r and z in simulations of dust and a locally adiabatic gas. The Mach number at radii of less than ~ 3 au is similarlyhigh in the simulation without (adi_Lr=9au; left panel) or with dust (adi_Z=0.02_Lr=9au; right panel). This is in spite of the streaming instability operating in the latter simulation, but not in the former one. At larger radii in the simulation that includes dust, the streaming instability induces large-scale perturbations with ${\mathcal{M}}_{\text{g},z}\approx {10}^{-3}$ away from the midplane. These perturbations are similar in shape to the ones caused by the vertical shear instability, but they are weaker, are bent inwards rather than outwards, and are not symmetric with respect to the midplane (compare with Fig. 1). In the inlay, which extends to 0.1 gas scale heights Hg above and below the midplane, small-scale perturbations with ${\mathcal{M}}_{\text{g},z}\approx {10}^{-2}$ in the midplane can be seen.

In the text

Fig. 10 Ratio of the dust scale height, as the root mean square of zd, to the gas scale height Hg (left panel) as well as mass-weighted root mean square of ${\mathcal{M}}_{\text{g},z}$ in the midplane (right panel). Both quantities are depicted as functions of r. To compute the root mean square, we average over 50  or 500 yr after the dust scale height has reached an equilibrium value in the simulation domains spanning 1 au ≤ r ≤ 10 au or 10 au ≤ r ≤ 100 au, respectively. The dashed line marks the boundary between these domains. The model with an initial dust-to-gas scale height ratio of ${(H_{\text{d}}/H_{\text{g}})}_{\text{init}}=0.1$, a dust-to-gas surface density ratio of Z = 2%, and the fiducial resolution (blue line) is shown together with the models that deviate from this fiducial one in that ${(H_{\text{d}}/H_{\text{g}})}_{\text{init}}=1$, that Z = 1% (green line),or that the initial and maximum resolution are doubled (purple line). While both the dust scale height and the Mach number depend on the radius, they are largely independent of the initial dust scale height, the surface density ratio, and the resolution.

In the text

Fig. 11 Rootmean square of ${\mathcal{M}}_{\text{g},z}$ in the midplane as a function of r (lower abscissa) and the dimensionless gas pressure gradient parameter Π in the midplane (upper abscissa). The root mean square is calculated using the mass-weighted mean over 500 yr after an equilibrium value of the dust scale height is reached. In the simulation adi_Z=0.02_Lr=90au_taustop=0.1 (blue line), in which the dimensionless dust stopping time τstop is fixed at 0.1, the Mach number increases with the pressure gradient parameter if Π ≲ 0.09, but is approximately constant if the pressure gradient is stronger. The Mach number in this one and the simulation adi_Z=0.02_Lr=90au (green line), in which the stopping time increases from 0.046 to 0.46, is comparable at all radii. That is, the turbulent strength does not depend significantly on the stopping time in this stopping time regime. The Mach number is lower overall, but depends more strongly on the radius in the simulation adi_Z=0.02_Lr=90au_a=3mm (orange line) in which the stopping time ranges from 0.0046 to 0.046. In this regime, the turbulent strength increases with both the strength of the pressure gradient and the stopping time.

In the text

Fig. 12 Radial gas velocity (solid line) as a function of r in the simulations adi_Z=0.02_Lr=9au and adi_Z=0.02_Lr=90au. (The dashed lines marks the boundary between the domains of the two simulations.) The velocity is computed as the mass-weighted average over the vertical domain extent and a time-span of 50  and 500 yr, respectively,after an equilibrium dust scale height has been attained in each simulation. The average radial velocity at r ≥ 3 au is plotted as a dotted line and amounts to 0.035 m s−1. This outward motion is caused by the dust drifting radially inwards and its angular momentum being transferred to the gas. However, it is evident that the mean of the velocity is less than its standard deviation, which is equal to 0.16 m s−1. This is a result of the turbulent motions caused by the streaming instability.

In the text

Fig. A.1 Error of the radial coordinate rcyl that is computed by our implementation of the Leapfrog algorithm for cylindrical geometries relative to the coordinate $r_{\text{Car}}=\sqrt{x^2+y^2}$ that is calculated by the already implemented algorithm for Cartesian geometries. We conduct two equivalent simulations of a particle in a point mass gravitational potential, one applying the former and one the latter algorithm. Both simulations end after ten Keplerian orbital periods at the initial position of the particle PK,0. As can be seen, the relative error, while varying, is at most 2 × 10−3.

In the text

Fig. B.1 Absolute error of the total momentum ptot of dust and gas relative to the initial total momentum ptot, init (upper panels) and mean absolute error of the simulated particle displacement Δrsim with respect to the analytic solution Δrana (see Eq. (B.7); lower panels). In the latter case, the mean is calculated by averaging over all particles, with the standard deviations plotted as error bars. (With one exception, the standard deviations are too small for the error bars to be visible, though.) Both the momentum and the displacement errors are computed after 3 td,stop. We show the errors as functions of the time-step Δt (left panels),of the grid cell edge length Δx (middle panels), and of the dust-to-gas density ratio ϵ (right panels). In the title of every column of panels, the fixed values of the other two quantities are given. Simulations in which particles cross the domain boundaries are marked with red crosses. The error in the total momentum can exceed one in these simulations, but is in general of the order of 10−10 to 10−9 otherwise. The error in the displacement increases linearly with the time-step since our drag algorithm is first-order accurate. On the other hand, it is independent of the cell size for sizes of at least 1 cs td,stop. At these resolutions, the particles are displaced by less than half a cell within 3 td,stop and therefore do not transverse the boundaries. The error for the fiducial dust-to-gas density ratio of ϵ = 0.9 amounts to 10−2; while for a ratio of ϵ = 0.01 it is greater by a factor of a few, it is as small as ~10−6 if the ratio is equal to ϵ = 100. This is despite a number of particles crossing the boundaries if the density ratio is much less or much greater than the fiducial value.

In the text