Streaming instability in a global patch simulation of protoplanetary disks

In the recent years, sub/mm observations of protoplanetary disks have discovered an incredible diversity of substructures in the dust emission. An important result was the finding that dust grains of mm size are embedded in very thin dusty disks. This implies that the dust mass fraction in the midplane becomes comparable to the gas, increasing the importance of the interaction between the two components there. We address this problem by means of numerical 2.5D simulations in order to study the gas and dust interaction in fully global stratified disks. To this purpose, we employ the recently developed dust grain module in the PLUTO code. Our model focuses on a typical T Tauri disk model, simulating a short patch of the disk at 10 au which includes grains of constant Stokes number of $St=0.01$ and $St=0.1$, corresponding to grains with sizes of 0.9 cm and 0.9 mm, respectively, for the given disk model. By injecting a constant pebble flux at the outer domain, the system reaches a quasi steady state of turbulence and dust concentrations driven by the streaming instability. For our given setup and using resolutions up to 2500 cells per scale height we resolve the streaming instability, leading to local dust clumping and concentrations. Our results show dust density values of around 10-100 times the gas density with a steady state pebble flux between $3.5 \times 10^{-4}$ and $2.5 \times 10^{-3} M_{\rm Earth}/\mathit{year}$ for the models with $\mathit{St}=0.01$ and $\mathit{St}=0.1$. The grain size and pebble flux for model $\mathit{St}=0.01$ compares well with dust evolution models of the first million years of disk evolution. For those grains the scatter opacity dominates the extinction coefficient at mm wavelengths. These types of global dust and gas simulations are a promising tool for studies of the gas and dust evolution at pressure bumps in protoplanetary disks.


Introduction
The interaction between the gas and the dust is a crucial step in the process of planet formation. Once grains collide in protoplanetary disk, they are able to stick and grow which leads them to decouple from the gas motion (Safronov 1972;Whipple 1972;Adachi et al. 1976;Weidenschilling 1977;Wetherill & Stewart 1993). Once grains settle to the midplane, they concentrate speeding up the dust coagulation and growth even further (Weidenschilling 1997;Stepinski & Valageas 1997;Weidenschilling 2000;Laibe et al. 2008;Brauer et al. 2008;Birnstiel et al. 2011). Depending on the radial profiles of temperature and density, these dense dust layers are prone to other types of instabilities, especially the gas and dust drag instabilities Zhuravlev 2020). The streaming instability (SI henceforth), one important subclass of these drag instabilities, has been studies extensively in the recent two decades, in particular for its role in explaining planetesimal formation (Youdin & Goodman 2005;Bai & Stone 2010;Yang & Johansen 2014;Carrera et al. 2015;Yang et al. 2017;Schreiber & Klahr 2018;Carrera et al. 2020). The growth rate of the SI depends on the local dust-gas-mass ratio and can reach values of around the orbital timescale for dust-to-gas mass ratios close to unity Pan & Yu 2020). Recent works focused on the SI with multigrain species (Laibe & Price 2014;Krapp et al. 2019;Zhu & Yang 2020;Paardekooper et al. 2020) demonstrating that its growth rate is in general reduced when considering multiple grain sizes. Lagrangian and fluid methods were used in the past and both have their advantages and disadvantages. As Lagrangian methods introduce a fixed number of grains, in stratified disk models this means that they can only resolve a certain height of the dust disk and one has to ensure for a good sampling to supress the noise level (Cadiou et al. 2019). On the other hand, they allow to follow individual grain motions becoming particularly suited to study larger grains which decouple from the gas motion.
Recent studies emphasized again the importance of the level of gas turbulence to determine where the SI can operate (Jaupart & Laibe 2020;Umurhan et al. 2020). So far most of the simulations have been performed in local box simulations, and only recently first simulations appeared using global unstratified simulations (Kowalik et al. 2013;Mignone et al. 2019) confirming the main characteristics of the SI. Very recently Schäfer et al. (2020) investigated the interplay between the vertical shear instability and the SI Article number, page 1 of 14 arXiv:2103.15146v1 [astro-ph.EP] 28 Mar 2021 in stratified global models, demonstrating the importance of the large scale gas motions for the dust concentrations.
In this work we propose a framework to study the SI in global stratified disk simulations with more realistic conditions of the radial pressure gradient profile and the pebble flux. For the first time to the extent of our knowledge, the SI is investigated in a high resolution stratified global simulations using spherical geometry. The challenge is twofold: first, the SI requires resolutions of several hundreds cells per gas scale height and second, the finite extent in the radial direction introduces a time limit to study the SI as the grains radially drift through the domain. In Section 2 we explain the numerical method and the disk setup, while in section 3 we present the results and compare them to local box simulation, emphasizing the role of the pebble flux compared to the classical total dust-to-gas mass ratio. Finally we test our results with constraints from typical T Tauri star disk systems and dust evolution models and give estimates on the optical depth at mm wavelengths. We present the discussion and conclusion in section 4 and 5.

