© The Authors 2025

1 Introduction

Observations reveal that protoplanetary disks are turbulent, with this turbulence generally being subsonic (e.g. Rosotti 2023). Up to supersonic turbulent speeds are only reached in the upper disk layers or very close to the central star (Carr et al. 2004; Najita et al. 2009; Doppmann et al. 2011). This renders turbulence measurements based on the broadening of molecular emission lines challenging since thermal and turbulent broadening are difficult to disentangle.

Nevertheless, owing in particular to the unprecedented spatial and spectral resolution of the Atacama Large Millimeter/submillimeter Array (ALMA), line broadening observations have yielded turbulent strengths in the outer regions of a handful of disks. In the disks surrounding HD 163296, TW Hya, MWC480, and V4046 Sgr, these strengths have been constrained to be at most ~1% of the sound speed (Flaherty et al. 2015, 2017, 2018; Flaherty et al. 2020; Teague et al. 2018). Yet, in DM Tau and IM Lup disks Mach numbers of 0.2–0.6 have been inferred (Flaherty et al. 2020; Rosotti 2023; Paneque-Carreño et al. 2024; Flaherty et al. 2024). Moreover, ALMA observations have enabled not only measurements of the magnitude of gas turbulence, but also the detection of structures in the gas kinematics (e.g. Pinte et al. 2023). Although structures like meridional flows (Teague et al. 2019; Yu et al. 2021) have been interpreted as indications of the presence of planets, they could also be an imprint of the vertical shear instability (Barraza-Alfaro et al. 2021, 2024).

Further observational constraints on the strength of gas turbulence can be obtained indirectly from the vertical and radial distribution of dust. These are degenerate with the strength of the drag coupling between dust and gas, however. The dust scale height, on the one hand, constitutes a means to measure gas turbulence since it is regulated by the vertical stellar gravity and turbulent diffusion. From dust-to-gas scale height ratios of at most 0.1, Pinte et al. (2016) and Villenave et al. (2022) infer turbulent speeds of ~1% of the sound speed or less in the disks around HL Tau and Oph163131, respectively (see also Ueda et al. 2021). Dust scale heights ranging between 0.1 and one gas scale height as well as Mach numbers as low as ≲ 0.05 and as high as 0.1–0.5 are derived in different parts of the ring-andgap structure of the HD 163296 disk (Ohashi & Kataoka 2019; Doi & Kataoka 2021; Liu et al. 2022) – complementary with the estimates for this disk obtained from line broadening.

Similar to the scale height, the radial width of the dust rings observed in a number of disks has been utilised to measure gas turbulence (Dullemond et al. 2018; Sierra et al. 2019; Rosotti et al. 2020; Facchini et al. 2020; Carvalho et al. 2024). If the rings are associated with gas pressure maxima, turbulent diffusion counteracts the tendency of the dust to drift towards the centre of the rings. This approach has yielded a Mach number of 0.1 in the disk surrounding HD 169142 (Sierra et al. 2019) and values that are comparable or smaller by up to an order of magnitude in the disks around AS 209, HD 163296, and Elias 2-24 (Rosotti et al. 2020; Zagaria et al. 2023; Carvalho et al. 2024). Differences in the turbulent strengths derived from the radial and vertical dust distributions might be indicative of anisotropic turbulence (Doi & Kataoka 2021; Villenave et al. 2022; Zagaria et al. 2023).

Recent years have seen a rapid development of our understanding of protoplanetary disks turbulence not only from an observational perspective. Theoretically, a variety of instabilities have been established as potential sources and a map of the dominant mechanisms to drive turbulence at different disk radii and heights is emerging (e.g. Lesur et al. 2023). Nevertheless, this picture is not yet comprehensive as we are only beginning to apprehend the intricacies of disk turbulence, including the interaction of different instabilities and how this interplay varies with length scale.

In Schäfer et al. (2020, hereafter SJB20) and Schäfer & Johansen (2022, hereafter SJ22), we explore the coexistence of two among the most promising instabilities, the vertical shear instability and the streaming instability. The former is a hydrodynamic instability that accesses the free energy associated with – as the name indicates – vertical shear (Barker & Latter 2015). Vertical shear arises, for instance, from baroclinity, that is to say a misalignment between the density and pressure gradients. Owing to temperature gradients, protoplanetary disks are indeed generally baroclinic. Since it is generally dampened by vertical buoyancy, the instability requires gas cooling to be sufficiently rapid (Nelson et al. 2013; Lin & Youdin 2015). This condition is met in the outer disk regions at radii ≳ 10 au (Lin & Youdin 2015; Malygin et al. 2017; Pfeil & Klahr 2019). We note that the linear analysis conducted by Latter & Papaloizou (2018) shows that the vertical shear instability experiences growth even at the longer cooling times relevant in optically thick regions, albeit with a reduced growth rate.

The vertical shear instability drives anisotropic turbulence, with the turbulent gas speeds being much higher in the vertical than in the radial dimension (Nelson et al. 2013; Stoll et al. 2017; SJB20). In models including only the drag coupling of the dust to the gas (not vice versa), the vertical gas motions stir up dust to a scale height that is comparable to the gas scale height even if the Stokes number of the dust – the ratio of the stopping time, the timescale of this drag coupling, to the dynamical timescale – is as large as 0.1 (Stoll & Kley 2016; Flock et al. 2017; Dullemond et al. 2022). However, if the drag that the dust exerts on the gas is taken into account, the tendency of the dust to sediment to the disk mid-plane ‘weighs’ down the gas. This effective buoyancy causes a reduction of both the vertical turbulent gas velocities and the dust scale height (Lin & Youdin 2017; Lin 2019; SJB20), the latter to ~10% of the gas scale height for the above Stokes number and the canonical dust-to-gas surface density ratio of 1% (Lin 2019).

The streaming instability, as originally discovered by Youdin & Goodman (2005), arises where dust and a (negative) radial gas pressure gradient are present, that is all but everywhere in the dust layer of protoplanetary disks. The pressure gradient is the source of free energy that the instability taps into (Youdin & Johansen 2007) as it causes isolated gas to orbit with a subKeplerian speed, in contrast with the Keplerian rotation of dust in isolation¹. The mutual drag between dust and gas rotating with different speeds results in a relative radial drift and instability.

In addition to arguably being the leading candidate among mechanisms to induce planetesimal formation (e.g. Johansen et al. 2014; Lesur et al. 2023), the streaming instability is by nature also a source of turbulence. The strength of both this turbulence and the dust diffusion it entails increases with the Stokes number, as do the length scales of turbulence (Johansen & Youdin 2007; Schreiber & Klahr 2018; SJB20; Yang & Zhu 2021; Baronett et al. 2024). This is in agreement with larger linear growth rates and wavelengths for higher Stokes number (Youdin & Goodman 2005). The instability drives largely isotropic turbulence when considering the radial and vertical dimensions (Johansen & Youdin 2007; SJB20; Yang & Zhu 2021; Baronett et al. 2024).

Both Schreiber & Klahr (2018) and Li & Youdin (2021) report that dust diffusion, on the other hand, is stronger radially than vertically, although this discrepancy disappears with increasing dust-to-gas ratio (Schreiber & Klahr 2018; Gerbig & Li 2023). In contrast, the work by Yang & Zhu (2021) and Baronett et al. (2024) reveals the opposite for Stokes numbers of order unity, and no anisotropy for smaller Stokes numbers. We note that of these studies only Li & Youdin (2021) includes the vertical stellar gravity. The dust scale height amounts to ~1% of the gas scale height (Bai & Stone 2010b; Yang & Johansen 2014; Li et al. 2018; SJB20) and is self-regulatory: If the dust would settle to a thinner layer, the density of this layer would be higher, and the instability would cause stronger turbulence and diffusion (Bai & Stone 2010b).

Besides models of only the vertical shear instability or the streaming instability, we explore two scenarios in which both instabilities are active in SJB20; SJ22, and in this paper. In the scenario SIwhileVSI, both instabilities begin to grow simultaneously. Here, turbulence in the mid-plane dust layer is predominantly caused by the streaming instability, and away from this layer by the vertical shear instability. In the scenario SIafterVSI, on the other hand, the vertical shear instability saturates before the streaming instability starts to operate, and remains the primary source of turbulence both in and away from the dust layer. Nonetheless, the streaming instability causes turbulence locally in dust overdensities also in this scenario. Expanding on our previous study focused on gas turbulence, we aim to explore dust diffusion induced by the two instabilities in isolation and in conjunction in this paper.

