The effect of the streaming instability on protoplanetary disc emission at millimetre wavelengths

In this paper, we investigate whether overdensity formation via streaming instability is consistent with recent multi-wavelength ALMA observations in the Lupus star forming region. We simulate the local action of streaming instability in 2D using the code ATHENA, and examine the radiative properties at mm wavelengths of the resulting clumpy dust distribution by focusing on two observable quantities: the optically thick fraction $ff$ (in ALMA band 6) and the spectral index $\alpha$ (in bands 3-7). By comparing the simulated distribution in the $ff-\alpha$ plane before and after the action of streaming instability, we observe that clump formation causes $ff$ to drop, because of the suppression of emission from grains that end up in optically thick clumps. $\alpha$, instead, can either increase or decline after the action of streaming instability; we use a simple toy model to demonstrate that this behaviour depends on the sizes of the grains whose emission is suppressed by being incorporated in optically thick clumps. In particular, the sign of evolution of $\alpha$ depends on whether grains near the opacity maximum at a few tenths of a mm end up in clumps. By comparing the simulation distributions before/after clump formation to the data distribution, we note that the action of streaming instability drives simulations towards the area of the plane where the data are located. We furthermore demonstrate that this behaviour is replicated in integrated disc models provided that the instability is operative over a region of the disc that contributes significantly to the total mm flux.