Methods and disk setup
In order to setup our model we follow the work of Nakagawa et al. (1986) which prescribes the dust and gas velocities (v and V , respectively) using cylindrical geometry (R, φ, Z) as where ρ and ρ d denote, respectively, the gas and dust density, Ω K = GM/R 3 is the Keplerian frequency with the gravity constant G, M is the mass of the star. We define the factor with the dust to gas mass ratio = ρ d /ρ and the dimensionless Stokes number St = t s Ω K with t s being the stopping time. Likewise, we introduce where P is the (gas) pressure. We note here the factor of 1/2 which was also adopted by . With this, the pure gas azimuthal velocity translates to v φ = (1 − η)RΩ. 1 We assume a local isothermal equation of state with the pressure defined by P = c 2 s ρ with c s (R) being the speed of sound. The scale height H of the gas is defined as H = c s /Ω K with radial dependence 1 We note that some investigators define η without the factor 1/2, see also Takeuchi & Lin (2002).
where P and Q are the radial profile exponents for the density and temperature T ∼ c 2 s . The initial profile of the gas density in the R − Z plane is set by where r = √ R 2 + z 2 is the spherical radius. The dust density ρ d is defined similarly to Eq. (8) with the dust scale height H d : Finally, the initial total dust to gas mass ratio is defined through the ratio of the vertically integrated surface densities, that is, Σ d,0 /Σ 0 , related to the midplane densities as for the gas and for the dust using Σ d,0 and H d,0 respectively. The parameter for the model are summarized in Table 1.

Numerical configuration
Our computations are performed in spherical geometry using the HLL Riemann solver and the 2 nd -order Runge Kutta integration in time to advanced conserved variables. Piecewise linear reconstruction with the MC limiter has been used. The gravity force is handled by adding the vector force on the right hand side of the momentum equation for both gas and particles. The radial boundary conditions for the hydro variables are zero-gradient without allowing for material to enter the domain. This is obtained by setting the normal velocity component to zero in case the velocity in the active domain is pointing inward. At the meridional boundary, ghost zones are filled by extrapolating the exponential profile of the gas density while the remaining variables are set to have zero-gradient. Similarly to the radial boundaries, gas is not allowed to enter the domain. Dust particles are advanced in time using the exponential midpoint method and the cloud-in-cell (CIC) weighting scheme is applied to determine the gas values at the grain position (Mignone et al. 2019). Gas and dust are coupled by mutual feedback terms accounting for drag force effects. As particles are stored locally on each processor, we reach best parallel performance when we use a decomposition X : 2 for the r : θ domain. For more information about the method we refer to our previous work (Mignone et al. 2019). We note that the resolution was chosen to resolve the grains drag and therefore the grid size should fulfill to resolve the grains stopping length. For our setup with a fixed Stokes number this becomes ∆x < StH. Further we have to resolve the SI for the regimes of St = 0.1 and St = 0.01 and we adopt a grid resolutions similar to that used by Yang et al. (2017) with around ∼ 1000 cells per H. The aspect of resolution is discussed in detail in Section 4.2.
We use a uniform grid spacing in the radial and in the θ direction. The total number of grid cells and domain extent are given in Table 2.