This paper is structured as follows: we introduce our numerical model and its parameters in Sect. 2. In Section 3, we examine gas turbulence induced by the vertical shear instability and the streaming instability individually and jointly, before commenting on the shape of the dust layer in Sect. 4. Sections 5 and 6 are dedicated to a detailed investigation of vertical and radial dust diffusion, respectively. This is followed by a comparison between our models and the one presented by Flock et al. (2020). We discuss implications and limitations of our work in Sect. 8, and conclude in Sect. 9.

2 Simulations

The numerical model employed in this study is the same as in SJB20 and SJ22. We provide a summary of the model below, and refer to SJB20 for a more detailed description.

We applied the FLASH Code²,³ (Fryxell et al. 2000) to simulate the gas and dust components of protoplanetary disks, the former on a Eulerian grid and the latter as Lagrangian particles. Our model includes the drag coupling between the two components as well as stellar gravity. In Table 1, we list all simulations and the parameters that set them apart. The simulation names are composed of the scenario modelled and the dust size a, the latter given in centimetres.

To begin with, we consider scenarios in which only either the vertical shear instability or the streaming instability is active. To ensure that the respective other instability does not operate, in the former scenario we consider only the drag exerted by the gas on the dust. (The drag of the dust onto the gas was included in all other scenarios.) In the latter scenario, the gas equation of state is adiabatic since under this condition the vertical shear instability is quenched by vertical buoyancy. (An isothermal equation of state was used in all other scenarios.) Moreover, we simulate the two instabilities in conjunction. In the scenario SIwhileVSI, both the streaming instability and the vertical shear instability are in operation from the onset, while in scenario SIafterVSI dust particles are introduced and the streaming instability therefore begins to grow only after 50 kyr. At this point, the vertical shear instability has attained a saturated state everywhere in the simulation domain (SJB20).

2.1 Domain size, resolution, and boundary conditions

We performed two-dimensional, axisymmetric simulations with a cylindrical geometry, which is a natural choice to represent protoplanetary disks. The simulation domains extend from 10 au to 100 au in the radial dimension and to either 1 or 2 gas scale heights above and below the mid-plane in the vertical dimension. They are thus shaped like isosceles trapezoids with curved legs as the gas scale height increases non-linearly with the radius (see Eq. (2)). We chose the larger vertical domain size for all models involving the vertical shear instability because we find this to be necessary for the instability to drive turbulence with a strength that is consistent with previous studies (SJB20). Since the streaming instability induces a dust scale height of ~1% of the gas scale height (Bai & Stone 2010b; Yang & Johansen 2014; Li et al. 2018; SJB20), the smaller vertical domain size is sufficient for our model of this instability only.

The base resolution of our simulations amounts to 10 grid cells per astronomical unit, or 8 and 150 cells per gas scale height in the mid-plane at the inner and outer radial domain boundaries, respectively (see Eq. (2)). On top of that, we used static and adaptive mesh refinement with up to six levels of refinement. Every refinement level corresponds to a doubling in resolution, with the maximum resolution thus being equal to 320 cells per astronomical unit or 270 and 4800 cells per scale height at the radial boundaries. Static refinement was applied to increase the resolution by a factor of two within 5 au and by an additional factor of two (for a total factor of four) within 1 au above and below the mid-plane. Adaptive mesh refinement, on the other hand, enhanced the resolution of blocks consisting of 10 × 10 cells if the number of particles in any cell within these blocks exceeded ten. Conversely, the resolution of a block was reduced if no particles remained in a cell. This approach to refinement enables us to simulate protoplanetary disks on a global scale, with a vertical domain extent that is sufficient to model the vertical shear instability, while simultaneously resolving diffusion in the dust mid-plane layer.

At the radial and vertical boundaries, gas and dust particles were allowed to move out of but not into the domains. To prevent the boundary conditions from artificially affecting our results, we excluded the innermost 10 au of the domains from all quantitative analysis. The vertical boundaries and the outer radial one are highly unlikely to influence our study of dust diffusion because the dust rapidly settles vertically and drifts radially away from these boundaries.

2.2 Gas

The gas is initially in vertical hydrostatic equilibrium, with the stellar gravity being balanced by a density gradient given by $${\rho }_{\mathrm{g}}={\rho }_{\mathrm{g}}(z=0)\exp [-\frac{\gamma G{M}_{\mathrm{S}}}{c_{\mathrm{S}}^2}(\frac{1}{r}-\frac{1}{\sqrt{r^2+z^2}})],$$(1)

where r and z are the radial and vertical coordinate, respectively, with r = 0 and z = 0 being the coordinates of the central star and the disk mid-plane, respectively. Furthermore, γ is the adiabatic index, G the gravitational constant, M_S = 1 M_⊙ the stellar mass, c_s = (γRT / μ)^1/2 the sound speed, R the ideal gas constant, T the temperature, and μ = 2.33 the mean molecular weight. The temperature initially does not vary with height. The gas scale height is equal to $$H_{\mathrm{g}}=\sqrt{\frac{c_{\mathrm{s}}^2{r}^3(2\gamma G{M}_{\mathrm{S}}-c_{\mathrm{s}}^{2}r)}{{(c_{\mathrm{s}}^{2}r-\gamma G{M}_{\mathrm{S}})}^2}}=0.846\mathrm{a}\mathrm{u}{\left(\frac{r}{10\mathrm{a}\mathrm{u}}\right)}^{5/4}.$$(2)

The initial radial gas density and temperature gradients can be expressed as $${\rho }_{\mathrm{g}}(z=0)=5.62\times {10}^{-12}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}{\left(\frac{r}{10\mathrm{a}\mathrm{u}}\right)}^{-9/4}\text{ and }$$(3) $$T=88.5\mathrm{K}{\left(\frac{r}{10\mathrm{a}\mathrm{u}}\right)}^{-1/2}.$$(4)

The latter is adopted from the minimum mass solar nebula model (Hayashi 1981), while the surface density gradient Σ_g ∝ r⁻¹ is consistent with observed ones (Andrews et al. 2009, 2010). The negative radial pressure gradient and the baroclinity resulting from the temperature gradient give rise to the streaming instability and the vertical shear instability, respectively (Youdin & Goodman 2005; Nelson et al. 2013).

As noted above, we employed an isothermal equation of state, $$P=\frac{RT}{\mu }{\rho }_{\mathrm{g}},$$(5)

in models of the vertical shear instability in isolation or in coexistence with the streaming instability since under this condition there is no vertical buoyancy counteracting the vertical shear instability (Nelson et al. 2013; Lin & Youdin 2015). (In these models, the adiabatic index γ = 1.) Conversely, in our model of only the streaming instability buoyancy quenches the vertical shear instability since the equation of state is adiabatic, $$P=K{\rho }_{\mathrm{g}}^{\gamma },$$(6)

where $K=RT{\rho }_{\mathrm{g}}^{1-\gamma }/\mu$ is the polytropic constant and the adiabatic index γ = 5/3. We note that the polytropic constant is constant in time only as it varies with density and temperature. The gas initially orbits with a sub-Keplerian velocity such that the centrifugal and the pressure gradient force balance the radial stellar gravity.

2.3 Dust particles

We adopted an approach to modelling dust aggregates in protoplanetary disks using Lagrangian particles that is commonly applied in simulations of the streaming instability (Youdin & Johansen 2007; Bai & Stone 2010a): The mass and momentum of each particle are equal to those of a huge number of aggregates, while the drag coupling to the gas is equal to that of a single aggregate. In every simulation, 10⁶ particles were initialised with a radially uniform distribution, while their vertical positions were randomly sampled from a Gaussian distribution whose scale height amounts to 10% of the gas scale height. The noise in this vertical distribution constitutes the seed for the streaming instability. The particles were introduced at the start of simulations in our model of only the streaming instability and in the scenario SIwhileVSI, but not until after 50 kyr in our model of the vertical shear instability in isolation and in the scenario SIafterVSI. Our results are independent of whether 5 × 10⁵ or 10⁶ particles were simulated and of the random seed, while the dependence on the initial scale height is examined in SJB20.