INTRODUCTION
According to core accretion theory, planets form through dust growth from initial m-sized grains up to the size of a planet (Safronov & Zvjagina 1969). Since a substantial change in the grain mass and size is required in this process, different physical processes are expected to happen during dust growth; therefore, the process is often divided into three main stages (e.g., Lissauer 1993;Papaloizou & Terquem 2006;Armitage 2007): grain growth, where the dust grains grow from m to cm, mainly through collision and sticking processes (Dominik & Tielens 1997;Birnstiel et al. 2012;Garaud et al. 2013;Dominik et al. 2016); planetesimal formation, from cm-sized dust grains to kmsized planetesimals; protoplanet formation, which leads to spherical objects of ∼ 10 3 km, that could be either rocky planets or gaseous planet cores (see, for example, Safronov & Zvjagina 1969;Kokubo & Ida 1996;Rafikov 2003;Ormel et al. 2010;Kobayashi et al. 2016;Kobayashi & Tanaka 2018).
The planetesimal formation stage represents a critical step in the growth process because the formation of km-sized objects is hampered by the so-called "metre-sized barrier" or "radial drift barrier" (Adachi et al. 1976;Weidenschilling 1977;Takeuchi & Lin 2002;Brauer et al. 2008;Pinte & Laibe 2014). The aerodynamic interaction between dust and gas in protoplanetary discs, in fact, causes ★ E-mail: ces204@ast.cam.ac.uk dust particles to lose angular momentum and to drift inwards. How rapidly the particles drift depends on the coupling between the gas and the dust, which is usually measured through the so-called Stokes number, s . Since s is proportional to the grain size, , we can relate the radial drift to , and it is possible to show that for grains characterised by ∼ cm − m drift can be as fast as drift ∼ 100 yr at 1 AU from the star (Armitage 2007;Birnstiel et al. 2016, for example). Consequently, unless other mechanisms interfere with the inward drift, the ∼cm-m sized grains are soon lost and they are no longer available to form planetesimals.
A potential way to avoid the metre-sized barrier is dust concentration in overdensities that rapidly collapse due to self gravity , and the streaming instability (Youdin & Goodman 2005) is a popular idea for how this might happen . Like radial drift, the streaming instability is caused by the interaction between dust and gas components within the disc, and, more precisely, by the backreaction of dust on the gas. When exerted by partially coupled particles ( s ∼ 1), the dust backreaction is strong enough that it can locally trigger on the disc midplane a powerful hydrodynamic instability (the streaming instability, Youdin & Goodman 2005), which, in appropriate conditions, promotes fast particle clumping, forming dust overdensity regions Carrera et al. 2015;Yang et al. 2017). Thus, streaming instability might be the key to solving the planetesimal formation issue: on one hand, the formed clumps would interact in a different way with the gas and this would slow down the drift velocity; on the other hand, local overdensities would facilitate the gravitational collapse and then planetesimal formation .
Clump formation through streaming instability can qualitatively be described as follows. The backreaction on the gas causes the gas component to orbit faster, and therefore it reduces the difference in azimuthal velocity between the gas and the dust, thus the headwind on the particles becomes weaker (Youdin & Goodman 2005;Squire & Hopkins 2020). If a dust overdensity were present, it would perturb the system equilibrium, causing a stronger backreaction, and therefore a reduced radial drift (which is determined by the difference in velocity between the two components). The system therefore faces an instability, as the initial dust overdensity increases more and more due to the new material drifting inward from outer orbits and stopping in the clump  described this effect as a "traffic jam"); the consequence is an exponential growth of clumps, whose density can become as high as 10 3 times the gas density, g (e.g, Bai & Stone 2010b,c). At the same time some particles are lost from the inner radius of the clumps, due to radial drift; therefore at some point a steady state is reached and the exponential growth saturates. The development of such exponentially growing instabilities is expected to be easier for marginally coupled particles , because their backreaction is stronger (the fastest radial drift is expected to occur for s ∼ 1).
The ability of streaming instability to form dust clumps depends on the local characteristics of the system, and in particular on the following three parameters: the grain Stokes number s , which regulates the gas-dust interaction; the pressure support, Π, which gives rise to the relative velocity and thus ultimately generates the streaming instability; and the local dust to gas mass ratio, , which must exceed a certain threshold in order to create dust clumps. From the theoretical point of view, several studies focused on exploring this parameter space, in order to find which combinations of ( s , Π, ) can trigger dust clumping via streaming instability. Bai & Stone (2010b,c) showed that, for a given grain size population, the stronger the pressure support, the higher is the critical to have streaming instability; whereas, for a given , smaller values of Π allow smaller particles to form clumps. More recently, Carrera et al. (2015) and Yang et al. (2017) studied in detail the critical required to trigger streaming instability for fixed s (single grain population) and Π, identifying the regions of the − s plane where streaming instability is expected to occur. Several studies also focused on the properties of streaming instability in multi-species simulations, both in the linear instability phase (Laibe & Price 2014;Drążkowska et al. 2016;Krapp et al. 2019;Zhu & Yang 2020) and in the non-linear phase (Bai & Stone 2010c;Schaffer et al. 2018). They found that, under certain conditions, the efficiency of linear growth can be increased/decreased when multiple grain species are considered with respect to the single grain models. As well as the physical properties of streaming instability, recent works studied in detail the influence of numerical parameters (size of the simulation box, resolution, etc.), algorithms and boundary conditions on streaming instability simulations (Yang & Johansen 2014;Li et al. 2018).
From an observational perspective, the recent Disk Substructures at High Angular Resolution Project (DSHARP) survey (Andrews et al. 2018), conducted with the Atacama Large Millimeter Array (ALMA), provided high resolution data for 20 protoplanetary discs, finding in most of them interesting substructures, such as rings, spirals, gaps, etc. These observations are allowing a more detailed analysis of the processes operating in protoplanetary discs, including planetesimal formation, thereby deepening our knowledge of these systems. Among other results, it has been found that the rings are marginally optically thick (Huang et al. 2018;Dullemond et al. 2018); this feature may be a mere coincidence, but it could also be the consequence of a common mechanism happening in rings, e.g. clump formation via streaming instability (Stammler et al. 2019).
Furthermore, Tazzari et al. (2020a,b) recently presented a new ALMA survey at 3 mm in Lupus star forming regions. By combining these results to the previous surveys at 0.9 mm and 1.3 mm (Ansdell et al. 2016(Ansdell et al. , 2018, they were able to study in detail the optical properties of observed systems and to extract important information about disc radii, the size-luminosity relation, the system optical depth, etc. Focusing on the spectral index and the optically thick fraction 1 , they noticed that the data concentrate in the region of spectral index ∼ 2.4 − 3 and with a range of values of the optically thick fraction, presenting a number of potential models that might explain the observed distribution in this plane. In this context, the aim of this paper is to investigate whether the formation of optically thick substructures through the action of streaming instability is consistent with the observation by Tazzari et al. (2020a,b). We therefore simulate systems where streaming instability is present (taking advantage of the broad parameter analysis performed in the studies mentioned above) and analyse their optical properties after the formation of dust clumps. We then compute the two observable quantities considered by Tazzari et al. (2020a,b) -the optically thick fraction and the spectral index -in order to perform a comparison between the data and the final distribution from our simulations.
The paper is organised as follows: in section 2 we briefly summarise the main parameters and equations describing streaming instability; in section 3 we describe the simulations performed and the method applied to analyse their optical properties; section 4 illustrates the process of clump formation in our simulations; the explanation of the consequences of the presence of substructures on optical properties is covered by section 5, where we also introduce a toy model which mimics the action of streaming instability, and we compare the local results to the disc-averaged observations of Tazzari et al. (2020b); in section 6 we define an integrated disc model and its optical properties are compared with the foregoing local results; finally, in section 7 we discuss the influence of the main model parameters and grain composition on the final results.

DUST-GAS INTERACTION AND STREAMING INSTABILITY
Protoplanetary discs contain both dusty and gaseous material, in an overall initial dust to gas ratio of about 1 : 100. Both the components orbit around the central star, but their azimuthal velocities are moderately different between each other; in fact, is the result of the balance of gravitational, pressure and centrifugal forces in the case of gas; whereas the dust is subject only to gravity and centrifugal force, and the absence of pressure force makes the dust azimuthal velocity slightly faster than that of the gas. This results in an aerodynamic interaction between the two components, exerted through drag forces, which (together with the consequent backreaction) plays the main role in triggering the streaming instability.
In this paper, we consider the Epstein drag force (valid for particles smaller than their main free path, Epstein 1924) where is the particle radius, g is the gas density and th is the gas thermal velocity. The ratio between the particle momentum and the drag force provides the interaction timescale between the gas and the dust grains of different sizes which can be conveniently expressed in units of the dynamical timescale Ω −1 obtaining the dimensionless stopping time, usually called Stokes number, which provides information about the coupling between gas and dust: s 1 means that the timescale for interaction is significantly smaller than the dynamical time, thus the particles, at first approximation, follow the gas motion; on the contrary, if s 1, the dust-gas coupling is weak and the particle motion is barely affected by the gas.
Equation (3) can also be reversed to compute the grain size corresponding to a given Stokes number where the th is written in terms of the isothermal sound speed s ∼ Ω k g , and g = 0 exp[− 2 /(2 2 g )] = Σ g /( √ 2 g ) is the gaussian profile assumed for the gas density (Σ g is the vertically integrated gas density).
The pressure gradient, which is the main cause of drag between dust and gas, is often characterised by the so-called pressure support parameter, defined as where g is the gas scale height, k is the keplerian velocity, s is the sound speed, and is related to the pressure gradient Finally, we introduce the dust to gas ratio (also called metallicity), defined as the ratio between the dust and the gas densities (Σ d and Σ g , respectively) which is required to be high in order to trigger streaming instability. Note that the requirement of high is just a local requirement, thus it is not in tension with the fact that the global dust to gas ratio is of the order of 1 : 100. Bai & Stone (2010a) implemented in A the dust-gas equations, written in a local reference frame located at a fixed radius 0 and rotating at angular velocity Ω (non-inertial reference frame). The continuity equation for the gas is where g and u g are the gas density and velocity, respectively; the Euler equation for the gas is g u g + ∇ g u g u g + I = where is the pressure, is the dust to gas ratio (for particle species ), , is the stopping time of particle species, and v k are the particle local mass density and velocity, respectively. The left-hand side includes the momentum time derivative and the effects due to advection and pressure gradient; the right-hand side includes the action of Coriolis force (first term), radial tidal momentum (second term), vertical gravity (third term) and backreaction (fourth term). The i-particle equation of motion is where the first term on the right-hand side is the effect of the pressure gradient, and all the other terms corresponds to those in Equation 9.
We have neglected the self-gravity of the dust in these simulations. Since, as we show later, the clumps formed in our simulations without self-gravity already contribute negligibly to the total flux of the disc, the inclusion of self-gravity, which acts to further enhance clumping, is not expected to have a significant effect on the observed properties.
We also neglected the effects of turbulence, as is common in streaming instability simulations (Liu & Ji 2020). Turbulence is generally expected to reduce particle clumping via streaming instability (Umurhan et al. 2020;Gole et al. 2020), therefore the usual noturbulence assumption corresponds to the most optimistic case for clumping via streaming instability.

METHODS
To understand whether the action of streaming instability is consistent with recent ALMA observations, we characterise the system observable quantities before and after the action of streaming instability. Our method can be described as a four-step process: (1) we perform hydrodynamics simulations of systems where streaming instability takes place; (2) we define a physical disc model, which allows us to translate the simulation results to physical systems; (3) we compute the radiative properties of these systems and (4) we compare two observable properties to those derived from ALMA observations by Tazzari et al. (2020a,b).

Numerical simulations
To simulate the action of streaming instability, we perform 2D hydrodynamics simulations of dust and gas using the hybrid code A , which treats the gas as a fluid on an Eulerian grid (Stone et al. 2008) and the dust as superparticles on that grid (Bai & Stone 2010a). We use a 2D (vertical and radial directions) shearing box approach, which allows to simulate a portion of the disc located at an arbitrary radius 0 and rotating at angular velocity Ω k ( 0 ).
In the following the simulation parameters are given in code units: the time unit is the dynamical time Ω −1 k and the length unit is the gas layer thickness g . The choice of the mass unit is arbitrary, because the equations describing gas and dust interaction (Equation 9, 10) depend on the dust to gas mass ratio, but are independent of the disc mass.
For all the 12 simulation we performed, we considered a box of  In each box we include tot = 7 · 10 5 particles, equally distributed in species = 28 particle species. The Stokes numbers of the different species are uniformly distributed in log space between a minimum and maximum Stokes number, min and max ; 6 of our simulations are characterised by min − max = 10 −4 −10 −1 , and 6 are characterised by min − max = 10 −3 − 1. The pressure support parameter is Π = 0.025 in all the simulations, while two values for the dust to gas mass ratio are considered ( = 0.02 and = 0.03).
All the simulations evolve for at least 1500 dynamical times, in some cases simulations had not reached a steady state in 1500 Ω −1 , thus the running time has been extended to 2000 Ω −1 (see Appendix A for more details).
We use reflecting boundary conditions in the vertical direction and periodic boundary conditions in the radial direction.
As initial conditions, all the particle species are distributed according to a Gaussian shaped distribution in the vertical direction, centred in the disc midplane, whose standard deviation (corresponding to the initial layer thickness) is d = 0.015 g . We use the parameter to specify the particle mass distribution per logarithmic particle radius bin or, equivalently, the number distribution ( ) where p and p are the particle mass and number, respectively. In our simulations we consider = 3, = 3.5, = 4 for each choice of and set of Stokes numbers.
In Table 1 we summurise the parameter used in our simulations.

Physical model
Studying the optical properties of simulated systems requires the definition of a disc model, allowing conversion of the dimensionless simulation results into physical units. As reference model, we consider the Minimum Mass Solar Nebula (MMSN) model by Chiang & Youdin (2010). Chiang & Youdin (2010) model is characterised by a flaring index = 2/7, thus the aspect ratio at the generic disc radius is where 0 is the box location (in the following we take 0 = 35 AU, unless otherwise specified) and the aspect ratio at 0 corresponds to the length unit unit = g ( 0 ). The time unit is defined as the inverse Keplerian velocity at the box location unit = Ω −1 k ( 0 ). Thus, it depends only on the box location and the star mass, which we assume to be the same as that of the Sun * = . We then choose the arbitrary mass unit by fixing the gas density at 1 AU and computing the gas density at box location 0 as where = 1.5 in Chiang & Youdin (2010) model. We do not fix a single value for the gas density, but we vary it in the range Σ g (1 AU) = 100 − 3000 g/cm 2 (cf 2200 g/cm 2 in Chiang & Youdin 2010 model), which corresponds to Σ g ( 0 ) = 0.5 − 14.5 g/cm 2 ; for each simulation, therefore, we obtain a population of discs characterised by different masses. For a disc characterised by inner and outer radii in = 0.1 AU and in = 70 AU, respectively, the chosen density distribution corresponds to discs of masses between 10 −3 M and 3 · 10 −2 M .
Finally, we define the temperature profile where the temperature at 1 AU is 0 = 120 K (unless stated otherwise) and the temperature profile is the same as that used in the review by Chiang & Youdin (2010), while the star luminosity is * = . The corresponding temperature at the box location (35 AU) for the chosen temperature profile is ( 0 ) ∼ 26 K.

Optical properties
To study the dust optical properties, we start by computing the opacity of dust grains; we use Birnstiel et al. (2018) dust code (which applies the Mie theory combined with optical datasets to compute the system optical properties), assuming (similarly to Tazzari et al. 2016) spherical compact grains, whose composition includes water, silicates (from Warren & Brandt 2008;Draine 2003), troilite (from Henning & Stognienko 1996) and organics (from Zubko et al. 1996). Note that the chosen composition is similar to that labelled as 'Zubko' in Birnstiel et al. (2018). 2 Through Birnstiel et al. (2018) code we obtain the dust material density p = 1.632 g/cm 3 , used to compute the dust grain size Then, we extract from Birnstiel et al. (2018) code the opacity corresponding to each simulated grain size and selected wavelength, obtaining the single grain opacity single ( ).
Since the real grain size distribution is expected to be continuous, each simulated represents a set of grain sizes between and + , thus we compute the opacity of each grain size as the average opacity Then we compute the size averaged opacity, which gives the system average response (see the Appendix B for further details) where ( ) is the mass of the dust grain. The main characteristics of avg as a function of the maximum grain size max are the presence of a steep increase when max ∼ /(2 ), followed by a decline towards larger max values.
We use the information on avg to obtain the optical depth where we assume the inclination angle to be = 0, and we use this information to obtain the specific intensity (assuming that the system is in local thermodynamical equilibrium) where is the Planck function evaluated at the local temperature ( ) (for the location of the considered box we obtain ( 0 ) ∼ 26 K). Finally, we compute the flux

Comparison with observations
After studying the systems' optical properties, we compare our simulations to the observations obtained by Tazzari et al. (2020a,b) in the Lupus star forming region. Therefore, we focus on the same plane used in Tazzari et al. (2020a,b) to study their data distribution: the optically thick fraction -spectral index plane. The optically thick fraction is defined as the ratio between the system flux and the flux that the system would emit if it was a black body where is the Planck function. Note that if a system were optically thick ( 1), it would have ∼ (see Equation 20) and therefore ∼ 1; on the contrary, an optically thin system ( 1) would be characterised by ∼ , thus ∼ . The spectral index, instead, measures the frequency dependence of the flux It is also useful to introduce the opacity index (see Appendix B for further information) which describes the variation of the opacity with frequency (assuming that avg ∝ ). ( max ) describes how avg varies with max ; we thus expect to have a peak for max ∼ /(2 ), i.e. where avg has a steep increase with grain size. Since the value of max corresponding to peak opacity scales with the wavelength, if we consider any two wavelengths 1 and 2 , so that 1 /(2 ) < max < 2 /(2 ), the opacity is much higher at wavelength 1 and hence the value of beta is large (see Appendix B for further details and plot regarding the spike in ).
In the limit of low optical depth ∝ , can be related to the opacity index which can be derived by using Equation 23 and Equation 24. In the Rayleigh-Jeans limit, = 2 + . Thus higher is associated with higher in the optically thin case; for optically thick emission in the Rayleigh-Jeans limit, = 2. We populate the − plane with the simulation results and, by adding on the same plane the data distribution by Tazzari et al. (2020a,b), we verify whether or not the simulations are consistent with observations. Since the data are available in ALMA bands 3 (100 GHz, 3.3 mm), 6 (230 GHz, 1.33 mm), and 7 (345 GHz, 0.88 mm), we compute the optical properties in these bands; specifically, we compute the optically thick fraction in band 6 B6 , while we use band 3 and band 7 to compute the spectral index B3−B7 (note that we need two values to compute , because it is defined as a derivative).

CLUMP FORMATION VIA STREAMING INSTABILITY
Considering, as an example, simulation T41-Z03-Q4 (see Table 1) in this section we study the properties of clumps formed in the simulated box due to the action of streaming instability.
In Figure 1 we show the vertically integrated density profiles of dust ( dust ( , ) is the 2D dust density, in vertical and radial directions) for a selection of 10 particle species from the 28 species used in the simulation. The three panels correspond to different timesteps of the system evolution: . The density profile of each particle species corresponds to a different colour in a rainbow palette, from the smallest particle species in red, to the biggest one in magenta. 3 The initial density profiles (left panel) are uniform for all the particle species, and clumps are gradually formed as the system evolves (central and right panels). By comparing the central and right panels, we notice that the biggest particle species (magenta, violet and blue lines) are already involved in clumps after 200 Ω −1 , while smaller species (cyan and green lines) require more time to participate in clumps; the smallest particle species (yellow and orange lines), instead, never participate in clumping. In fact, particles characterised by Stokes number close to 1 are the most affected by the drag forces, thus they are expected to clump via streaming instability more rapidly; while the smallest particles are highly coupled to the gas, thus they follow gas evolution and do not clump.

IMPACT OF STREAMING INSTABILITY ON OBSERVATIONS: LOCAL MODEL
In this section we study the observable properties of the simulated boxes undergoing streaming instability. In the following, unless spec- . Dust radial density profiles as a function of radius in the simulated box for simulation T41-Z03-Q4. The rainbow colour palette goes from red, which corresponds to the smallest particle species, to magenta, which corresponds to the biggest particle species. The three plots correspond to three different simulation timesteps: the initial condition = 0 Ω −1 (left panel), when all the particle species have an uniform density Σ d ∼ 0.003; an intermediate snapshot ified otherwise, we adopt the disc physical model described in section 3.2.

Distribution in the − plane before and after particle clumping
To analyse the effects on the optical properties caused by clump formation via streaming instability, we examine the optically thick fraction (Equation 22) and the spectral index (Equation 23) of the integrated emission from the box, both before and after the formation of clumps. Since the simulations are scale free (being characterised by the dimensionless numbers set out in Table 1), we can map a single simulation on to a variety of outcomes through considering a range of gas surface densities in the box; we considered 10 values for Σ g at the box location, logarithmically distributed between 0.5 g/cm 2 and 14.5 g/cm 2 .
In Figure 2 we show the distribution in the − plane for 3 simulations (T41-Z03-Q4, T30-Z02-Q4, T30-Z03-Q3; see Table 1). In all the three panels the blue squares represent the initial distribution of the selected systems, whereas the magenta diamonds correspond to the distribution of the 10 discs at the end of the simulation, i.e. after clumps have formed. The black arrows link together the initial and final condition of the same system. The intensity of the blue/magenta symbols indicates the gas surface density, with the darkest (lightest) shading indicating the highest (lowest) surface density values. For reference, the data distribution is shown in the same plane with the green stars; we notice that clump formation either drives the models towards the area occupied by the data (left panel), or keeps the models in the data area (central panel), or does not affect significantly the distribution (right panel). See section 5.2 and section 6 for a detailed comparison between simulations and data.
Focusing first on the initial conditions (blue squares), we notice that even if they are different in the three considered simulations (due to the difference in the parameter choices), they show some common features. In all cases, the system characterised by the highest gas surface density (darkest blue square) represents the highest value for and the lowest one for : at fixed metallicity this system has the highest dust surface density and hence optical depth and the value of is therefore close to 2 as expected for optically thick emission. Conversely, for lower gas surface density, the lower optical depth results in systems with lower and higher .
We consider first the evolution of when clumps form. Clumping involves depositing some of the grains in regions where the optical depth is higher than it was initially; if the clumps are optically thin then this does not affect the over-all flux produced by the box. However, if the clumps are optically thick then it means that some of the emission from grains in the uniform initial conditions is now hidden and therefore the over-all flux declines. Consequently clumping in general reduces , though in the case of very optically thin initial conditions (as in the right hand panel where the high maximum Stokes number and relatively low means that the emission is dominated by large, low opacity grains) may be hardly affected by clumping.
The effect on is more subtle and depends on the grain size (and hence optical properties) of the grains that are removed from the uniform background and deposited in optically thick clumps where their emission is concealed. Recalling that it is the largest grains that are deposited in the clumps, the direction of evolution of depends on whether these largest grains have higher or lower opacity index, , than the rest of the grain population.
We illustrate this effect via a 'toy model' (see also Appendix C for more details) in which we modify the initial conditions of a system to mimic clump formation by modifying the particles' distribution. Since in section 4 we observed that only particle species characterised by high Stokes numbers participate in clumping, we split the grains in two groups: half of the grains (the smallest ones) are left in the uniform background, whereas the density profile of the remaining from the highest density (Σ g ( 0 ) = 14.5 g/cm 2 , corresponding to disc ∼ 3 · 10 −2 M ), represented by the darkest blue/magenta markers, to the lowest density (Σ g ( 0 ) = 0.5 g/cm 2 , corresponding to disc ∼ 10 −3 M ), represented by the lightest blue/magenta markers. The green stars illustrate the data distribution by Tazzari et al. (2020a) . grains (the biggest ones) is modified as follows where indicates the -th species; Σ ,0 is the initial uniform density and represents the radial fraction of the 2D box in which we concentrate the dust = ( max − min ), where ∈ [0, 1]. Note that we are re-distributing only the biggest particles, thus Σ d will not be 0 outside the main peak, since the smallest particles are retained in the background; thus the total dust density for each clump model is given by In the left panel of Figure 3 we show the dust density profiles obtained by redistributing the initial uniform dust density (black line) of simulation T30-Z03-Q4, considering Σ = 0.5 g/cm 2 and 11 values for parameter , logarithmically distributed between 0.01 and 1 (the clump height increases from the lowest value in black to the highest value in pale pink). In the right panel of Figure 3 we show the distribution in the − plane (as before, we consider Σ g = 0.5 − 14.5 g/cm 2 ) obtained by computing the observable quantities for each clump width: the colours in this plot correspond to the dust distribution on the left with the same colour, while the black arrows link together the initial and final conditions of the same system. To explain the behaviour of the spectral index, we focus on the lowest density model (lower right in the plot) and the highest density model (upper left in the plot), which show opposite behaviours in .
We first recall that the spectral index (in the optically thin limit) can be related to the opacity index as follows thus we expect to increase (decrease) when increases (decreases). We therefore show in Figure 4 the opacity index as a function of the maximum grain size max (black line). The two panels correspond to the highest and lowest density models on the left and right respectively. In both the panels, two coloured areas highlight the corresponding to clumping (intense violet/blue area) and non-clumping (light blue/violet area) particle species in the simulation (separated at a Stokes number of s = 0.036, which corresponds to different grain size in the two panels). We notice that, in the lowest density case, the clumping grains are close to the opacity resonance and hence have higher values than the rest of the population. Once these grains go into the clumps (which are optically thicker than the uniform background), their contribution to the system's emission is downweighted. Therefore the overall decreases and thus we expect to decrease -as, indeed, happens in the lowest density model in the right panel of Figure 3. On the contrary, when we consider the opacity index for the highest density model (right panel in Figure 4), we observe that clumping species are here characterised by lower values than the non-clumping species. Thus we expect both and to increase after particle clumping; which is consistent with the behaviour of the highest density model in the right panel of Figure 3.
These behaviours can be related to Figure 2. 4 In the central panel, the high maximum Stokes number implies relatively large grains and  Figure 3. Left panel: re-distribution of the initial homogeneous dust density of simulation T30-Z03-Q4 (for the case with local gas density Σ g = 0.5 g/cm 2 ) to create artificial clumps characterised by various heights and widths, assuming that only half of the particle species clump (the biggest ones); a colour code is used, so that dust concentration increases moving from black to pale pink. Right panel: observable quantities for all the density distributions considered in the left panel (the colour of each diamond indicates at which dust density re-distribution that colour corresponds) and different values of Σ g .  the lowest density model in this paper (Σ g = 0.5 g/cm 2 ), while the plot in the left panel is obtained considering the highest density model in this paper (Σ g = 14.5 g/cm 2 ). Note that grain size distribution is assumed to be a power-law with q=4.
hence clumping increases the spectral index (though this effect is weaker at the lowest gas surface densities where the grain size is smaller at fixed Stokes number). In the left hand panel, conversely, the lower maximum Stokes number places the largest grains in the simulation close to the opacity resonance. Clumping therefore reduces in this case.

Comparison with data
In Figure 5 we show the distributions in the − plane for all the simulated systems and for the data distribution (green stars) obtained by Tazzari et al. (2020a,b).
We split the simulations into two groups characterised by the same Stokes number range: in the left panels we show the initial (upper panel) and final (lower panel) distribution for simulations  Table 1). The left panels illustrate the distribution in the observable plane for simulations characterised by max = 0.1, at the beginning (upper panel) and at the end (lower panel) of the system evolution. The two panels on the right show the same distributions, but for simulations characterised by max = 1. The values of the main parameters are indicated in the legend, according to the following rule: the squares correspond to = 3, the diamonds to = 3.5, the dots to = 4; in the left panel, the shades of red (blue) are used for systems characterised by = 0.02 ( = 0.03); in the right panel, the shades of pink (brown) are used for systems characterised by = 0.02 ( = 0.03). In all the panels, the data distribution is illustrated through the green stars.
characterised by s = 10 −4 − 10 −1 ; the right panels show the initial (upper panel) and final (lower panel) distribution for simulations characterised by s = 10 −3 − 1. In each panel, the different colours refer to different combinations of and , as indicated in the plot legend.
It is worth noticing that the initial distributions (upper panels) partially cover the data distribution. However, most cases characterised by max = 10 −1 have too high optically thick fractions in the highest density models, which are unable to cover the lower left area where data lie; secondly, the lowest density models present too high values for . Similarly, some of the max = 1 models present higher than those of the data (see, in particular, the yellow and pink models).
If we then compare the initial distributions to the final ones (lower panels), we note that the action of streaming instability pushes the simulated systems into the region occupied by data. Indeed, the optically thick fraction is lowered by clump formation, and this effect enables the high density models to re-cover the low data that previously were not matched by simulations. Moreover, the spectral index of the low density models is significantly reduced for cases with max = 0.1; this happens because in these models the clumping species are those characterised by values close to the resonance, thus by removing these grains from the background emission, the value must decrease. Therefore, for these cases, the action of streaming instability changes the system optical properties so that they become consistent with the data; thus the streaming instability can be considered a candidate to explain the data distribution.
We further underline that there are two simulations (T30-Z02-Q3 and T30-Z03-Q3) whose observable quantities barely evolve when clumps form. This behaviour is related to the fact that, in these cases, both the clumping and the not clumping species are relatively big (due to the particular combination of s and ), thus the creation of clumps does not alter significantly the system opacity.

IMPACT OF STREAMING INSTABILITY ON OBSERVATIONS: INTEGRATED MODEL
After studying the local observational properties of clumps, we define in this section an integrated disc model to study the global optical properties of systems undergoing streaming instability.
In the integrated disc model we assume azimuthal symmetry, and we define a disc characterised by the following gas density profile where Σ g (1AU) varies between 100 g/cm 2 and 3000 g/cm 2 (note that these values at 1 AU correspond to the values used in the local model at the box location) and we define in = 0.1 AU and out = 70 AU as inner and outer radii, respectively. The corresponding dust density profile can be obtained as Σ d = Σ g , where is chosen according to the considered simulation (see Table 1). In this model, the Stokes number is assumed to be constant across the disc (found to be a reasonable approximation in the Birnstiel et al. 2012 twopopulation model).
We divide the defined disc in rings = 100 'rings' from in = 0.1 AU to out = 70 AU, equally spaced on a linear scale. We then use the Σ g and g local values to map simulations into physical units. If the streaming instability is triggered in a particular ring, then the quantities calculated for that ring are those from the end state of the simulation; otherwise the initial conditions are used.
Once the dust density is obtained, both the opacity and the optical depth are computed at each ring following the same procedure as that used in the local model; in fact, each ring behaves as a single box. In each ring we compute the flux ring , then we sum over all the rings to obtain the total flux It is worth noting that for realistic emissivity profiles, the weighting by surface area ensures a relatively large contribution from the outer disc.
As final step, we compute the observable quantities, which can be simply obtained by using Equation 22 and Equation 23 where TOT is used instead of .
To explore different scenarios, we consider 4 different disc models: (1) all the disc is involved in streaming instability; (2) the inner disc is involved in streaming instability; (3) the outer disc is involved in streaming instability; (4) a ring in the disc is involved in streaming instability.
Following the method outlined above, we determine the distribution in the − plane for all the four integrated disc models. Since the flux is dominated by the outer disc, we found that in models (2) (where the instability occurs only within 35 AU) and (4) (where the instability occurs in a ring of width 10 AU located at 35 AU), the action of streaming instability hardly affects the flux and, therefore, the observable quantities; for the same reason, streaming instability modifies the observable quantities in models (1) and (3), whose behaviour is similar to each other. Therefore, in the following we show the results only for the 'outer disc' model, where streaming instability is assumed to take place for > SI and we choose SI = 35 AU. As in the local model, we compare the distribution in the − plane obtained from all the performed simulations to that obtained from the data, and we show the results in Figure 6. As in Figure 5, the upper (lower) panels correspond to the initial (final) distributions of system characterised by max = 0.1 and max = 1, on the left and on the right, respectively. The legend is the same as that used in Figure 5.
As in the local model, the initial conditions partially match the data distribution, in fact, the two distributions overlap in the area characterised by ( , ) ∼ (0.6, 2.5); nevertheless, the simulated systems are initially unable to cover the central area where most of the data are located ( , ) ∼ (0.4, 2.5). After the action of streaming instability, however, the optically thick fractions decreases, because, as previously explained, the surface density of particles in the optically thin background is reduced by clump formation, reducing the flux and allowing the final distribution to cover the data area that the initial conditions were unable to match. Note that the models do not cover few data located at ( , ) ∼ (2.4, 0.8); this can be explained by noticing that increases and decreases as the disc radius is reduced, as shown by Fig. 5 and Fig. 7 in Tazzari et al. (2020b). Since our integrated disc models are characterised by a constant disc size ( out = 70 AU), we expect the low-high-data to be reproduced by considering smaller discs (indeed, we verified that they are matched by models with out = 35 AU).
It is also worth noticing that in Tazzari et al. (2020a) the spectral index is correlated with radius, suggesting that the large spectral indices of transition discs can be explained if they have large radii for their masses. This result is consistent with the finding by Andrews et al. (2018) that transition discs are large for their flux; which is also reconcilable with the fact that they have a cavity in their mmemmision. Such a cavity in the inner disc, however, is not expected to modify significantly the results obtained in this section, as the outer disc is dominant in determining the observable consequences of the action of streaming instability.
We caution that the results might depend on the particular choice of the integrated disc model parameters, see section 7.2 for a detailed discussion.
Finally, we underline that if we compute the mm fluxes from our models, we find that they are broadly consistent with those observed, if we adopt the same temperature profile and disc sizes as those observed.

Influence of model parameters on the distribution (local model)
Although the physics of streaming instability depends only on the parameters chosen in the simulations, the observable quantities also depend on the specific disc model used to convert the simulation dimensionless code units to physical units: the choice of the mass  Table 1) applying the outer disc integrated model. The upper panels show the distribution in the plane for the initial condition, for simulation characterised by max = 0.1 and max = 1 in the upper left panel and upper right panel, respectively; the lower panels show the distribution after the system evolution (in the left panel for max = 0.1 and in the right panel for max = 1). The values of the main parameters are indicated in the legend (following the same rule as in Figure 5). In all the panels, the data distribution is illustrated through the green stars.
unit influences the grain size, hence the opacity; the local gas density depends on the choice of the box location and the chosen density profile; the temperature profile influences the Planck function. In Figure 7 we test the effect on the distribution in the − plane of changing the temperature, changing 0 to 60 K in blue, 120 K in purple and 240 K in green (corresponding to temperatures of 13 K, 26 K, 52 K at the box location.); 5 as an example, we consider simulation T41-Z03-Q4. The shades of each colour represent different local 5 Note that we only calculate the effect of changing temperature on the radiative properties and do not model the effect of changing Π(∝ 0.5 ) on the simulations; according to Bai & Stone (2010b), a factor two increase of gas density, logarithmically distributed from Σ g = 0.5 g/cm 2 to Σ g = 14.5 g/cm 2 (the darker the colour, the higher the density). For comparison, the green stars represent the data distribution.
Changing 0 , has no effect on the optically thick fraction, since it only depends on the dust density profile and the opacity (i.e. the grain sizes and their composition), that are independent of 0 (which only modifies the Planck function). Regarding the spectral index (Equation 23), we expect it to be unaffected by 0 if the emission is in the Rayleigh-Jeans limit (in which case = 2 + ); for low 0 , Π over the canonical value causes a modest increase in the dust to gas ratio required to trigger the streaming instability. simulation T41-Z03-Q4 considering different temperature scaling: 0 = 60 K in blue, 0 = 120 K in purple, 0 = 240 K in green (local temperatures of 13 K, 26 K, 52 K). The shades of colours correspond to different values for local Σ g (logarithmically distributed in the interval 0.5 − 14.5 g/cm 2 ). The green stars represent the data distribution by Tazzari et al. (2020a) . however, the emission is not in the Rayleigh-Jeans regime and this means that is lower for lower temperatures due to the fact that the higher frequency emission is beyond the Wien peak, resulting in a lower value of . We see from Figure 7 that the effect on the final distribution is relatively small and that the overall distribution is not considerably affected by the choice of 0 .