Lagrangian particle setup
In order to sample the dust density we introduce individual particles to the simulation. Those particles represent a swarm of grains. Similar to the work by Yang et al. (2017), we start the simulation assuming a dust scale height which is reduced compared to the gas scale height. This increases locally the dust to gas mass ratio and enhances the growth rate of the SI.
We present our results using particle sample runs with approximately 80 particle per cells for model ST1 and 5 particles per cell for model ST2. For our models we determined that a sampling of around 5 cells or more is needed for consistent results. More details on the sampling and a benchmark can be found in appendix A and appendix B.
The total dust mass in the domain is Using a total number of grains N tot = 1.04×10 7 , we sample the dust mass in our domain with particle swarms with masses of 1.902 × 10 −7 Earth masses. Once the mass of the particle is fixed, the number of particles in each cell can be determined with N cell = (ρ d ∆V )/m g . A vertical profile of the detailed sampling over height is shown in the appendix A.
The initial profile of the gas and dust density in the R-Z plane is shown in Fig. 1. As seen in the contour plot, the particle method can only resolve a certain vertical extent which depends on the total number of grains. The initial gas and dust velocities are plotted in the bottom panel of Fig. 1, following the equilibrium solutions, Eq. (1), (2), (3) and Eq. (4).
In this equilibrium solution at the midplane, the grains slowly drift radially inward while the gas slowly drifts radially outward.

Damping zones
To prevent numerical effects from the radial domain boundary we need to implement buffer zone in which the vari-ables of density and velocity are relaxed back to their initial value. To implement those, we apply a wave killing zones close to the radial inner and outer boundary to reduce the interaction with the boundary. Similar buffer zones were already tested and implemented in our previous work Mignone et al. (2019). In this zones, the variables are relaxed to the initial equilibrium values, where Q(x x x, t) represents either the radial or azimuthal fluid velocity, Q 0 (x) is the corresponding equilibrium value, T = 1.0 and is a tapering function with w = 0.05. A second important point is to prevent the dust dragging the gas material out. As we constantly inject new grains at the radial outer zone we also have to replenish the gas material. To study the SI in a quasi steady state configuration, we implement a density relaxation which relaxes the loss of gas mass to the initial value. For each grid cell we apply setting the parameter T = 1.0. In appendix C we show that the influence of parameter T in this regime of the density relaxation remains small.

Injection and destruction of particles
At the outer boundary, new particles must be constantly injected as the radial drift empties these buffer zones. Refilling has to be done carefully as in this model the vertical settling quickly changes the structure in the buffer zones.
To prevent that the solution is quickly shifted away from the equilibrium solution we damp the vertical velocity as long as the particles remain in the outer buffer zone. More specifically, at every timestep we set in the radial outer buffer zone, using χ = 10 −5 . We found that this efficiently prevents the particle from settling already to the midplane before they have entered the active domain. In addition, we inject only dust when the dust density drops below a certain factor f , which is ρ d < f ρ d,0 with f = 0.5 for St = 0.1 particles and f = 0.25 for St = 0.01. In this way, we are able to resupply the disk with a constant pebble flux without any accumulation at the buffer zone edges. Particles are constantly injected in r ∈ [10.6, 10.7] au as they radially drift inward. For r < 10.6 au they start to settle eventually triggering the SI. At around 10.5 au and inwards, the structure and the dynamics of the dust and gas remains self-similar. Particles are removed from the computational domain once they cross the inner boundary at 9.3 au. To avoid dust accumulation at the inner buffer zone we reduce the azimuthal velocity through at every timestep in the inner buffer zone using χ = 10 −8 . This reduces their rotational velocity thus increasing their radial drift to avoid any concentration close to the buffer zone edge at 9.4 au.

Particle size
In our models we fix the Stokes number of the grains which makes a comparison to previous works easier. As the gas density in our domain does not vary much, we can determine approximately the dust grain sizes. For our disk setup we are in the Epstein regime and we can determine the size using where, for the grain density, we employ ρ grain = 2.7 g cm −3 . Using the midplane density at 10 au we determine the grain sizes of 0.88 mm for model ST2 and 8.8 mm for model ST1.

Results
After a few orbits, the dust settles to the midplane and radially drifts further inwards. At the same time, the SI is triggered leading to dust concentration and clumping. Dust grains are resupplied in the outer buffer zones, enabling a constant pebble flux. In the following we investigate for the dust concentrations and the dust scale height in our models.