All dust particles possess the same mass m_d = 2.53 × 10²⁴ g. Since the particles are initially uniformly distributed in the radial dimension, their mass can be calculated as $$m_{\mathrm{d}}=\frac{1}{N_{\mathrm{d}}}{\int }_{L_r}2\pi r{\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{d}r.$$(7)

Here, N_d = 10⁶ is the total particle number, L_r = 90 au the radial domain size, and Σ_d the dust surface density. Because we set a constant ratio of dust to gas surface density and the latter is inversely proportional to the radius, the particle mass is radially constant. The surface ratio was chosen to amount to 2%. We note that this value is twice as high as the canonical dust-to-gas ratio in the interstellar medium, but lies within the range of dust-togas mass ratios obtained from disk observations (e.g. Miotello et al. 2023). The enhanced surface density ratio does not affect gas turbulence and dust diffusion in our model of the vertical shear instability in isolation since the drag of the dust onto the gas is not considered. Similarly, simulations of the streaming instability with a lower surface density ratio of 1% yield comparable turbulent strengths and dust scale heights to the 2% case (SJB20). Nonetheless, it is possible that dust clumping, which is stronger for higher surface density ratios, influences diffusion (Li & Youdin 2021). Furthermore, a higher dust-to-gas ratio results in a stronger effective buoyancy suppressing the vertical shear instability in the scenarios in which both instabilities coexist (Lin & Youdin 2017; Lin 2019; SJB20).

We conducted two simulations of every scenario, one with a dust size of a = 3 mm and one with a = 3 cm. These sizes are consistent with the maximum ones derived from the opacity spectral index of the thermal dust emission from protoplanetary disks (Macías et al. 2019, 2021; Carrasco-González et al. 2019; Tazzari et al. 2021a,b; Maucó et al. 2021). Because both dust sizes are smaller than the gas mean free path length everywhere in our model, that is to say Epstein drag is applicable, they can be converted to Stokes numbers in the disk mid-plane as $$\begin{array}{rl}\mathrm{St}(z=0)& =\frac{a{\rho }_{\mathrm{s}}}{c_{\mathrm{s}}{\rho }_{\mathrm{g}}(z=0)}{\mathit{\Omega }}_{\mathrm{K}}(z=0)\\ & =6\times {10}^{-3}\frac{1}{\sqrt{\gamma }}\left(\frac{a}{3\mathrm{m}\mathrm{m}}\right)\left(\frac{r}{10\mathrm{a}\mathrm{u}}\right).\end{array}$$(8)

Here, ρ_s = 1 g cm⁻³ is the dust material density and Ω_K the Keplerian orbital frequency. This conversion relation is illustrated in Fig. 2 of SJ22. The initial orbital velocity of the dust particles is Keplerian.

3 Gas turbulence

While the focus of this study is on dust diffusion, in this section we briefly discuss the gas turbulence that induces this diffusion. To this end, we provide a synopsis of the findings of SJB20 and complement these with new results.

In SJB20, we report that the vertical shear instability causes anisotropic turbulence as the gas speeds are higher in the vertical than in the radial dimension (Nelson et al. 2013; Stoll et al. 2017). The Mach number of the former speeds amounts to ℳ_g,z ≈ 0.1 (Flock et al. 2017). The streaming instability, on the other hand, is a source of largely isotropic turbulence with a Mach number between 0.001 and 0.01 (Johansen & Youdin 2007; Flock & Mignone 2021; Li & Youdin 2021; Yang & Zhu 2021; Baronett et al. 2024).

These findings are reflected in Fig. 1. Here, we depict the evolution of the kinetic energy associated with the radial and vertical components of the gas velocity as well as the sum of the two components in all our scenarios. We note that the figure shows the evolution beginning with the time at which the dust particles are initialised, that is to say after the vertical shear instability has saturated in our model including only this instability and in the scenario SIafterVSI.

While the kinetic energy is overall greatest in our model of the vertical shear instability in isolation, it is more than an order of magnitude larger when considering only the vertical component than when taking into account only the radial one. In comparison, the streaming instability in isolation gives rise to turbulence with the least kinetic energy, with the energy of the radial gas motions being marginally larger than that of the vertical motions. It is interesting to note a spike in the kinetic energy after ~2 kyr in this model. We do not find a peak or a valley in the maximum or mean dust-to-gas density ratio at the same time, that is to say the abrupt increase in turbulent strength is not related to the instability suddenly inducing stronger or weaker dust concentration (see also SJ22).

In the scenario SIwhileVSI, both the streaming instability and the vertical shear instability are active from the outset. The gas kinetic energy is initially comparable to that in the model of the streaming instability in isolation, with the distribution among the radial and vertical velocity components being roughly equal. This shows that the streaming instability grows more rapidly in energy at least close to the mid-plane where the density is highest. Subsequently, though, the energy in the scenario SIwhileVSI increases to exceed that in the model of the streaming instability only. Moreover, considerably more energy is associated with the vertical than the radial velocity component. That is, at later times the kinetic energy is dominated by the vertical shear instability. This is despite the turbulence in the dust mid-plane layer being primarily caused by the streaming instability (SJB20).

The vertical shear instability is the dominant source of gas turbulence in the scenario SIafterVSI, both in and away from the mid-plane (SJB20). Thus, the total kinetic energy and the energies of the radial and vertical motions are similar to those in our model of the isolated vertical shear instability. Over time, though, particularly the total energy and the energy in the vertical velocities decline. We interpret this as a consequence of the settling dust causing an effective vertical buoyancy that weakens the vertical shear instability (Lin & Youdin 2017; Lin 2019).

It is interesting to note that the kinetic energies in the scenarios SIwhileVSI and SIafterVSI attain comparable values ~10 kyr after dust initialisation. This is not true for the strength of dust diffusion especially in the vertical dimension, though, as can be seen from Fig. 5 of SJ22 and is detailed below.

In Appendix A, we present and analyse kinetic energy power spectra. We note, though, that their meaningfulness is limited by the grid in our simulations being neither periodic nor uniform.

Fig. 1 Total gas kinetic energy Eg,kin in the region between r = 20 au and 30 au and restricted to one gas scale height off the mid-plane as a function of time t after the initialisation of the dust particles at td,init. Differently coloured lines represent simulations of the streaming instability and the vertical shear instability in isolation as well as of the scenarios SIwhileVSI and SIafterVSI, respectively, with the dust size being equal to 3 cm in all simulations. Like the model including only the vertical shear instability, the scenario SIafterVSI as well as the scenario SIwhileVSI at late times are characterised by a comparatively high total energy, with the energy of the vertical motions being substantially higher than that of the radial motions. In comparison, the total energy is less in the model of the streaming instability and in the scenario SIwhileVSI at early times, but largely equally distributed among the radial and vertical velocity components.

4 Dust layer thickness

In Figure 2, we show the ratio of dust to gas density in our simulations of the vertical shear instability in isolation and of the scenario SIafterVSI with a dust size of 3 cm. The scale height of the dust, that is to say the height above or below the disk mid-plane to which the dust is elevated owing to large-scale turbulent diffusion, amounts to ~10% of the gas scale height in both simulations.

However, a stark contrast in the thickness of the dust layer is evident. The vertical shear instability alone causes a minuscule amount of small-scale diffusion both in the vertical and in the radial dimension, resulting in a dust layer that is not resolved in our model despite the application of adaptive mesh refinement. In comparison, in the scenario SIafterVSI the dust layer is well-resolved because the streaming instability induces diffusion within the layer in this scenario.

It is important to note here that, although at first glance the much denser dust layer in our simulation of the vertical shear instability in isolation appears to provide better conditions for gravitational collapse and planetesimal formation, the drag exerted by the dust onto the gas is not included and the streaming instability thus not active in this model. However, neither this drag nor the streaming instability can be neglected when studying dust concentration and planetesimal formation.

In light of this finding, in the following we explore dust diffusion with a focus on its scale-dependence in all of our models. We begin by discussing the diffusion in the vertical dimension in the following section, with Sect. 6 being dedicated to the radial diffusion.

5 Vertical diffusion

To measure vertical diffusion on large and small scales, we employ the root mean square $$z_{\mathrm{d},\mathrm{r}\mathrm{m}\mathrm{s}}=\sqrt{⟨z_{\mathrm{d}}^2⟩}$$(9)

and standard deviation $${\sigma }_{z_{\mathrm{d}}}=\sqrt{⟨z_{\mathrm{d}}^2⟩-{⟨z_{\mathrm{d}}⟩}^2}$$(10)

of the vertical positions of the dust particles z_d, that is the root mean square distances to the mid-plane and to the average vertical particle position, respectively. In other words, the former gives the dust scale height relative to the disk mid-plane, the latter with respect to the middle of the dust layer. Figure 3 shows the evolution of these two quantities after the dust has sedimented for 3-cm-sized dust in our models of only the vertical shear instability and only the streaming instability.

As illustrated in the top and bottom left panels of Figure 2, the vertical shear instability induces an almost threadlike, waveshaped dust layer (see also Flock et al. 2017; Flock et al. 2020; Dullemond et al. 2022). Accordingly, the time-averaged root mean square of the vertical dust positions is more than an order of magnitude larger than the mean standard deviation. In contrast, the streaming instability gives rise to a dust layer that possesses a Gaussian shape centred on the disk mid-plane (see, e.g., the top panel of Fig. 3 of Schäfer & Johansen 2022; Bai & Stone 2010b; Li & Youdin 2021). Thus, the root mean square and standard deviation of the vertical dust positions are comparable.

Youdin & Lithwick (2007) introduce an expression for the ratio of dust to gas scale height as a function of the dimensionless dust diffusion coefficient in the vertical dimension δ_z, $$\frac{H_{\mathrm{d}}}{H_{\mathrm{g}}}=\sqrt{\frac{{\delta }_z}{{\delta }_z+\mathrm{S}\mathrm{t}}},$$(11)

where H_d is the dust scale height. Rearranging this equation, we can compute a diffusion coefficient both from the root mean square z_d,rms, $${\delta }_z(z_{\mathrm{d},\mathrm{r}\mathrm{m}\mathrm{s}})=\frac{\mathrm{S}\mathrm{t}}{{\left(\frac{H_g}{z_{\mathrm{d},\mathrm{r}\mathrm{m}\mathrm{s}}}\right)}^2-1},$$(12)

and from the standard deviation σ_zd, $${\delta }_z({\sigma }_{z_{\mathrm{d}}})=\frac{\mathrm{S}\mathrm{t}}{{\left(\frac{H_{\mathrm{g}}}{{\sigma }_{z_{\mathrm{d}}}}\right)}^2-1},$$(13)

of the vertical particle positions. The resulting diffusion coefficients are depicted in the left panel of Figure 4 and listed in the third and fourth columns of Table 2 for simulations of all four of our scenarios, one each with a dust size of 3 mm and of 3 cm.

In line with what we discuss above and what can be seen in Fig. 2, we find the vertical shear instability in isolation to induce considerably stronger large-scale than small-scale diffusion. This discrepancy is considerable in our simulation of centimetresized dust, where the diffusion coefficient derived from the root mean square of the dust positions is equal to 10⁻³, but the one derived from the standard deviation to a minuscule 10⁻⁶. The difference is less significant for the millimetre-sized dust that is more tightly coupled to the gas, sediments more slowly, and is thus more affected by the (large-scale) turbulence caused by the vertical shear instability. Here, both the large-scale diffusion coefficient amounts to 3 × 10⁻³ and the small-scale one to 5 × 10⁻⁴.

The streaming instability alone, in contrast, in our model causes diffusion that is largely independent of scale. The diffusion coefficient computed from the root mean square vertical position of the centimetre-sized dust is marginally higher than the one obtained from the standard deviation, while for the millimetre-sized dust the two coefficients are in fact equal. All four diffusion coefficients are of the order of 10⁻⁵, with the two coefficients of the larger dust being greater by a factor of a few than those of the smaller dust. This is in agreement with both linear analysis and previous non-linear simulations showing that the streaming instability causes stronger turbulence and diffusion if the Stokes number of the dust is higher (Youdin & Goodman 2005; Johansen & Youdin 2007; Schreiber & Klahr 2018; SJB20; Yang & Zhu 2021; Baronett et al. 2024). Nonetheless, Schreiber & Klahr 2018 report a positive correlation also between diffusion coefficient and scale. We caution, though, that their simulations do not include the vertical stellar gravity.

Since the streaming instability is the dominant source of turbulence and diffusion in the dust layer of the scenario SIwhileVSI (SJB20), the diffusion coefficients in this scenario are comparable to the ones arising if only the streaming instability is considered. Nevertheless, for the centimetre-sized dust the diffusion coefficient calculated from the standard deviation of the vertical positions, and even more so the one derived from the root mean square, are greater in the SIwhileVSI simulation than in the simulation of the streaming instability only. We interpret this as the streaming instability sufficiently diffusing the dust for it to be susceptible to the large-scale diffusion induced by the vertical shear instability.

Similar to in our model of the vertical shear instability alone, in the scenario SIafterVSI small-scale diffusion is significantly weaker than large-scale diffusion. This is because in this scenario turbulence and diffusion are mainly driven by the vertical shear instability also in the dust layer (SJB20). Nonetheless, since the scenario SIafterVSI includes the drag exerted by the dust on the gas, the vertical shear instability is suppressed by an effective buoyancy resulting from the tendency of the dust to settle to the disk mid-plane (Lin 2019; SJB20). The diffusion coefficients of the millimetre-sized dust are thus reduced compared to when only the vertical shear instability is simulated. On the other hand, we find the streaming instability to contribute to the turbulence in the dust layer (SJB20), and indeed to compensate for the lack of small-scale diffusion caused by the vertical shear instability. This is evident both from the middle and bottom right panels of Fig. 2 and from the coefficient of small-scale diffusion of the centimetre-sized dust, that is to say the coefficient obtained from the standard deviation of the vertical dust positions, being similar to the value in the scenario SIwhileVSI.

Fig. 2 Ratio of dust ρd to gas (volume) density ρg as a function of radius r and height z, the latter in units of gas scale heights Hg (upper panels). Zoom-in on dust-to-gas density ratio ρd/ρg as a function of radius r and height z, with grid structure and dust particles being plotted as black lines and grey dots, respectively (lower panels). Simulations of the vertical shear instability only and of the scenario SIafterVSI with a dust size of 3 cm are depicted in the upper and left panel as well as lower and right panel, respectively. While the dust is stirred up to a marginally greater scale height with respect to the mid-plane in the former simulation, the dust layer is much thinner. Its width indeed amounts to less than a grid cell, while the layer in the latter simulation extends over multiple cells both in the radial and in the vertical dimension.

Fig. 3 Root mean square zd,rms (blue solid line) and standard deviation σzd (orange solid line) of the vertical dust particle positions as functions of time, considering only the time after dust settling. The mean of each quantity is shown as a dashed line. The left and right panel depict simulations of the vertical shear instability and the streaming instability with 3 cm-sized dust. The average root mean square is more than an order of magnitude greater than the standard deviation in former simulation, while they are similar in the latter simulation.

Fig. 4 Dimensionless coefficients of vertical dust diffusion δz calculated from the time-averaged root mean square zd,rms (circles) and standard deviation σzd (squares) of the vertical dust particle positions (left panel). Dimensionless radial dust diffusion coefficients δr computed considering all dust particles that are initially located within Δrinit = 0.1 au (circles) or Δrinit = 0.01 au (squares). In both panels, open and filled symbols represent simulations with a dust size of 3 mm and 3 cm, respectively. Simulations of only the vertical shear instability, only the streaming instability, the scenario SIwhileVSI, and the scenario SIafterVSI are depicted using blue, orange, green, and purple markers, respectively. We note that the vertical diffusion coefficients in the vertical shear instability simulation with the larger dust size fall outside of the ordinate range and are therefore plotted in an inset with a different ordinate.

6 Radial diffusion

To examine radial dust diffusion induced by the streaming instability and the vertical shear instability in isolation and in conjunction, we first need to distinguish the scales on which turbulence is the main source of diffusion from the scales on which it is primarily caused by the inward drift that arises even in the absence of turbulence. To this end, we track the separation of pairs of dust particles through time, specifically of all pairs among 10 000 particles randomly selected in every simulation. Figure 5 shows the evolution of the separation of pairs of centimetre-sized dust particles in simulations of the vertical shear instability and the streaming instability, respectively. It can be gathered that, if their initial separation is large, two dust particles approach each other over time. This is because the outer particle drifts more rapidly inwards than the inner particle since the radial drift speed scales with the dust Stokes number, which in turn is proportional to the dust size (see Eq. (8) and Sect. 4.2 of SJ22).

If it initially is less than ~0.1 au – this corresponds to 0.08 gas scale heights at the inner and 1.5 scale heights at the outer radial domain boundary – on the other hand, in the simulation of only the vertical shear instability the pair separation remains largely constant. Consistent with what can be seen from Fig. 2 and what we discuss in the previous section with respect to the vertical dimension, this demonstrates that the vertical shear instability induces minimal diffusion in the radial dimension. In contrast, the diffusion induced by the streaming instability in isolation leads to the pair separations gradually increasing with time, with this increase being largely independent of the initial separation.

To investigate the scale-dependence of radial diffusion, in what follows we compute and compare diffusion coefficients on scales of 0.02 au and of 0.004 au. As is evident from Fig. 5, the former is the largest scale on which the effect of the inward drift on the relative diffusion of the pairs of centimetre-sized dust particles in our simulation of the vertical shear instability is negligible – since turbulence-induced diffusion is stronger, the transition to drift-induced diffusion occurs on greater scales in all other simulations. Our choice of the smaller scale, on the other hand, is constrained by the simulation resolution, as the grid cell edge length at the highest level of refinement amounts to 3.125 × 10⁻³ au.

The dust diffusion coefficient in the radial dimension can be calculated as (Youdin & Lithwick 2007; Johansen & Youdin 2007) $$D_r=\frac{1}{2}\frac{\mathrm{\Delta }{\sigma }_r^2}{\mathrm{\Delta }t},$$(14)

where $\mathrm{\Delta }{\sigma }_r^2/\mathrm{\Delta }t$ denotes the temporal change of the standard deviation of the radial dust particles positions⁴, and be nondimensionalised as ${\delta }_r=D_r/\left(H_{\mathrm{g}}^2{\mathit{\Omega }}_{\mathrm{K}}\right)$. To obtain coefficients for the larger and the smaller scale established in the previous paragraph, all particles that are initially located within 0.02 au and 0.004 au, respectively, around the radial midpoint of our simulation domains at 55 au are taken into account. We perform a linear regression of the squared standard deviation of their radial coordinates as a function of time, only taking into consideration the time after the dust sedimentation has concluded. The two dimensionless diffusion coefficients derived in this manner for every simulation are shown in the right panel of Fig. 4 and listed in the fifth and sixth columns of Table 2.

As is indicated by Figure 5, both the streaming instability and the vertical shear instability, separately or jointly, induce radial dust diffusion that depends not or only marginally on scale. This is in agreement with what Schreiber & Klahr (2018) observe in their two-dimensional simulations of the streaming instability. Nevertheless, we note that the radial scales we consider differ only by a factor of five, and a dependence may emerge in models resolving a broader range of scales.

As noted above, in our model the centimetre-sized dust experiences negligible diffusion by the vertical shear instability the dimensionless diffusion coefficients in this case are of the order of 10⁻⁹ and in fact negative. The coefficients in the simulation including millimetre-sized dust, on the other hand, amount to 4 × 10⁻⁴, similar to the coefficient of the small-scale vertical diffusion. A consistent picture of dust diffusion caused by the vertical shear instability emerges, reflected in Fig. 2: the diffusion is strongest on large scales in the vertical dimension, causing an undulating dust layer. Within this layer, though, the diffusion is isotropic and considerably weaker, by an order of magnitude for millimetre-sized dust while negligible if the dust is centimetre-sized.

In our model of the streaming instability, on the other hand, whether or not diffusion is anisotropic is dependent on the dust size. The simulation of millimetre-sized dust yields diffusion coefficients of 10⁻⁵ both in the radial and the vertical dimension. In contrast, we obtain vertical coefficients of ~5 × 10⁻⁵ but radial coefficients of 3 × 10⁻⁴, larger by a factor of 6, from the simulation of centimetre-sized dust. This is consistent with the findings by Schreiber & Klahr (2018), who show that dust is in general more strongly diffused in the radial than in the vertical dimension if its Stokes number amounts to 0.1 (see their Figs. 10 and 26), but the diffusion is mostly isotropic if the Stokes number is equal to 0.01 (see their Fig. 20). (These Stokes numbers approximately correspond to the dust sizes in our simulations, as can be gathered from Eq. (8).) A similar trend in the relation between anisotropy and Stokes number emerges in the study by Li & Youdin (2021), although the ratio between radial and vertical diffusion coefficients amounts to 5 in both their simulations of dust with a Stokes number of 0.01 and with a Stokes number of 0.1 (see their Fig. 6). We note that both Schreiber & Klahr (2018) and Gerbig & Li (2023) find this anisotropy to vanish in dust overdensities where the dust-to-gas ratio is high. Additionally, Yang & Zhu (2021) and Baronett et al. (2024) find diffusion to be isotropic for the Stokes numbers we study, and to be stronger in the vertical than in the radial dimension if the Stokes numbers is of order unity.

All radial diffusion coefficients, both of the millimetre-sized and the centimetre-sized dust as well as both on the smaller and on the larger scale, are doubled to tripled in the scenario SIwhileVSI compared to our model of the streaming instability in isolation. This indicates that, even though the streaming instability is the main source of turbulence in the dust layer causes turbulence locally in dust overdensities (SJB20), the vertical shear instability contributes to the dust diffusion. Furthermore, the contribution is independent of dust size and scale. This is noticeable since this instability induces virtually no diffusion of the centimetre-sized dust when regarded alone. Similar to what we discuss with regards to vertical diffusion in this scenario, we interpret this as the streaming instability giving rise to smallerscale diffusion that is picked up on by larger-scale diffusion owing to the vertical shear instability.

In comparison with the scenario SIwhileVSI, in the scenario SIafterVSI the radial diffusion coefficients are again enhanced by a factor of 3–7. They amount to 10⁻⁴ for the millimetre-sized dust, less than when the vertical shear instability alone is considered but higher than when only the streaming instability is simulated. In the scenario SIafterVSI, the vertical shear instability is the dominant source of turbulence in the dust layer, but the streaming instability causes turbulence locally in dust overdensities (SJB20). It is thus plausible that by taking into account all dust particles in a certain radius range, we infer coefficients that represent an ‘average’ over the dust diffusion induced by the streaming instability inside overdensities and by the vertical shear instability outside of them.

Nonetheless, in our simulation of the scenario SIafterVSI with centimetre-sized dust, the radial diffusion coefficients are as large as ~2 × 10⁻³, slightly higher even than the large-scale diffusion coefficient in the vertical dimension. This implies that the streaming instability and the vertical shear instability jointly diffuse dust of this size in the radial dimension. As in the scenario SIwhileVSI, this can be explained by the streaming instability inducing a base level of diffusion that renders the dust susceptible to the large-scale diffusion by the vertical shear instability.

Fig. 5 Difference between average and initial radial separation Δr of dust particle pairs as a function of their initial separation in simulations of the vertical shear instability (left panel) and the streaming instability (right panel) with 3 cm-sized dust. We note that the scale of the abscissa is linear between 10−2 and −10−2 and logarithmic otherwise. If their initial separation is less than ~0.1 au (equivalent to 0.08 and 1.5 gas scale heights at the inner and outer radial boundaries), the distance between two particles remains roughly constant if only the vertical shear instability is simulated, while the streaming instability causes the distance to gradually increase with time (differently coloured solid lines). The dashed grey lines mark initial separations of 0.004 au and 0.2 au; we choose these scales to measure small-scale and large-scale radial diffusion induced by the two instabilities individually and jointly. If the initial separation exceeds ~0.1 au, on the other hand, two particles move closer together with time in both simulations. This is because the inward drift speed of the particles increases with the radius, and the separation of particles at different radii thus decreases.

Fig. 6 Dust-to-gas density ratio ρd/ρg as a function of radius r and height z in our two-dimensional simulation of the vertical shear instability with 3 mm-sized dust. The dust layer is similar in thickness, both with respect to mid-plane of the disk and to the mid-plane of the layer itself, to the layer of millimetre-sized dust in the model by F20 (see their Fig. D1).

7 Comparison with three-dimensional model of vertical shear instability

Since the inability of the vertical shear instability to diffuse centimetre-sized dust in the radial dimension as well as on small scales in the vertical dimension is unexpected in our eyes, we aim to validate of our model of this instability. To this end, in this section we compare it with the model of the instability presented in Flock et al. (2020, hereafter F20). We note that the dust included in the latter authors’ model is at most 1 mm in size, and a direct comparison is thus possible only with our simulation of dust with a size of 3 mm. – Since F20 consider a dust material density of 3 g cm⁻³ while we assume a material density of 1 g cm⁻³, the Stokes number of the dust is indeed identical in the two models if the properties of the gas are. Nonetheless, a verification of our simulation of millimetre-sized dust can be extended also to the simulation of centimetre-sized dust.

F20 investigate a three-dimensional radiation hydrodynamical simulation of an irradiated T Tauri star-disk system extending from 20 to 100 au. The grid resolution is uniform with around 70 cells per gas scale height. Half a million dust particles each of 1 mm and of 100 μm were embedded to study their motions. Since the drag exerted by the particles on the gas was neglected, the streaming instability is not active. Similar to us, these authors observe mixing of the dust owing to the vertical shear instability as well as radial drift. While at first the mm-sized dust is concentrated in a narrow layer similar to what is shown in Fig. 6, as the instability grows the dust layer evolves to be broader in both the radial and the vertical dimension.

Figure 6 depicts the dust layer in our simulation of solely the vertical shear instability and mm-sized dust. As is evident when comparing to Fig. D1 of F20, both the scale height relative to the disk mid-plane and the thickness of this layer is comparable in our and in their simulation.

This observation by eye is corroborated by the dust diffusion coefficients we infer. Figure 7 depicts and Table 2 lists the dimensionless coefficients measured at radii of 30 au and 55 au in the simulation by F20 as well as at r = 55 au in our simulation. These coefficients are computed in the same manner as the ones shown in Fig. 4, that is they represent large-scale and smallscale diffusion in the vertical and the radial dimension. Indeed, the two vertical diffusion coefficients (see the left panel of Fig. 7) calculated at r = 55 au are similar in our two-dimensional and their three-dimensional simulation. This concordance indicates that the results obtained from our model are robust, particularly with respect to the dust scale height, despite it being axisymmetric.

The coefficients of radial diffusion (right panel of Fig. 7), on the other hand, are smaller by almost an order of magnitude in the simulation by F20 than in our simulation. This potentially is due to a three-dimensional effect not captured in our twodimensional model. Nonetheless, the strength of radial diffusion is largely independent of scale both when simulating two and three dimensions.

At a radius of 30 au in the simulation performed by F20, a giant vortex is induced by the Rossby wave instability. Dust diffusion in this vortex is independent of both direction and scale, with the radial diffusion being stronger at this radius than at r = 55 au where its source is the vertical shear instability. Furthermore, the vertical diffusion coefficients evince that the dust layer possesses a Gaussian shape centred on the disk mid-plane rather than the typical wave-like shape caused by the vertical shear instability.

Fig. 7 Left panel: average root mean square zd,rms and standard deviation σzd of vertical particle positions. Right panel: radial dust diffusion coefficients Dr for Δrinit = 0.1 au (circles) and Δrinit = 0.01 au (squares). In both panels, blue, orange, and green symbols represent quantities at r = 55 au in our two-dimensional simulation of the vertical shear instability with a dust size of 3 mm as well as at r = 55 au and at r = 30 au for the dust with a size of 1 mm in the three-dimensional simulation conducted by F20. While the vertical diffusion coefficients at r = 55 au in the latter simulation and in ours are comparable, the radial coefficients are larger by almost an order of magnitude in our simulation. Nonetheless, vertical diffusion is similarly dependent on scale while radial diffusion is scale-independent both in the two- and the three-dimensional simulation. At r = 30 au, the Rossby wave instability induces a giant vortex in the simulation by F20, which gives rise to diffusion that is similar in strength in the radial and vertical dimensions as well as on small and large scales.

8 Discussion

8.1 Implications for modelling gas turbulence and dust diffusion in protoplanetary disks

Quantitative analyses of protoplanetary disk turbulence are most commonly based on the α-model formulated by Shakura & Sunyaev (1973). This model was devised to describe turbulence as an effective viscosity that facilitates stellar accretion by enabling angular momentum transport away from the star. Turbulence does not necessarily entail outward angular momentum transport, though – while the turbulence induced by the vertical shear instability indeed does (Nelson et al. 2013; Stoll & Kley 2014), the streaming instability on average causes gas to move outwards, not inwards (SJB20).

Furthermore, the α-model is not designed to capture anisotropy and scale-dependence. Accurately describing, for instance, the strength of anisotropic turbulence induced by the vertical shear instability (Nelson et al. 2013; Stoll et al. 2017; SJB20) thus requires reporting at least an α-parameter each for the radial and for the vertical dimension. We note that, moreover, two different definitions of the α-parameter exist. In studies of the vertical shear instability, it is most commonly defined based on the Reynolds stress to quantify angular momentum transport (e.g. Nelson et al. 2013; Stoll & Kley 2014, 2016; Stoll et al. 2017). Yet, it can also be defined as the square of the Mach number when considering a mixing length approach and assuming the eddy turnover timescale to be similar to the inverse of the orbital frequency.

Similar to measuring gas turbulence, translating from its strength to the strength of the dust diffusion it entails is challenging. Taking our model of the streaming instability as an example, both gas turbulence and diffusion of millimetre-sized dust are largely isotropic, in agreement with previous models (Johansen & Youdin 2007; Schreiber & Klahr 2018; SJB20). In contrast, the centimetre-sized is diffused more strongly in the radial than in the vertical dimension (see also Schreiber & Klahr 2018; Li & Youdin 2021). In addition to vertical diffusion being counteracted by dust sedimentation, Li & Youdin (2021) speculate that dust accumulation in clumps – which is strong for the combination of dust-to-gas surface density ratio and dust size in this case (SJ22) – affects both the radial and the vertical diffusion coefficient. On the one hand, the radial coefficient is enhanced since dust particles embedded in clumps radially drift more slowly than isolated particles. On the other hand, the vertical coefficient is reduced because clumps are concentrated in the disk mid-plane.

Finally, our study highlights that the shape of a dust layer in the vertical dimension cannot necessarily be expressed using a single scale height. The scale height relative to the mid-plane can be substantially greater than the one relative to the mid-plane of the dust layer, most notably by more than an order of magnitude in our simulation of the vertical shear instability including centimetre-sized dust. This instability gives rise to strong largescale dust diffusion resulting in a wave-shaped layer, but comparatively weak small-scale diffusion. Nonetheless, the streaming instability causes similarly strong small-scale and large-scale diffusion, resulting in the two scale heights being comparable.

Dust diffusion has been argued to impede the formation of planetesimals. Linear analyses including an α-model prescription of turbulence reveal reduced growth rates of the streaming instability (Umurhan et al. 2020; Chen & Lin 2020), and indeed planetesimal formation has been demonstrated to be inhibited by driven Kolmogorov-like turbulence in numerical simulations of the non-linear instability (Gole et al. 2020). Gerbig et al. (2020) and Klahr & Schreiber (2020) derive criteria for gravitational collapse that expand on the Roche criterion by requiring dust self-gravity to overcome not only the stellar tidal forces but also diffusion.

Nevertheless, instabilities like the streaming instability, the vertical shear instability, and the convective overstability are not only sources of diffusion, but also drive planetesimal formation either directly or by causing pressure bumps and vortices in which dust accumulates (SJB20; Raettig et al. 2021; Lehmann & Lin 2022; SJ22). Indeed, in SJ22 we find that dust concentration is stronger and planetesimal formation therefore possible for smaller dust sizes and lower dust-to-gas surface density ratios when simulating the vertical shear instability and the streaming instability in concert than when modelling the streaming instability alone. This is despite this work showing that diffusion is stronger when both instabilities are active.

8.2 Limitations

The most significant limitation of our study is the assumption of axisymmetry. While this assumption is justified when considering the linear regimes of the vertical shear instability (Nelson et al. 2013; Barker & Latter 2015) and the streaming instability (Youdin & Goodman 2005), the non-linear regimes exhibit deviations from this symmetry (Nelson et al. 2013; Stoll & Kley 2014; Kowalik et al. 2013). We address this shortcoming by comparing our two-dimensional model of the vertical shear instability to the three-dimensional one presented by F20. This comparison shows that our model reproduces the strength of vertical dust diffusion and the resulting dust scale heights, though not the strength of radial dust diffusion. Moreover, our model cannot replicate the formation of vortices owing to the vertical shear instability (Richard et al. 2016; Latter & Papaloizou 2018; F20), with these vortices affecting dust diffusion as noted in Sect. 7.

For simplicity, we employ an isothermal equation of state when simulating the vertical shear instability (and an adiabatic equation of state to quench it when simulating the streaming instability in isolation). As detailed above, the more sophisticated radiation hydrodynamic model applied by F20, which provides less ideal conditions for the vertical shear instability to develop, nonetheless yields similar results. However, even this model does not include the coupling between dust diffusion owing to the vertical shear instability and gas cooling via gas-dust collisions regulating the growth of the instability (Fukuhara et al. 2021; Fukuhara & Okuzumi 2024; Pfeil et al. 2023; Pfeil et al. 2024). Self-consistent models by Fukuhara & Okuzumi (2024) and Pfeil et al. (2024) indeed demonstrate that, for millimetre-sized and centimetre-sized as we consider in this work which is only weakly coupled to the gas, the vertical shear instability possibly cannot grow rapidly enough and induce sufficient diffusion to prevent run-away dust settling.

Further limitations of our study include not modelling other sources of gas turbulence and dust diffusion in protoplanetary disks besides the streaming instability and the vertical shear instability as well as considering only uniform dust sizes and neglecting dust evolution. Previous studies have explored the interplay between the magnetorotational instability and the vertical shear instability (Cui & Lin 2021; Cui & Bai 2022), between the magnetorotational instability and the streaming instability (Johansen et al. 2007, 2011; Balsara et al. 2009; Tilley et al. 2010; Yang et al. 2018), as well as between the convective overstability and the streaming instability (Raettig et al. 2021). Dust size distributions have been included in models of both the streaming instability (e.g., Bai & Stone 2010b; Schaffer et al. 2018, 2021; Krapp et al. 2019; Zhu & Yang 2021; Yang & Zhu 2021) and the vertical shear instability (e.g., Stoll & Kley 2016; Flock et al. 2017; Flock et al. 2020), with Pfeil et al. (2023); Pfeil et al. (2024) additionally simulating dust growth.

9 Conclusion

Building on our previous work (Schäfer et al. 2020; Schäfer & Johansen 2022), in this study we analyse gas turbulence and dust diffusion in two-dimensional, axisymmetric global models of the vertical shear instability and the streaming instability with a dust-to-gas surface density ratio of 2%. The focus of our study lies on diffusion and its the dependence on both direction and scale.

The vertical shear instability, on the one hand, is a source of anisotropic turbulence (Nelson et al. 2013; Stoll et al. 2017; Schäfer et al. 2020) and dust diffusion, both are stronger in the vertical than in the radial dimension. A central finding of our study is that, while vigorous large-scale vertical diffusion leads to the dust layer possessing the shape of a wave (Flock et al. 2017; Flock et al. 2020; Dullemond et al. 2022), the instability induces considerably weaker diffusion on small scales in the vertical dimension as well as in the radial dimension. We measure dimensionless diffusion coefficients of the order of 10⁻³ when considering the large-scale vertical diffusion, but otherwise coefficients of ~5 × 10⁻⁴ in a simulation of millimetre-sized dust and vanishingly low ones in a simulation of centimetre-sized dust. Therefore, the dust layer is too thin to be resolved in the latter simulation despite the application of mesh refinement based on dust concentration.

The streaming instability, on the other hand, diffuses dust in a scale-independent manner, with the dust layer thus exhibiting a Gaussian shape (Bai & Stone 2010b; Li & Youdin 2021; Schäfer & Johansen 2022). Both turbulence (Johansen & Youdin 2007; Schäfer et al. 2020; Yang & Zhu 2021; Baronett et al. 2024) and diffusion of millimetre-sized dust are further independent of direction, while centimetre-sized dust is diffused more strongly in the radial direction than in the vertical one (see also Schreiber & Klahr 2018; Li & Youdin 2021; Yang & Zhu 2021; Gerbig & Li 2023; Baronett et al. 2024). The dimensionless diffusion coefficients amount to ~10⁻⁵ for the smaller dust size, and to 5 × 10⁻⁵ in the radial dimension and to 4 × 10⁻⁴ in the vertical dimension for centimetre-sized dust.

We further consider models in which both instabilities are active. As we show in Schäfer et al. (2020), the streaming instability is the dominant source of turbulence in the dust layer in the scenario SIwhileVSI, but the vertical shear instability in the scenario SIafterVSI. While dust diffusion is thus overall stronger in the latter scenario, both scenarios exhibit contributions of both instabilities to this diffusion. Most notably, the streaming instability compensates for the lack of small-scale diffusion owing to the vertical shear instability. This way, the undulating mid-plane layer internally acquires a width that is set by the diffusion caused by the streaming instability. Whether planetesimals can form in such an undulating layer of high density will require future computer simulations that include both large vertical extents (for the vertical shear instability to thrive), good resolution of the layer width (for the streaming instability), and the inclusion of the self-gravity of the dust particle component.

Acknowledgements

We thank the anonymous referee for their constructive feedback that helped improve this paper. To analyse and visualise the simulations presented in this paper, the Python packages⁵ (Turk et al. 2011), Matplotlib⁶ (Hunter 2007), and NumPy⁷ (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 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 SCIENCE HPC Center at the University of Copenhagen. U.S. and A.J. are thankful for funding from the Danish National Research Foundation (DNRF Chair Grant DNRF159). A.J. further gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant 2019.0442), the Göran Gustafsson Foundation, and the Carlsberg Foundation (Semper Ardens: Advance grant FIRSTATMO). M.F. acknowledges support from the European Research Council (ERC), under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 757957).