Influence of model parameters on the distribution (integrated model)
We also test the influence of the parameter choices in the case of the integrated model, focusing, as an example, on the 'outer disc' model. The main parameters to be studied are out and SI .
In the left panel of Figure 8, we plot three distributions obtained using the method described in section 6 for different disc outer radii: out = 50 AU (blue squares), out = 70 AU (purple squares), out = 100 AU (green squares). The value of SI is modified so that the portion of disc which is involved in streaming instability is always half of the total disc ( SI = out /2). The considered values for the gas density at 1 AU is the same regardless of the disc outer radius (Σ g (1AU) = 100−3000 g/cm 2 ). The three obtained distributions are very similar, and the only difference is that smaller discs are slightly optically thicker; in fact, the outer disc portion (which dominates the flux) in smaller discs is optically thicker than in a bigger disc (see Equation 31), therefore overall the disc is expected to be optically thicker.
In the right panel of Figure 8 we test the effect of changing SI when out = 70 AU is fixed. In particular, we consider SI = 15 AU (blue squares), SI = 35 AU (purple squares), SI = 55 AU (green squares). We notice that when SI , is lower, and thus more of the disc participates in the streaming instability, the optically thick fractions tend to be slightly lower and the bunching towards a smaller range of (see left hand panel of Figure 2) also becomes more pronounced.
Overall, the distributions shown in Figure 8 are not significantly influenced by the choice of either out or SI in the range considered. Therefore, we can state that the result obtained in Figure 6 is reasonably model independent, provided that the streaming instability is operating over a region of the disc that contributes significantly to the total emission at mm wavelengths.