Dust concentration and streaming instability
In the following we analyse the maximum dust density at the center of the domain at 10 au in a small radial patch of 0.1 au. Results, plotted in Fig. 2, show that after roughly 10 orbits the dust concentration reaches up to hundred times the initial value. Grains with St = 0.1 (top panel in Fig. 2) show a strong concentration, attaining a (temporally-averaged) maximal dust to gas mass ratios of max = 24 ± 14. The average concentration level at the midplane saturates at the time-averaged value of = 4 with a large scatter. Grains with St = 0.01 (bottom panel in Fig. 2) show slightly lower concentrations with (temporallyaveraged) maximum concentrations of around 10 ± 3 of the dust to gas mass ratio. The spatial averaged midplane value of the dust to gas mass ratio is 2.8. Fig. 3 shows the evolution of the vertically-integrated dust surface density for both models. After 10 orbits, the SI transforms the smooth surface density into dust fragments of low and high dense filaments. Fig. 3, top, shows narrow, close to horizontal stripes which indicate the fast inward radial drift for model ST1. At ∼ 80 and ∼ 100 orbits, two large dust accumulation become visible in the surface density which leads to a reduction of the radial drift. This is to be expected since large dust clumps can shield each other from the gas headwind. In these large dust clumps, the maximum dust concentration can lead to dust to gas mass ratios of ∼ 100, see Fig. 2, top, close to 90 orbits. In the bottom panel of Fig. 3 we display the temporal evolution of the dust surface density for model ST2. Here the SI also leads to overdense structures in the surface density. Due to the slower inward radial drift, the dust clumps show a broader structure over time. The concentrations in the surface density remain on a similar level as model ST1 although we do not observe the large dust accumulations. Fig. 4 shows snapshots of the dust density after ∼ 100 orbits for both models in the meridional plane. The dust layer is very thin, with a vertical extent of only 0.01 au and attaining dust densities of around 10 −11 g cm −3 . Fig. 4 (top panel) shows dust clumps for model ST1 on top of the narrow dust layer that remains concentrated in a region corresponding to roughly 1% of the gas scale height. The turbulent structures have sizes of around 1/10 of the gas scale height in radius, while there is a sharp density contrast along the vertical direction. Likewise we show (in the bottom panel of the same figure) the snapshot for model ST2. Here, the turbulent structures look much finer compared to model ST1, also with a smoother density contrast in the vertical direction.

Vertical dust scale height and effective α
To calculate the dust scale height we follow first the approach by Yang et al. (2017) and determine the standard deviation of the vertical position of the grains with z = R cos θ. Fig. 5 shows the evolution of the dust scale height over time which remains between 0.2 to 0.3 % of the gas scale height for both models.
To verify the scale height determination technique we follow another approach outlined in Flock et al. (2020) and plot the averaged vertical dust density profile (see Fig. 6). The plot provides the time averaged profile from 20 orbits until the end of the simulation and spatially averaged at 10 au ± 0.15au. The dust density is normalized by the gas density (which remains effectively flat in the vertical direction with only 5 per mill deviation). Both profiles fit best with a value H p /H of ∼ 0.003, particularly inside the first 0.01 au from the midplane. Above 0.01 au from the midplane the profile becomes more shallow, probably because of sudden bursts of dust clumps, as it is seen in Fig. 1

. By
Article number, page 5 of 14 A&A proofs: manuscript no. main  determining the dust scale height we can effectively determine the level of turbulence which would be equivalent to produce such a profile.
Following the calculation from our previous work Flock et al. (2020) and Dubrulle et al. (1995) we can estimate the turbulent diffusivity α with assuming a Schmidt number Sc of unity. Inserting the values for the Stokes numbers we derive an effective α of about 10 −6 for model ST1 and 10 −7 for model ST2.