Appendix A Gas kinetic energy power spectra

To further examine the direction- and scale-dependence of the turbulence induced by the streaming instability and the vertical shear instability in isolation and in conjunction, we compute power spectra of the gas kinetic energy. We note that these spectra are limited in their expressiveness because, firstly, the domains of the simulations we conducted are non-periodic and, secondly, the grid resolution is variable because of the application of mesh refinement. To remedy the latter, we superimpose a covering grid with the base resolution on the domain, interpolating to this resolution where necessary.

Fig. A.1 Power spectra of the gas specific kinetic energy eg,kin, non-dimensionalised using the square of the sound speed cs, as functions of the wavenumber $k=\sqrt{k_r^2+k_z^2}$. The top left, top right, bottom left, and bottom right panels, respectively, show simulations of the vertical shear instability and of the streaming instability as well as of the scenarios SIafterVSI and SIwhileVSI with a dust size of 3 cm. The kinetic energy is calculated taking into account different velocity components, with the resulting spectra being plotted using differently coloured lines. All spectra are averaged over 100 yr and normalised to the spectrum of Kolmogorov turbulence Eg,kin ∝ k−5/3. The spectra are similar in our simulations of only the vertical shear instability, of the scenario SIafterVSI, and of the scenario SIwhileVSI when considering scales greater than ~1 au, with the kinetic energy associated with the vertical velocity being higher by as much as an order of magnitude than the energy associated with the radial velocity. In comparison, the energies in the simulation of the streaming instability in isolation are largely equal for the two velocity components but generally lower; only on scales of ~1 au and smaller are they comparable to the energy in the radial motions owing to the vertical shear instability. At the same scales, the power spectra in the simulation of the scenario SIwhileVSI seem a mixture of the spectra in the simulations of the streaming instability and the vertical shear instability individually in terms of magnitude and direction-dependence of the kinetic energy.