Composition and porosity
In this paper we analysed the optical properties of the simulated systems by considering the composition for dust grains similar to that used in Tazzari et al. (2016) (which corresponds to the composition labelled as 'Zubko' in Birnstiel et al. 2018). However, it is important to discuss how the results would be affected by different compositional choices. We first note, by combining Equation 4 and Equation 19, that at fixed dust to gas ratio, the Stokes number depends on the maximum grain size, optical depth and opacity via: We found that, in our models, the streaming instability was successful in improving the match to the properties of observed discs on account of the changes in optical depth and spectral index effected when the maximum grain size is close to the opacity resonance. For the compositions adopted here, when the optical depth fraction is several tens of per cent (as in observed systems) and is close to the opacity resonance, then the Stokes number is in the range where the streaming instability is strong. We can estimate the effect of using, for example, the DSHARP composition (see the case labelled as 'default' in Birnstiel et al. 2018). From the top panel in Figure 10 of Birnstiel et al. (2018), we observe that the DSHARP opacity is lower than the 'Zubko' one (approximately by a factor of 20) while the peak (middle panel) is located at higher values (approximately by a factor of 4). This means that for the DSHARP composition the Stokes number corresponding to the resonance at the same optical depth would be around a factor 5 less than what we have simulated in this paper. It remains to be seen whether the streaming instability is sufficiently strong at these lower Stokes numbers to have a significant effect on disc radiative properties.
Likewise we note that if, instead, porous grains are employed, the peak in associated with the opacity resonance virtually disappears (see middle panel of Figure 10 in Birnstiel et al. 2018) and thus the streaming instability would have little impact on mm emission from discs. In the absence of the streaming instability, Tazzari et al. (2020b) showed that the data in Lupus requires lower values than can be produced by porous grains. Thus since we find that the streaming instability does little to the radiative properties of discs in this case, we are led to disfavour the hypothesis of porous grains.