Comparison to previous local box simulations
In what follows, we compare our simulation results with previous stratified local box models of the SI. Two param- eters are important to characterise the evolution of the SI, namely, the Stokes number and the value of σ d = Σ d /(ρ 0 ηr) which can be understood as an average dust-to-gas mass ratio at the midplane layer. Most of the dust mass is concentrated in a regions of 1% of the gas scale height around the midplane, which motivates the need to quantify the SI using different parameters such as done by Sekiya & Onishi (2018) (2017) found a dust vertical scale height for St = 0.01 grains and σ d = 1 of about H p = 0.014, roughly 3 times higher than in our models. Carrera et al. (2015) presented a model with σ d = 0.5 and his particle scale height was around H p = 0.005, very similar than the value we found. Sekiya & Onishi (2018) found that the strength of the SI scales with σ d , a result also predicted from the analytical works of ; Pan & Yu (2020) who demonstrated that the growth rate depends on the dust to gas mass ratio. As our models adopt a lower value of σ d , it might be that the strength of the SI is reduced. Yang et al. (2017) found long lasting dust concentrations appearing after hundreds to thousands of orbits while Sekiya & Onishi (2018) noticed that such concentrations appear for values of σ d ≥ 1. Model ST1 showed two events of secondary dust concentration reaching values of 100 times the gas density which reduced the radial drift of the dust clump, albeit this concentration was not enough to reach the critical Roche density, see appendix D. On the other hand, model ST2 showed no major dust concentration. We also note that such dust accumulations have been observed on timescales of several hundreds of orbits (Yang et al. 2017). The radial drift in combination with our limited radial domain extent does not allow us to trace individual grains for such a long time (we cannot follow individual grains for more than around 10 orbits in model ST1).

Dust Clumping
We thus conclude that our simulation results show similar values of dust clumping as observed in local box simulations. We find a smaller value of the dust scale height and no secondary long lasting dust accumulation, both possibly because of our choice of σ d .
Finally we performed a local box simulation with the same setup as presented in Yang et al. (2017) in Appendix E, to compare our new numerical method to previous models. We report that the results of dust concentration and particle vertical mixing in our local box runs are very similar as found in previous works.

Resolving the streaming and other instabilities
Both grid resolution and particle sampling are important to correctly represent and resolve the dust and gas interactions.
Grid resolution is crucial in order to resolve the fastest growing modes of the SI.  pointed out that the wavenumber of the fastest growing mode for the SI follows roughly kηr ∼ 1/St. In our model ST1, assuming k ∼ 1/∆x and η ∼ (H/R) 2 = 0.0049 we obtain H 2 /(R∆r) ∼ 45, therefore allowing us to resolve wavenumbers (ideally) up to 45 using 640 zones per scale height. Here the fasting growing mode corresponding to kηr ∼ 10 is well resolved. Model ST2 can capture wavenumbers kηr up to 180 and so also resolves the fastest growing mode corresponding to kηr ∼ 100. Yang et al. (2017) showed SI operating with St = 0.01 with a range of resolutions from kηr = 32 up to 256 while the strongest concentrations appeared when using resolutions close to the fastest growing wavelength.
Accurate particle sampling is fundamental to capture correctly the dust feedback and to resolve the dust distribution (Mignone et al. 2019). Sekiya & Onishi (2018) adopt a particles over grid cells ratio N par /N cell ≈ 1.39 while Yang et al. (2017) used a ratio of unity. In our models we employ N par /N cell ≈ 79.3 for model ST1 and N par /N cell ≈ 9.9 for model ST2, both much larger than the previous models. In the appendix of Yang et al. (2017) he compared results obtained using N par /N cell = 10 and N par /N cell = 1 without finding any significant difference.
From this perspective, our models provide the necessary grid resolution and particle sampling in order to capture the basic SI properties as well as to resolve for the dust feedback.
Another type of instability which might have an important effect is the vertical shear SI (Ishitsu et al. 2009) which is driven by the vertical gradient of the velocity shear between the dust and the gas. A recent work by Lin (2021) emphasizes the role of this instability in determining the vertical scale height of the dust. However, we point out that scales of the order of 10 −3 H have to be resolved to capture the vertical shear SI. Future high resolution simulations reaching ten thousand of cells per H are needed to verify the importance of new types of instabilities, like the settling instability or the vertical shear streaming instability.

Gas transport
Without damping, the gas surface density is quickly reduced creating strong radial pressure gradients which affected the gas surface density structure. Owing to the limited domain extent, we could not investigate this interesting effect. Specially due to the small vertical domain this effect is enhanced as there is no resupply of gas material from the upper layers. For this models we applied the gas density relaxation to study in more detail the SI in quasi steady state. Future simulations should include a much larger vertical extent to examine the effect of the dust drag onto the gas, particularly at regions where dust is expected to accumulate, such as the water ice line, where the dust drag can become very important for the gas motion (Gárate et al. 2020).