In Figure A.1, we depict the kinetic energy power spectra obtained at the end of our simulations of only the vertical shear instability and the streaming instability as well as of the scenarios SIwhileVSI and SIafterVSI, all including dust with a size of 3 cm. The spectra are normalised to the spectrum of Kolmogorov turbulence e_g,kin ∝ k^−5/3, where e_g,kin is the gas specific kinetic energy and the wavenumber $k=\sqrt{k_r^2+k_z^2}$. This spectrum is characteristic of freely evolving three-dimensional turbulence, though, while two-dimensional turbulence exhibits an inverse energy cascade e_g,kin ∝ k^−5/3 from the injection scale towards larger scales and an enstrophy cascade – enstrophy can be defined as the square of the vorticity – with e_g,kin ∝ k⁻³ towards smaller scales (Kraichnan 1967; Lyra & Umurhan 2019).

Consistent with previous work (Nelson et al. 2013; Stoll et al. 2017; SJB20) and with what we discuss in Sect. 3, the vertical shear instability gives rise to anisotropic turbulence, with the kinetic energy associated with vertical motions being greater by up to an order of magnitude than that associated with radial motions on all scales. In comparison, the turbulence induced by the streaming instability is both weaker (SJB20) and mostly isotropic (Johansen & Youdin 2007; SJB20; Schreiber & Klahr 2018; Yang & Zhu 2021; Baronett et al. 2024). Only on scales of 1 au and smaller are the kinetic energies in the simulation of the streaming instability comparable to the energy in the radial motions in the simulation of the vertical shear instability.