CONCLUSIONS
In this paper we simulated the action of streaming instability by performing 12 2D shearing box simulations with multiple grain sizes using the A code. By comparing the results from our simulations to observations in Lupus by Tazzari et al. (2020a,b), we found that the action of streaming instability is overall consistent with the integrated emission from discs at mm wavelengths.
Analysis of dust density profiles after clump formation shows that the largest particles participate significantly in clumping, whereas the smallest ones remain almost uniformly in the background. Streaming instability thus affects the emission properties of discs in cases where  obtained using different outer disc radii: out = 50 AU (blue squares), out = 70 AU (purple squares), out = 100 AU (green squares); in all the three cases we rescale SI so that in all models half of the disc is involved in streaming instability ( SI = out /2). The range of gas density at 1 AU is the same for all the three cases: Σ g (1AU) = 100 − 3000 g/cm 2 . The right panel shows three distributions obtained by fixing the outer radius out = 70 AU and changing the radius dividing the area undergoing streaming instability from the area where streaming instability does not take place: blue squares correspond to SI = 15 AU; purple squares to SI = 35 AU; green squares to SI = 55 AU. In both the plots, the colours of the squares vary their intensity according to the density of the corresponding disc (the higher is the disc density, the darker is the colour).
grains towards the upper end of the grain size distribution exhibit a steep dependence of opacity, and its wavelength dependence, on grain size. Such a steep dependence is associated, in the case of compact grains, with the opacity resonance at a grain size around a few tenths of a mm.
We explored the effect of the streaming instability on the location of models in the − (optical depth fraction versus spectral index) plane. The instability always reduces the optical depth fraction, because the effect of clumping is to reduce the emission from those grains that are translated from optically thin to optically thick regions (in agreement with the finding by Stammler et al. 2019 that the streaming instability reduces the system's optical depth). The spectral index can evolve in either direction depending on whether the maximum grain size is above or below the size corresponding to the opacity resonance. The net effect is to drive modeled disc systems towards a relatively narrow range of spectral indices (around ∼ 2.5) which agrees well with the observed distribution of discs in the − plane (see Figure 5).
Finally we remark that although we have investigated the specific scenario of clump generation by streaming instability, our results are likely to apply, at least qualitatively, to any situation where clumping predominantly affects the largest grains. Clumping is likely strongest for large particles in any effective clumping mechanism; as is the case for trapping by spiral arms in self-gravitating discs (as long as s 1, Booth & Clarke 2016), pressure maxima created by planets, and vortices. Our primary result (that clumping reduces the optically thick fraction and brings systems to a relatively narrow range of spectral index values) is therefore likely to be of general applicability.

APPENDIX A: CONVERGENCE
Since we are interested in the changes of optical properties before and after the action of streaming instability, we must ensure that at the end of our the running time our simulations have converged. We therefore computed for each simulation the observable quantities as a function of time, and we considered the simulation as converged when and stopped to increase/decrease. In Figure A1 we show, as an example, the variation of the spectral index with time for simulation T30-Z03-Q4: in the first 800 dynamical times increases considerably, then it stabilises on a nearly constant value. We performed this convergence test for all the simulations, finding that the  Figure A1. Spectral index as a function of time for simulation T30-Z03-Q4, considering the standard disc model with local gas density Σ g = 14.5 g/cm 2 . most of simulations have converged after 1500 Ω −1 , apart from a few simulations (simulations T41-Z02-Q4, T41-Z03-Q4, T30-Z02-Q3, T30-Z03-Q3) which required a longer evolution to reach convergence (2000 Ω −1 ).

APPENDIX B: OPACITY
Considering as an example simulation T41-Z03-Q4, we explain here the details of the method applied to compute the size averaged opacity for our simulations. The method can be split in three phases: (1) we compute the single grain opacity, using Birnstiel et al. (2018) code; (2) we interpolate over the grain size interval that each simulated grain represents (Equation 17); (3) we compute the size averaged opacity (Equation 18). To test whether our computation works properly we compute the opacity index as a function of the maximum grain size max and we compare it to the same quantity obtained with Birnstiel et al. (2018) code.
In the left panel of Figure B1, the magenta line shows obtained by considering a simulation characterised by 7 particle species and without operating the grain interpolation. The thick grey line show the result obtained from Birnstiel et al. (2018) code (as well as in the other two panels). This computation clearly unable to return the correct opacity, as the low number of species and the absence of grain size interpolation causes strong oscillations. By increasing the number of particle species to 28 (central panel), we note a considerable improvement in the computation of opacity, since the majority of oscillations disappear; nevertheless, at high max , seems to be underestimated. Once we introduce also the size grain interpolation (right panel), we are finally able to recover precisely the correct shape of − max . This analysis allowed us to choose the correct number of particle species to insert in our simulation (28 species) and to verify that Equation 17 is needed in order to compute the opacity properly.
In Figure B2 we applied the method outlined above to the cases = 3 (magenta line), = 3.5 (cyan line) and = 4 (violet line). The thin lines refer to our computation, the thick translucent ones refer to the corresponding result obtained through Birnstiel et al. (2018) code. Even though for small max , the three cases behave similarly,  Figure B1. Opacity index as a function of the maximum grain size for different models, obtained from simulation T41-Z03-Q4. The grey line is the result obtained from Birnstiel et al. (2018) code; the magenta line is the result of our computation. The left panel shows the result when we consider the single grain opacity and 7 particle species; the middle panel shows the result when we interpolate the grain distribution using 28 particle species; the right panel shows the result when we consider that each grain represents a set of grains and we extrapolate the distribution to include the low mass grains.  Figure B2. Opacity index as a function of the maximum grain size, for different grain size distributions: q=3 (magenta), q=3.5 (cyan), q=4 (violet). The thick translucent lines corresponds to the values obtained through Birnstiel et al. (2018) code; the thin solid lines corresponds to the values computed using our sets of dust grain distributions. they differ significantly among each other for the peak height and the value for big max . In particular, for big grains, the opacity variation with frequency is steeper and steeper as increases.

APPENDIX C: TOY MODEL
In section 5.1 we used a toy model to illustrate how the optical properties of clumping particles affect the increase/decrease in . Here, we show that the toy model is effective in reproducing the distribution in the − plane of simulated systems. We compute for each simulation the level of clumping for each particle species the mean density ratio (i.e. we compute the dust density for each particle species and average that over all the particles), then we use the result to identify the clump area as the area where the final mean density ratio is higher than the initial one, and finally we obtain by using Equation 28. The result is shown in Figure C1.
We then use the values of computed above to mimic clump formation of each particle species 6 through Equation 28 and Equation 27 and we compute the corresponding distribution in the − plane. We test the toy model by applying it to the initial conditions of simulations plotted in Figure 2 and we show the results in Figure C2 (the panels, from left to right, correspond to the panels in Figure 2). By comparing Figure C2 to Figure 2, we can notice that the behaviour of the observable quantities in the toy model is similar to that obtained through the actual simulations; this confirms that both the definition used to obtain in Figure C1 and the method of mimic clump formation by redistributing the particles are effective in reproducing the observable quantities. This paper has been typeset from a T E X/L A T E X file prepared by the author.