Regions of planetesimal formation
Over the recent years, several works have shown that the generation of planetesimal via the SI in a smooth disk profile remains difficult. First, a large amount of solid material (Z > 0.02), larger than the typical ISM value, has to be provided in a single dust species of certain Stokes number (Yang et al. 2017;Johansen et al. 2014), leaving a relative narrow range of parameter (Umurhan et al. 2020;Chen & Lin 2020), even more narrow when including the effect of turbulence (Jaupart & Laibe 2020). However such favorable conditions -a low amount of turbulence and a narrow range in mass distribution of large grains -is not expected from dust evolution models (Brauer et al. 2008;Birnstiel et al. 2011). Another difficulty arises in multi-grain simulations including different grain sizes, which showed the reduction of the SI grow regime (Krapp et al. 2019).
On the other hand, the locations of pressure maxima in the disk remain plausible regions for planetesimal formation, owing to the large amount of dust concentrations and the conditions for the SI to operate are favourable (Auffinger & Laibe 2018; Abod et al. 2019;Carrera et al. 2020).

The importance of the pebble flux
With this work we intend to emphasize the importance of the radial flux of pebbles rather than adopting the (more common) total gas-to-dust mass ratio when modeling the evolution of the SI. The pebble flux controls the transport of the solid material in the disk, it can show us where dust grains get concentrated and trapped and it is important for the accretion of solid material onto planets (Ormel & Liu 2018) which requires the understanding of their vertical distribution . The radial flux of pebbles can be determined using dust growth and evolution models evolution (Birnstiel et al. 2012;Takeuchi & Lin 2002;Drążkowska et al. 2016;Drazkowska et al. 2021).
Using the pebble flux simulator 2 we determine the pebble flux over time using the same disk profile employed in our simulations. The method and further references of the tool can be found in Drazkowska et al. (2021). The results are shown in Fig. 7. The maximum pebble flux is reached at around 10 4 years and it matches the value obtained in our model ST2, approximatelyṀ d ∼ 3.5 × 10 −4 M Earth /year (see Fig. 7 top). The pebble flux and grain sizes in model ST1 lie above the predicted values from the dust evolution models.
We then conclude that our model ST2 using (St = 0.01, σ d = 0.36) presents realistic conditions of the dust amount when compared to models of dust evolution. However such initial conditions are not favorable for secondary dust clumping events by the SI (Sekiya & Onishi 2018) needed to account for planetesimal formation (Yang et al. 2017). This is an important aspect which should be investigated in more detail in forthcoming simulations of the SI.

Total optical depth at mm wavelengths
An important question remains whether the dust layer observed in protoplanetary disks are optical thin or thick at a given wavelength. For this, we calculate the opacity of the two grain sizes at the wavelength of λ = 1.3 mm corresponding to the ALMA Band 6 observations. To calculate the opacity we use the optool 3 which is using the DIANA dust properties (Toon & Ackerman 1981;Woitke et al. 2016) and including the distribution of hollow spheres method (Min et al. 2005) to calculate the dust opacity. For the specific settings we use amorphous pyroxene (70% Mg) with a mass fraction of 87 % and 13 % of amorphous carbon (Zubko et al. 1996;Preibisch et al. 1993) and a water ice mantel with a mass fraction of 20 % and a porosity of 20%. We calculate the opacity for the two grain sizes 2 Pebble predictor tool on Zenodo 3 https://github.com/cdominik/optool/ using a narrow-size bin (0.8mm to 1mm) and (0.8cm to 1cm) for model ST2 and model ST1 respectively. The corresponding absorption and scatter opacity at λ = 1.3 mm are κ abs = 3.009 cm 2 /g, κ scat = 21.957 cm 2 /g and for model ST1 these are κ abs = 0.635 cm 2 /g and κ scat = 1.166 cm 2 /g. In Fig. 8 we plot the radial profile of the total optical depth τ = Σ d κ calculated for both models using the vertical integrated dust density.
The profiles show that for model ST1, the total optical depth remains around unity, while the optical depth from pure absorption opacity remains mostly optically thin. For model ST2, the grain size is closer to the corresponding wavelength. Here the optical depth is larger and it remains mostly above unity. The scatter opacity is much larger for this grains which leads to a total optical depth of around 10. More and more observations of protoplanetary disks at mm wavelength confirm the important effect of scattering (Sierra & Lizano 2020). Also grain sizes of around mm size are consistent with the observations (Carrasco-González et al. 2019).
Overall, the variations in τ abs caused by the SI fluctuate between 0.4 to 4 (thus a factor of ∼ 10) in model ST2 and between 0.02 and 2 (a factor ∼ 100) for model ST1. We point out again that these structures are on spatial scales of tens of H, which translates to scales of 0.1 au at the distance of 10 au from the star. Current radio interferometer capabilities of ALMA reach a spatial resolution of 5 au for the dust emission at mm wavelength in the most nearby star disk systems.