These findings are in line with the vertical shear instability diffusing dust more strongly on large vertical scales than the streaming instability (see Sect. 5). On the other hand, they render the inefficacy of the vertical shear instability to diffuse centimetre-sized dust in the radial dimension (see Sect. 6) and on small scales in the vertical dimension even more striking.

The power spectra in the scenario SIafterVSI resemble those in the simulation of the vertical shear instability in isolation. That is, they do not reflect that, while the turbulence elsewhere is predominantly caused by the vertical shear instability, the streaming instability is the main source of turbulence locally in dust overdensities (SJB20). This is likely because the contribution of the kinetic energy in these overdensities to the power spectra is negligible, particularly since we apply a covering grid with the base resolution to obtain the spectra.

On scales larger than 1 au, the turbulence in the scenario SIwhileVSI as well exhibits power spectra comparable to the ones of the turbulence induced by the vertical shear instability. On smaller scales, though, the spectra appear to be a superposition of those in our simulation of this instability and in our simulation of the streaming instability. Indeed, we find the turbulence to be driven by the vertical shear instability away from the dust layer in this scenario, but mostly by the streaming instability in this layer (SJB20).

References

All Tables

All Figures

Fig. 1 Total gas kinetic energy Eg,kin in the region between r = 20 au and 30 au and restricted to one gas scale height off the mid-plane as a function of time t after the initialisation of the dust particles at td,init. Differently coloured lines represent simulations of the streaming instability and the vertical shear instability in isolation as well as of the scenarios SIwhileVSI and SIafterVSI, respectively, with the dust size being equal to 3 cm in all simulations. Like the model including only the vertical shear instability, the scenario SIafterVSI as well as the scenario SIwhileVSI at late times are characterised by a comparatively high total energy, with the energy of the vertical motions being substantially higher than that of the radial motions. In comparison, the total energy is less in the model of the streaming instability and in the scenario SIwhileVSI at early times, but largely equally distributed among the radial and vertical velocity components.