Conclusions
In this work, we have presented a new generation of models to investigate the dust and gas drag instabilities in global stratified simulations of protoplanetary disks. High resolution, 2D global hydrodynamical simulations have been performed, including the dust back-reaction on the gas modeling the conditions of a protoplanetary disk around a one solar mass star. Our numerical method is based on the hybrid fluid-particle framework recently developed by Mignone et al. (2019), where the dust component is modeled by Lagrangian particles. We adopt 2D spherical ge- ometry covering the meridional domain (r : θ) with a grid resolution up to 1280 cells per gas scale height to resolve for the streaming instability. The dust grains are modeled with a constant Stokes number of St = 0.1 and St = 0.01 which corresponds to grain sizes of 880 micron and 8.8 mm, respectively, at 10 au. The dust grains radially drift through the domain and they undergo streaming instability, leading to the formation of large dust concentrations. By resupplying dust grains at the outer radial domain, we reach a quasi-steady state of pebbles flux and operating streaming instability. Our main results may be summarized as follows: -The streaming instability leads to dust clumping, with maximum values between 10 to 100 in terms of the dust to gas mass ratio. The average dust to gas mass ratio at the midplane remains between 2 and 4. -For St = 0.1 we observe the appearance of large dust clumping reaching dust to gas mass ratios above 100 which can effectively reduce the radial drift as grains shield itself from the gas drag. We finally wish to emphasize the important role of the pebble flux when determining the amount of dust in simulations of the streaming instability. The maximum pebble flux in the disk is reached during the first million years of disk evolution. T-Tauri star disk models with a pebbles flux of around 3.5 × 10 −4 Earth masses per year and Stokes numbers 0.01 St 0.1 are closed to what is expected from dust evolution models. For this range of parameters, σ d remains below unity making secondary dust clumping for planetesimal formation difficult (Sekiya & Onishi 2018). This novel class of global dust and gas simulations constitute a promising tool for forthcoming studies targeting gas and dust evolution in protoplanetary disks, especially for situations where the density and pressure are strongly changing (such as at pressure maxima).

A. Dust sampling
In our models, the dust density is sampled by individual particles. In Fig. .1 we show the number of particles per cell along the vertical direction for both models ST1 and ST2 at 10au. The black and blue dotted lines shows the initial dust density profile and the theoretical sampling. The dust density can only be sampled until a given height, where one has at least one particle per cell, as indicated with the red dotted line in Fig. .1.

B. Benchmark -Dust sampling
Here we investigate the effect of the particle sampling. For this we perform a series of runs, based on model ST1, using different numbers of particles, ranging from 0.5 up to 160 particles per cell. Fig. .2 shows the dust scale height and the maximum dust concentration max at 10 au over time.
The results indicate that a sampling of 5 particles per cell or more is enough to show converging results. Below this value sudden dust concentrations occur which also trigger larger dust scale heights, possibly due to Kelvin-Helmholtz Type instabilities.

C. Gas damping
As we relax the gas density in our domain we have to verify that the gas damping does not strongly affect the nonlinear evolution of the streaming instability. To this end, we perform a test run setting the damping factor T = 10.0 in Eq. 13 for the gas damping, leading to a reduced gas damping rate. The results are summarized in Fig. .3. The lower gas damping does not strongly effect the main results.

D. Roche Density
A common approach to investigate whether the SI could produce directly planetesimals through dust clumping which then collapse due to the self-gravity, is by determin-