In the text

Fig. 2 Ratio of dust ρd to gas (volume) density ρg as a function of radius r and height z, the latter in units of gas scale heights Hg (upper panels). Zoom-in on dust-to-gas density ratio ρd/ρg as a function of radius r and height z, with grid structure and dust particles being plotted as black lines and grey dots, respectively (lower panels). Simulations of the vertical shear instability only and of the scenario SIafterVSI with a dust size of 3 cm are depicted in the upper and left panel as well as lower and right panel, respectively. While the dust is stirred up to a marginally greater scale height with respect to the mid-plane in the former simulation, the dust layer is much thinner. Its width indeed amounts to less than a grid cell, while the layer in the latter simulation extends over multiple cells both in the radial and in the vertical dimension.

In the text

Fig. 3 Root mean square zd,rms (blue solid line) and standard deviation σzd (orange solid line) of the vertical dust particle positions as functions of time, considering only the time after dust settling. The mean of each quantity is shown as a dashed line. The left and right panel depict simulations of the vertical shear instability and the streaming instability with 3 cm-sized dust. The average root mean square is more than an order of magnitude greater than the standard deviation in former simulation, while they are similar in the latter simulation.

In the text

Fig. 4 Dimensionless coefficients of vertical dust diffusion δz calculated from the time-averaged root mean square zd,rms (circles) and standard deviation σzd (squares) of the vertical dust particle positions (left panel). Dimensionless radial dust diffusion coefficients δr computed considering all dust particles that are initially located within Δrinit = 0.1 au (circles) or Δrinit = 0.01 au (squares). In both panels, open and filled symbols represent simulations with a dust size of 3 mm and 3 cm, respectively. Simulations of only the vertical shear instability, only the streaming instability, the scenario SIwhileVSI, and the scenario SIafterVSI are depicted using blue, orange, green, and purple markers, respectively. We note that the vertical diffusion coefficients in the vertical shear instability simulation with the larger dust size fall outside of the ordinate range and are therefore plotted in an inset with a different ordinate.

In the text

Fig. 5 Difference between average and initial radial separation Δr of dust particle pairs as a function of their initial separation in simulations of the vertical shear instability (left panel) and the streaming instability (right panel) with 3 cm-sized dust. We note that the scale of the abscissa is linear between 10−2 and −10−2 and logarithmic otherwise. If their initial separation is less than ~0.1 au (equivalent to 0.08 and 1.5 gas scale heights at the inner and outer radial boundaries), the distance between two particles remains roughly constant if only the vertical shear instability is simulated, while the streaming instability causes the distance to gradually increase with time (differently coloured solid lines). The dashed grey lines mark initial separations of 0.004 au and 0.2 au; we choose these scales to measure small-scale and large-scale radial diffusion induced by the two instabilities individually and jointly. If the initial separation exceeds ~0.1 au, on the other hand, two particles move closer together with time in both simulations. This is because the inward drift speed of the particles increases with the radius, and the separation of particles at different radii thus decreases.

In the text

	Fig. 6 Dust-to-gas density ratio ρd/ρg as a function of radius r and height z in our two-dimensional simulation of the vertical shear instability with 3 mm-sized dust. The dust layer is similar in thickness, both with respect to mid-plane of the disk and to the mid-plane of the layer itself, to the layer of millimetre-sized dust in the model by F20 (see their Fig. D1).
In the text

Fig. 7 Left panel: average root mean square zd,rms and standard deviation σzd of vertical particle positions. Right panel: radial dust diffusion coefficients Dr for Δrinit = 0.1 au (circles) and Δrinit = 0.01 au (squares). In both panels, blue, orange, and green symbols represent quantities at r = 55 au in our two-dimensional simulation of the vertical shear instability with a dust size of 3 mm as well as at r = 55 au and at r = 30 au for the dust with a size of 1 mm in the three-dimensional simulation conducted by F20. While the vertical diffusion coefficients at r = 55 au in the latter simulation and in ours are comparable, the radial coefficients are larger by almost an order of magnitude in our simulation. Nonetheless, vertical diffusion is similarly dependent on scale while radial diffusion is scale-independent both in the two- and the three-dimensional simulation. At r = 30 au, the Rossby wave instability induces a giant vortex in the simulation by F20, which gives rise to diffusion that is similar in strength in the radial and vertical dimensions as well as on small and large scales.

In the text

Fig. A.1 Power spectra of the gas specific kinetic energy eg,kin, non-dimensionalised using the square of the sound speed cs, as functions of the wavenumber $k=\sqrt{k_r^2+k_z^2}$. The top left, top right, bottom left, and bottom right panels, respectively, show simulations of the vertical shear instability and of the streaming instability as well as of the scenarios SIafterVSI and SIwhileVSI with a dust size of 3 cm. The kinetic energy is calculated taking into account different velocity components, with the resulting spectra being plotted using differently coloured lines. All spectra are averaged over 100 yr and normalised to the spectrum of Kolmogorov turbulence Eg,kin ∝ k−5/3. The spectra are similar in our simulations of only the vertical shear instability, of the scenario SIafterVSI, and of the scenario SIwhileVSI when considering scales greater than ~1 au, with the kinetic energy associated with the vertical velocity being higher by as much as an order of magnitude than the energy associated with the radial velocity. In comparison, the energies in the simulation of the streaming instability in isolation are largely equal for the two velocity components but generally lower; only on scales of ~1 au and smaller are they comparable to the energy in the radial motions owing to the vertical shear instability. At the same scales, the power spectra in the simulation of the scenario SIwhileVSI seem a mixture of the spectra in the simulations of the streaming instability and the vertical shear instability individually in terms of magnitude and direction-dependence of the kinetic energy.

In the text