Gas accretion damped by dust back-reaction at the snow line

Context. The water snow line divides dry and icy solid material in protoplanetary disks. It has been thought to signiﬁcantly affect planet formation at all stages. If dry particles break up more easily than icy ones, then the snow line causes a trafﬁc jam because small grains drift inward at lower speeds than larger pebbles. Aims. We aim to evaluate the effect of high dust concentrations around the snow line onto the gas dynamics. Methods. Using numerical simulations, we modeled the global radial evolution of an axisymmetric protoplanetary disk. Our model includes particle growth, the evaporation and recondensation of water, and the back-reaction of dust onto the gas. The model takes into account the vertical distribution of dust particles. Results. We ﬁnd that the dust back-reaction can stop and even reverse the net ﬂux of gas outside the snow line, decreasing the gas accretion rate onto the star to under 50% of its initial value. At the same time, the dust accumulates at the snow line, reaching dust-to-gas ratios of (cid:15) (cid:38) 0 . 8, and it delivers large amounts of water vapor towards the inner disk as the icy particles cross the snowline. However, the accumulation of dust at the snow line and the decrease in the gas accretion rate only take place if the global dust-to-gas ratio is high ( ε 0 (cid:38) 0 . 03), the viscous turbulence is low ( α ν (cid:46) 10 − 3 ), the disk is large enough ( r c (cid:38) 100au), and only during the early phases of the disk evolution ( t (cid:46) 1 Myr). Otherwise the dust back-reaction fails to perturb the gas motion.


Introduction
Protoplanetary disks are composed of gas and dust.In the classical picture, a gas disk evolves through viscous evolution driven by outward transport of angular momentum (Lynden-Bell & Pringle 1974) and orbits at a sub-Keplerian speed due to its own pressure support.In another view, dust particles couple to the gas motion according to their size (Nakagawa et al. 1986;Takeuchi & Lin 2002) and small grains quickly follow the motion of the gas, while large boulders are decoupled from it.The mid-sized grains, or pebbles, feel a strong headwind that causes them to drift towards the gas pressure maximum (Whipple 1972;Weidenschilling 1977) which, in a typical disk, is directed towards the star.
At interstellar dust-to-gas ratios of 1%, the force exerted by the dust into the gas is mostly negligible.Yet in regions such as dead zones (Kretke et al. 2009;Pinilla et al. 2016), the outer edges of gaps carved by planets (Dipierro & Laibe 2017;Kanagawa et al. 2018), snow lines (Brauer et al. 2008a;Estrada et al. 2016;Dra ¸żkowska & Alibert 2017;Stammler et al. 2017;Hyodo et al. 2019), and pressure bumps in general (Pinilla et al. 2012), particles can accumulate and grow to larger sizes, reaching concentrations where the dust back-reaction may be strong enough to alter the dynamics of the gas (Taki et al. 2016;Onishi & Sekiya 2017;Kanagawa et al. 2017;Gonzalez et al. 2017;Dipierro et al. 2018).
In particular, the water snow line acts as a traffic jam for the dust if there is a change in the fragmentation velocity between silicates and ices (Birnstiel et al. 2010;Dra ¸żkowska & Alibert 2017;Pinilla et al. 2017).Previous results showed that the icy particles outside the snow line can grow to larger sizes (Gundlach et al. 2011) and drift faster to the inner regions.After crossing the snow line, the ice on the solid particles evaporates, leaving only dry silicates behind.Then the silicates in the inner regions fragment to smaller sizes and drift at lower speeds, creating a traffic jam.The traffic jam effect can concentrate enough material to trigger the formation of planetesimals through streaming instability (Schoonenberg & Ormel 2017;Dra ¸żkowska & Alibert 2017;Dra ¸żkowska & Dullemond 2018).
In this paper we study the dynamical effect of the snow line on gas dynamics by considering the effect of the dust backreaction onto the gas.We want to find the conditions under which the dust can slow down or revert the gas accretion rate and to test if further structures can appear beyond the snow line.
We used one-dimensional simulations that consider gas and dust advection, dust growth, and the back-reaction effects.To treat the global evolution of the disk we used the model of Birnstiel et al. (2012), which includes the size evolution of solids by using representative species.We implemented the modifications introduced by Dra ¸żkowska & Alibert (2017) which model the evaporation and recondensation of water at the snow line.
The paper is structured as follows.In Sect.2, we describe the gas and dust velocities considering the back-reaction and we present our model for the snow line.In Sect.3, we present the setup of our simulations and list the parameter space that we explored.In Sect.4, we show the conditions in which the accumulation of dust at the snow line results in strong backreaction effects that are capable of damping the accretion of gas to the inner regions.In Sect.5, we discuss the general effects of the back-reaction, when it should be considered, and what observational signatures might reveal dust-gas interactions in the inner regions.We summarize our results in Sect.6.We

Gas and dust evolution
The evolution of gas and dust can be described with the advection-diffusion equations as in Birnstiel et al. (2010): where r is the radial distance to the star, Σ is the surface density, v r is the radial velocity, and D d is the dust diffusivity.The subindex "g" and "d" denote the gas and dust, respectively.
An expression for the velocities can be obtained from the momentum conservation equations for both components (Nakagawa et al. 1986;Tanaka et al. 2005;Kanagawa et al. 2017;Dipierro et al. 2018), in which the gas experiences the stellar gravity, the pressure force, the viscous force, and the drag from multiple dust species, while each dust species only experiences the stellar gravity and the drag force from the gas.

Dust dynamics
Solving the momentum conservation equations, the radial and azimuthal velocities of the dust are given as per Weidenschilling (1977), Nakagawa et al. (1986), and Takeuchi & Lin (2002): where, for convenience, the dust azimuthal velocity is written relative to the Keplerian velocity, v K , as ∆v d,θ = v d,θ − v K .The same convention is used for the gas azimuthal velocity ∆v g,θ .
The Stokes number St is the dimensionless stopping time that measures the level of coupling of a dust species to the gas motion and is defined as: where Ω K is the Keplerian angular velocity and t stop is: with a the particle size, ρ s the material density of the solids, ρ g the gas density.The isothermal sound speed c s is: where k B is the Boltzmann constant, T the gas temperature, m H the hydrogen mass, and µ the mean molecular weight.From Eqs. ( 3) and ( 4), it can be inferred that small particles (St 1) move along with the gas, while large particles (St 1) are decoupled from it.Particles with St ∼ 1 experience the headwind from the gas with the strongest intensity and drift most efficiently towards the pressure maximum; in turn, these particles will also exert the strongest back-reaction onto the gas.
At the midplane, the Stokes number can be conveniently written as: The size of the particles is not static in time (Birnstiel et al. 2010(Birnstiel et al. , 2012) ) as the dust grows until it reaches the fragmentation barrier, where the particles are destroyed by high velocity collisions among themselves (Brauer et al. 2008b), or until the drift limit, where they drift faster than they can grow.The fragmentation barrier dominates the inner regions of the protoplanetary disk, and the maximum Stokes number that dust grains can reach before fragmenting is: where v frag is the fragmentation velocity which depends on the dust composition, and α t is the turbulence parameter for the dust fragmentation (Birnstiel et al. 2009).
Following Birnstiel et al. (2012), the drift limit can be approximated by: with v K the Keplerian velocity, and P the isothermal gas pressure at the midplane: with the gas scale height h g = c s /Ω K .Additionally, we assume that dust diffuses with with = Σ d /Σ g the vertically integrated dust-to-gas ratio, and ν the turbulent viscosity of the gas (Shakura & Sunyaev 1973): controlled by the viscous turbulence parameter α ν .We notice that, as in Carrera et al. (2017), our model considers two different turbulence parameters: α t for the dust turbulence (that controls the dust fragmentation, Eq. ( 9)) and α ν for the viscous turbulence (that controls the gas viscosity, Eq. ( 13)).
The (1 + ) −1 factor in Eq. ( 12) comes from considering that the dust concentration diffuses with respect to the gas and dust mixture, instead of the gas only.We neglect the (1 + St 2 ) −1 factor from Youdin & Lithwick (2007) since the particle sizes in our simulations remain small (St 2 1).

Gas dynamics
The gas velocities, considering the dust back-reaction onto the gas, have the following form: The gas velocity depends on the viscous velocity v ν , the pressure velocity v P , and the back-reaction coefficients A and B (Gárate et al. 2019).
The information related to the dust back-reaction is contained in the coefficients A and B, which in a dust free disk have values of A = 1 and B = 0.
In the absence of dust, the gas moves with the viscous velocity (Lynden-Bell & Pringle 1974): Similarly, if there is no dust, the gas orbits at sub-Keplerian speeds due to the pressure support, this pressure velocity is given by: ( The back-reaction coefficients are defined as follows: where X and Y are the following sums defined by (Tanaka et al. 2005;Okuzumi et al. 2012;Dipierro et al. 2018): where (m) is dust-to-gas ratio of the dust species with mass m, and the Stokes number can be related to the particle mass through m = 4πa 3 ρ s /3 and Eq. ( 8).
In principle, Eqs. ( 14) and (15) describe gas motion assuming that the dust is well mixed with the gas in the vertical direction (dust-to-gas ratio constant with the distance to the midplane).However, since dust grains settle towards the midplane (Dubrulle et al. 1995), the gas velocity in the midplane layers will be more affected by the dust back-reaction than the surface layers (Kanagawa et al. 2017;Dipierro et al. 2018).In Sect.2.2.1, we show how to calculate the corrected gas velocity derived from the net mass flux.
We discuss a physical interpretation for the back-reaction coefficients A, B in Sect.2.2.2 and provide an approximated expression valid for the case with a single dust species.

Effect of the vertical structure on the net mass flux
The corrected gas radial velocity vg,r can be obtained from the net mass flux, which is defined as: where z is the distance to the midplane.The vertical profile of the radial velocity v g,r (z) depends on: the vertical density distributions of gas and dust (ρ g (z), ρ d (m, z)), and the vertical profiles of the viscous and pressure velocities (v ν (z), v P (z)).
Assuming that the gas and dust are in vertical hydrostatic equilibrium, their respective density profiles are: where the vertical scale height of the dust species with mass m is defined in Birnstiel et al. (2010) as: From these profile we can obtain the dust-to-gas ratio of every particle species at height z with: and plug it into Eqs.( 20) and ( 21) to obtain the back-reaction coefficients A(z) and B(z) at every height.At this point we can actually generalize our expression for the gas radial velocity (Eq.( 14)) to every height v g,r (z) = A(z)v ν (z) + 2B(z)v P (z).Now the final step would be to define the vertical profiles of v ν (z) and v P (z), however, we find that assuming v ν (z) = v ν (as given in Eq. ( 16)), and v P (z) = v P (as given in Eq. ( 17)) is a good approximation.The net flux is then calculated with Eq. ( 22).We test the validity of our approximation in Appendix B and further discuss its physical interpretation.
It is worth noting that under this assumption, the radial gas velocity takes the following form: where Ā and B are the back-reaction coefficients corrected for the vertical structure (i.e., derived from the mass flux in Eq. ( 22)).
We repeat the same process to obtain the corrected radial velocity for the dust, by taking the net mass flux for each species of mass m: The gas and dust velocities derived from the net mass flux (v g,r and vd,r ) are used to transport the gas and dust in the advectiondiffusion Eqs.(1) and (2).

Understanding the back-reaction coefficients
While the back-reaction coefficients may seem rather obscure to interpret at first glance, they can be better understood as: (1) a "damping" factor (coefficient A) that slows the radial viscous evolution and reduces the pressure support and (2) a "pushing" factor (coefficient B) that tries to move the gas against the radial pressure gradient and adds some degree of pressure support to the orbital motion.
A quick estimate of these coefficients can be obtained if we consider the case of a single (well-mixed) particle species (Kanagawa et al. 2017;Dipierro et al. 2018;Gárate et al. 2019): From here we can see that both coefficients have values between 0 and 1 and that if the particles are small (St 2 ), then A ≈ ( + 1) −1 and B ≈ St ( + 1) −2 .
A&A 635, A149 (2020) In the case where the gas velocity v g,r is dominated by the viscous term such that Av ν > 2Bv P , the global evolution of gas and dust can be approximated as a damped viscous evolution.In Appendix A we further develop this idea and present a semianalytical test comparing the evolution of a simulation with backreaction and dust growth to the standard viscous evolution of a disk with a modified α ν parameter.
An equivalent expression for the gas and dust velocities, including the contribution of the back-reaction coefficients (Eqs.( 18) and ( 19)), can be found in Kretke et al. (2009) under the assumption that the particle sizes follow a single power-law distribution.
Further analysis on the effects of back-reaction considering different particle size distributions can be found in Dipierro et al. (2018), where the velocities given in Eqs. ( 14) and ( 15) are equivalent to their Eqs.( 11) and ( 12), while the integrals X and Y are equivalent to their λ i (Eq.( 17)).

Evaporation and recondensation at the snow line
To include the snow line in our simulations, we follow the model given by Dra ¸żkowska & Alibert (2017), which evolves four different species: a mix of hydrogen and helium, water vapor, silicate dust, and water ice that freezes over the silicate grains.
The gas phase is the sum of both hydrogen-helium and water vapor.It is traced by the surface density Σ g and advected according to Eq. ( 1).The water vapor, with surface density Σ vap , is advected with the same velocity as the gas but it also diffuses according to the concentration gradient.The mean molecular weight of the gas phase is then: where µ H 2 = 2.3 and µ vap = 18 are respectively the mean molecular weights of the hydrogen-helium mixture and the water vapor, and Σ H 2 = Σ g − Σ vap is the surface density of the standard hydrogen-helium mixture.
The dust grains are assumed to be a mixture of silicates and ices traced by Σ d , evolved according to Eq. ( 2), and have a material density of: where ρ sil = 3 g cm −3 and ρ ice = 1 g cm −3 are the densities of the silicates and ices, respectively, and Σ sil = Σ d − Σ ice is the surface density of the silicates.The composition of the dust grains determines the fragmentation velocity, where icy grains are stickier and can grow to larger sizes than the silicate grains.As in Dra ¸żkowska & Alibert (2017), we assume that the particles have the fragmentation velocity of ices v frag = 10 m s −1 (Wada et al. 2011;Gundlach et al. 2011;Gundlach & Blum 2015) if there is more than 1% of ice in the mixture, and the fragmentation velocity of silicates v frag = 1 m s −1 (Blum & Wurm 2000;Poppe et al. 2000;Güttler et al. 2010) otherwise.
The limit between evaporation and recondensation of water is given by the equilibrium pressure: P eq = P eq,0 exp(−A/T ), (33) with P eq,0 = 1.14 × 10 13 g cm −1 s −2 and A = 6062 K (Lichtenegger & Komle 1991;Dra ¸żkowska & Alibert 2017).The evaporation and recondensation of water are set to maintain the pressure of the water vapor at the equilibrium pressure (Ciesla & Cuzzi 2006), with: When the water vapor pressure is below this threshold (P vap < P eq ) the ice evaporates into vapor as follows: and vice-versa, if the vapor pressure is higher then it recondenses into ice with: where the factor next to ±(P vap − P eq ) transforms the pressure difference at the midplane into surface density.
As shown by Birnstiel et al. (2010) and Dra ¸żkowska & Alibert (2017), at the snow line, a traffic jam of dust is created because of the difference in the fragmentation velocities of silicates and ices.Recondensation also contributes by enhancing the amount of solids when the vapor diffuses and freezes back beyond the snow line (Stammler et al. 2017).

Simulation setup
We use the code twopoppy (Birnstiel et al. 2012) to study the global evolution of a protoplantary disk for 0.4 Myr around a solar mass star, advecting the gas and the dust according to the back-reaction velocities described in Sects.2.1 and 2.2, with the snow line model of Dra ¸żkowska & Alibert (2017) summarized above in Sect.2.3.

Two-population dust model
In twopoppy, the dust is modeled as a single fluid composed of two populations, an initial small particle particle population, and a large particle population with the size limited by the growth barriers (Eqs.( 9) and ( 10)), with a factor correction: St max = min(0.37× St frag , 0.55 × St drift ).
The dust velocity and the back-reaction coefficients are then calculated considering the mass fraction of the two populations.Birnstiel et al. (2012) found that the mass fraction of the large population if f m = 0.97 for the drift limited case, and f m = 0.75 for the fragmentation limited case.

Disk initial conditions
The gas surface density and temperature profile are defined by the following power laws: with r 0 = 1 au, Σ 0 = 1000 g cm −2 , T 0 = 300 K, p = 1 and q = 1/2.The disk surface density initially extends until r = 300 au.The disk size is intentionally large to provide a continuous supply of material during the simulation, and to make the interpretation of the back-reaction effects easier.We discuss the effect of the disk size in the outcome of the dust accumulation at the snow line in Sect.4.4.
We start the simulations with an uniform dust-to-gas ratio ε 0 such that Σ d = ε 0 Σ g , assuming that the solid material is composed of a mixture of 50 % ice and 50 % silicate (Lodders 2003, Table 11).The water vapor is introduced in the simulation as the ice evaporates.
The dust phase has a turbulence parameter of α t = 10 −3 , and an initial size of a 0 = 1 µm.

Grid and boundary conditions
The region of interest in our simulation extends from 0.1 to 300 au, with n r = 482 logarithmically spaced radial cells.To avoid any possible effects of the boundary conditions in our region of interest, we added 20 additional grid cells in the inner region between 0.05 and 0.1 au, and 58 additional grid cells in the outer region between 300 and 600 au.In total, our simulation consists of 560 grid cells from 0.05 to 600 au.
The additional cells at the inner region were added to avoid measuring the accretion rate onto the star too close to the inner boundary.The additional cells in the outer region were added to give the gas enough space to spread outwards without being affected by the outer boundary conditions.
At the inner boundary, we assume a constant slope for the quantity Σ g,d • r.At the outer boundary we have an open boundary condition for the gas and set a constant dust-to-gas ratio (but because of the additional grid cells, the gas never expands all the way to the outer boundary).To calculate the gas and dust velocities and take into account dust settling (Eqs.( 22) and ( 28)), we construct a local vertical grid at every radius with n z = 300 points, logarithmically spaced between 10 −5 and 10 h g .

Parameter space
The two most important parameters that control the strength of the back-reaction are the global dust-to-gas ratio ε 0 , and the gas viscous turbulence α ν .
We focus our study on three simulations with "Low", "Mid", and "High" global dust-to-gas ratios, with the respective values for ε 0 summarized in Table 1.
For the sake of clarity, through the paper, we use a single value for the viscous turbulence, with α ν = 10 −3 .This turbulence is low enough for the back-reaction effects to start affecting the gas dynamics (i.e., the term 2 Bv P becomes comparable to Āv ν in the gas velocity, Eq. ( 14)).
For completeness, in Appendix C we further extend our parameter space1 to include different values for the viscous turbulence α ν , although for simplicity, we keep the dust turbulence constant, with α t = 10 −3 .

Dust accumulation and gas depletion at the snow line
The evolution of gas is initially only dominated by the viscous accretion, but as time passes and dust grows, the back-reaction effects start to become dynamically important to the gas.At the water snow line, the Stokes number changes by two orders of magnitude (Fig. 1).In the inner disk, the particles can only grow to small sizes given by the fragmentation limit of silicates, while in the outer regions the dust size is limited by the fragmentation of water ice or the drift limit.The simulations with the higher dust-to-gas ratio show an increment in the Stokes number at the snow line location that is caused by the higher concentration of water vapor, which increases the fragmentation limit (by increasing the mean molecular weight and decreasing the sound speed; see Eqs. ( 7), (9), and (31)).
In the "Low ε 0 " simulation (Fig. 2, top panel), the change in particle size alone causes a traffic jam at the snow line location as the small dry silicates drift slower than the large icy particles, which results in a higher concentration of dust in the inner regions.Outside the snow line, the dust-to-gas ratio remains low, so the back-reaction from the large particles is not strong enough to perturb the gas.In this scenario, the gas surface density remains very close to the initial steady state.
Further effects can be seen in the "Mid ε 0 " simulation (Fig. 2, middle panel).First we notice an increment in the gas density profile at the snow line location, caused by the additional water vapor delivered by the icy grains (Ciesla & Cuzzi 2006).The water vapor and the dust are also more concentrated towards the snow line in this case as the higher dust-to-gas ratio damps the viscous velocity (| Āv ν | < |v ν |) more efficiently, slowing the diffusion of both gas and small particles.At the same time, the additional water vapor also increases the gas pressure which, in turn also increases the drift velocity of the large icy particles towards the snow line, resulting in higher dust concentrations.We also observe a small decrease in the gas surface density outside the snow line, caused by the dust back-reaction that slows down the gas velocity, reducing the supply to the inner regions.This effect becomes more pronounced for higher dust-to-gas ratios.
The back-reaction of dust onto the gas causes notorious perturbations in the "High ε 0 " simulation (Fig. 2, bottom panel).As in the "Mid ε 0 " simulation, the solids also accumulate at the snow line location, but now the icy dust particles outside the snow line exert a stronger push onto the gas, and reverse the gas (g/cm2) 0 : 5% Fig. 2. Surface density radial profiles of gas (red) and dust (blue) around the snow line.The dashed lines mark the initial conditions, and solid lines mark the simulation after 0.4 Myr.The dotted line marks the snow line at 0.4 Myr.Top, middle, and bottom panels: correspond to the cases with "Low", "Mid", and "High" ε 0 , respectively.accretion of the outer regions.This results in a depletion of gas outside the snow line (between r > 2.5 au), reaching a minimum density of ∼50% of its initial value.
Furthermore, the drop in gas density outside the snow line reduces the pressure gradient.Consequently, the drift speed of the large icy particles is also slowed down, allowing for an extended accumulation of dust in the outer regions.This process of gas depletion and dust accumulation is expected to continue as long as dust is supplied from the outer regions.
In the inner regions inside 1 au, the gas is depleted to ∼65% of its initial value.Only the additional water vapor supplied by the dust crossing the snow line prevents a further depletion of gas.The evolution of this simulation is illustrated in Fig. 3, where we can see the initial traffic jam caused by the change (g/cm2) 0 : 5% t: 0.4 Myr Fig. 3. Surface densities of gas (red), dust (blue), vapor (green) and ice (purple) of the "High ε 0 " simulation (ε 0 = 0.05), at different times.As time passes, dust accumulates around the snow line and the gas surface density is perturbed by the back-reaction. in particle size (t = 0.01 Myr), followed by a further concentration of solids once the vapor accumulates in the snow line (t = 0.1 Myr) and, finally, the depletion of gas outside the snow line, accompanied by the extended accumulation of dust (t = 0.4 Myr).
From Fig. 4, we see that the dust-to-gas ratios can reach extremely high values depending on the simulation parameters.The "Low ε 0 " simulation reaches a concentration of ≈ 0.1 in the inner regions (where the particles are small) because of the traffic jam but no further accumulation occurs outside the snow line.
In the "Mid ε 0 " case, the dust-to-gas ratio reaches a high value of ≈ 0.85 at the snow line and ≈ 0.4 at 1 au.The dust is more concentrated towards the snow line in this case because the back-reaction slows down the viscous diffusion (Eq.( 14)), yet, 10 0 10 1 r (AU) 10 2 10 1 10 0 0 : 1% 0 : 3% 0 : 5% as time passes, the dust should spread more evenly towards the inner regions.
The most extreme case is the "High ε 0 " simulation, where the dust accumulates both inside and outside the snow line.The dust accumulates in the inner regions due to the traffic jam caused by the change in particle size and the pressure maximum caused by the water vapor, reaching concentrations between ≈ 0.5 and 1.0.Outside the snow line the dust back-reaction depletes the gas and reduces the pressure gradient, creating another concentration point between 2.5 and 4 au where the dust-to-gas ratio reaches values of ≈ 1.0−2.0.The recondensation of vapor also contributes by enhancing the concentration of solids outside the snow line (Dra ¸żkowska & Alibert 2017; Stammler et al. 2017).

Accretion damped by the back-reaction
The radial velocity of the gas now depends not only on the viscous evolution, but also on the pressure gradient and the dust distribution (Eqs.( 14)-( 21)).Therefore, for high dust-to-gas ratios and large particles sizes, the gas flow may be damped and even reversed.
Figure 5 shows the gas velocities of the different simulations.In the "Low ε 0 " simulation the dust-to-gas ratio is higher in the inner regions (where grain sizes are small) and lower at the outer Gas accretion rate over time, measured at 0.5 au.The accretion rate decreases over time, dropping to a 85% of the initial value for the "Low ε 0 " simulation and to a 30−45% for the higher ε 0 cases.regions (where particle sizes are large).This trade-off between concentration and size means that the dust back-reaction does not dominate the evolution of the gas and that the gas velocity is only damped with respect to the steady state viscous velocity by a factor of a few.
The gas velocity is roughly v g,r ≈ 0.85 v ν inside the snow line and v g,r ≈ 0.80 v ν outside the snow line, where the transition is caused by the change in both particle size and dust-to-gas ratio.
This damping in the viscous velocity also leads into a similar decrease in the gas accretion rate onto the star, from Ṁ = 8 × 10 −9 M yr −1 to 6.8 × 10 −9 M yr −1 (Fig. 6).Once the dust supply is depleted, the accretion rate should return to its steady state value.
In the "High ε 0 " simulation, where the dust concentrations are high inside and outside the snow line, we can see the full effects of the dust back-reaction.In the inner regions (r < 2.5 au), the particles are small (St ∼ 10 −4 ), so the gas velocity is dominated by the term Āv ν , which corresponds to the viscous velocity damped by a factor of Ā ≈ (1 + ) −1 .In the outer region (r > 2.5 au) where the particles are large (St 10 −2 ), the velocity is dominated by the pressure velocity term 2 Bv P , which moves the gas outward, against the pressure gradient (Eq.( 14)).This reversal of the gas velocity causes the observed depletion in the gas surface density.Figure 7 shows the damping and pushing terms of the gas velocity to illustrate how the gas motion is affected by the dust back-reaction.
Since the gas inner disk is disconnected from the outer disk at the snow line in terms of mass transport, the accretion rate into the star is considerably reduced.As solid particles accumulate around the snow line, and the inner regions become more and more depleted of gas, the accretion rate reaches a value as low as Ṁ = 2.5 × 10 −9 M yr −1 .The only reason why the gas is not further depleted in the inner regions is because of the water vapor delivered by the icy dust particles crossing the snow line (Ciesla & Cuzzi 2006).
Meanwhile, the mass outside the snow line is transported outwards at a rate of ∼10 −9 −10 −8 M /yr.No instabilities seem to appear in gas surface density in the outer regions as the mass transported to the outer disk is only a small fraction of the total disk mass.Once the dust supply is exhausted, the back-reaction push will stop being effective and the gas accretion rate should retake the standard viscous evolution.
The behavior of the "Mid ε 0 " simulation is consistently in between the "Low ε 0 " and "High ε 0 " cases, with the gas flux practically frozen (v g,r ≈ 0) in the outer regions (r > 4 au).7. Gas velocity profile of the "High ε 0 " simulation after 0.4 Myr (black), and the decomposition of the two velocity terms Āv ν (red) and 2 Bv P (blue) (see Eq. ( 14)).In the inner regions, the pushing term 2 Bv P is negligible, as the particles Stokes number is too small, and the total velocity is dominated by the damped viscous velocity Āv ν .In the outer regions the term 2 Bv P overcomes the viscous evolution, and pushes gas against the pressure gradient.
All simulations show that the back-reaction push is particularly strong in a narrow region outside the snow line (between r ≈ 2.5 and 4 au), where the concentration of icy particles increases because of the recondensation of water vapor.
In Sect.5.3 we comment on the effects of dust settling on the accretion rate at different heights.

Depletion of H 2 and He inside the snow line
From the gas velocities, we see that in the cases where the back-reaction is effective it can stop or reverse the accretion of gas outside the snow line, causing the inner regions to become relatively depleted of gas.In particular, the dust back-reaction reduces the supply of H 2 , He to the inner regions as outside the snow line, this is the dominant gas component.
At the same time, the icy grains cross the snow line and deliver water vapor to the inner regions.Therefore, the gas will present a lower H 2 , He mass fraction in the inner disk than in the outer disk.
The total amount of water delivered to the inner regions depends on the initial dust-to-gas ratio 0 , while the dust backreaction affects how it is distributed.
Figure 8 shows that even in the "Low ε 0 " case, the mass fraction of H 2 , He is reduced to a 90%.
For the "Mid ε 0 " and "High ε 0 " cases the dust back-reaction onto the gas reduces the supply of light gases to the inner regions, creating environments dominated by water vapor inside the snow line, with a H 2 , He mass fraction between 40%−65%.The depletion is more concentrated towards the snow line because the damping term of the gas velocity ( Ā v ν ) slows down the viscous diffusion of water vapor.After the dust supply is exhausted, the region inside the snow line will be gradually refilled with gas from the outer regions in the viscous timescale (t ν ≈ 0.5 Myr at 4 au), and the H 2 , He mixture will be replenished to become the dominant component once more.

What happens without the back-reaction?
So far we have studied the impact on the dust back-reaction into the gas and dust density profiles, and in the gas velocity.So, how different is the situation when the back-reaction effect is ignored?In Fig. 9, we turn off the back-reaction effects (v g,r = v ν , ∆v g,θ = −v P ) and ignore the collective effect of dust on its diffusivity (D d = ν).The simulation with ε 0 = 0.01 shows only minor differences, corresponding to a faster dust accretion.This is an indication that for low dust-to-gas ratios, the back-reaction onto the gas is not important.For the simulations with ε 0 ≥ 0.03, we observe that without the back-reaction effect, the dust only concentrates in the inner regions due to the traffic jam caused by the change in particle sizes at the snow line.Accordingly, the water vapor delivered by the icy particles also increases the total gas content.
Figure 10 shows how the dust-to-gas ratio profile is affected by the dust back-reaction.Only when the back-reaction is considered the solid particles can pile up outside the water snow line due to the perturbed pressure gradient and the slower dust motion.For the simulations with ε 0 ≥ 0.03, the dust backreaction increases the dust-to-gas ratio by over an order of magnitude outside the water snow line.This is in agreement with the previous results of Dra ¸żkowska & Alibert (2017) and Hyodo et al. (2019), where the dust back-reaction was incorporated as the collective drift of the dust species.

The importance of the disk profile and size
How much the dust can perturb the gas surface density depends on the dust-to-gas ratio and the dust sizes, and also on how long the back-reaction is, effectively, acting.
In the "High ε 0 " case, the dust first creates a small depletion into the gas outside the snow line, the pressure slope changes and allows for large particles to further accumulate.Yet this scenario assumes that icy particles are being constantly delivered towards the snow line, while in reality the supply has a limit given by the disk size.
We made a test simulation with ε 0 = 0.05 as in the "High ε 0 " case, but this time starting with a self-similar profile (Lynden-Bell & Pringle 1974), following: with a cut-off radius of r c = 100 au.
From Fig. 11 we can see the evolution of this simulation until 1 Myr.Although we still observe that dust accumulates at the (g/cm2) 0 : 5% Fig. 9. Comparison of the surface density profiles when the backreaction is considered (solid lines) and ignored (dashed lines), after 0.4 Myr.For the cases with ε 0 ≥ 0.03, the gas surface density is reduced when the back-reaction is considered in the inner regions and the dust concentration is extended.snow line, reaching dust-to-gas ratios between = 0.7−0.8, and that the back-reaction push still creates a small dip in the gas surface density outside the snow line, the supply of solids is not enough to perturb the gas over extended periods of time.In this disk of limited size, no extended dust accumulation outside the snow line is observed.
The effect that still remains present is the decrease of the accretion rate (Fig. 12).As long as dust is delivered at the snow line, the accretion rate of gas is damped, and the mass fraction of the H 2 , He mixture is decreased in the inner regions.We find that between 0.4 and 0.5 Myr the dust concentration reaches its maximum value at the snow line (roughly the time required for the dust in the outer regions to grow and drift through the disk), and 10 0 10 1 r (AU) 10 2 10 1 10 0 0 : 1% 0 : 3% 0 : 5% BR -Off BR -On Fig. 10.Comparison of the dust-to-gas ratio profiles when the backreaction is considered (solid lines) and ignored (dashed lines) after 0.4 Myr.When the back-reaction is ignored, the dust accumulates only inside the snow line.
the accretion rate reaches its minimum of 3.0 × 10 −9 M yr −1 , where only 60% of the accretion flow corresponds to H 2 , He.
After 1 Myr the dust is completely depleted, the disk surface density roughly recovers the self similar profile and the accretion rate rises back again.

When is dust back-reaction important?
So far, we have seen that when the back-reaction is effective, it can enhance the dust concentration at the snow line (Fig. 4), damp the gas accretion rate (Fig. 6), and deplete the inner regions from hydrogen and helium (Fig. 8).
All of these effects can be traced back to the push exerted by the dust back-reaction onto the gas (Eq.( 14)) that reduces the pressure gradient (which enhances dust accumulation) and slows down the flux of material from outside the snow line to the inner regions.
As a rule of thumb, the gas dynamic is altered whenever the pressure velocity term is comparable to the damped viscous velocity ( Āv ν ∼ 2 Bv P , Eq. ( 14)), which occurs roughly when the particles have large Stokes number and high dust-to-gas ratios such that St /( + 1) ∼ α (Kanagawa et al. 2017;Dipierro et al. 2018).
In an inviscid disk (α ν ≈ 0), the gas velocity is dominated by the term 2 Bv P and the gas moves against the pressure gradient (Tanaka et al. 2005).On the other side, if the disk is highly turbulent (α ν St ), then the gas evolves with a damped viscous velocity Āv ν .In Appendix D, we include an equivalent criterion to determine the effect of the back-reaction, based on the angular momentum exchange between the dust and gas.
Through this paper we found that a high global dust-to-gas ratio of ε 0 0.03 and a low viscous turbulence of α ν 10 −3 (see Appendix C) are necessary for the back-reaction push to perturb the combined evolution of gas and dust.
We also showed that the duration and magnitude of these effects depends on the disk size, as the dust accumulation and the perturbation onto the gas stop once the solid reservoir is exhausted (Fig. 11).In particular, for a disk with cut-off radius of r c = 100 au the dust drifts from the outer regions to the snow line in 0.4 Myr.Afterwards, the back-reaction effects decay in a viscous timescale of the inner regions (roughly another 0.5 Myr).   . 11. Surface density profiles of gas (red) and dust (blue) at different times (solid lines).The initial condition corresponds to the self-similar profile (dashed lines).Top: simulation initially behaves in the same way as the power law profile until 0.1 Myr.Mid: at 0.4 Myr the dust supply gets exhausted before the back-reaction push can further deplete the gaseous disk.Bottom: after 1 Myr, the gas profile looks very similar to its initial condition, but most of the dust has been accreted.Moreover, part of the dust accumulated at the snow line will be converted into planetesimals through streaming instability (Youdin & Goodman 2005; Dra ¸żkowska & Alibert 2017) which, in turn, will reduce the dust-to-gas ratio and smear out the back-reaction effects.
We should keep in mind, however, that the results presented in this paper only occur if the snow line acts as a traffic jam for dust accretion, which is caused by the difference in the fragmentation velocities of dry silicates and icy aggregates.Yet recent studies suggest that there is no difference between the sticking properties of silicates and ices (Gundlach et al. 2018;Musiolik & Wurm 2019;Steinpilz et al. 2019), implying that the traffic jam should not actually form in the first place.M g M H2 M LBP Fig. 12. Accretion rate over time for the simulation with self-similar profile and ε 0 = 0.05.The gas accretion rate (red) decreases as the dust back-reaction damps gas velocity, and rises again after the dust is depleted.The accretion rate of H 2 , He (black) is even lower, as the gas supply of the outer regions is reduced at the snow line.The accretion rate of the standard self-similar solution (dotted line) is plotted for comparison.

Other scenarios where the back-reaction might be important
Similar traffic jams and dust traps can occur in different regions of the protoplanetary disk.Given high dust concentrations and large particles sizes, the dust back-reaction may perturb the gas in locations such as dead-zones (Kretke et al. 2009;Ueda et al. 2019;Gárate et al. 2019), the outer edge of gaps carved by planets (Paardekooper & Mellema 2004;Rice et al. 2006;Weber et al. 2018), and the edge of a photo-evaporative gap (Alexander & Armitage 2007).
In numerical models of protoplanetary disks, the backreaction effects should be considered when estimating the gas accretion rate (which is reduced by the interaction with the dust, Kanagawa et al. 2017), the planetesimal formation rate (which would be enhanced for higher dust concentrations, Dra ¸żkowska & Alibert 2017), or the width of a dusty ring in the outer edge of a gap carved by a planet (Kanagawa et al. 2018;Weber et al. 2018;Dra ¸żkowska et al. 2019).
The effects of the back-reaction could actually become effective at later stages of the disk lifetime, provided that other mechanisms continue to trap the dust delivered from the outer regions, for example, if a planet forms from the planetesimal population at the water snow line (Dra ¸żkowska & Alibert 2017), it would carve a gap that can effectively trap dust particles (Pinilla et al. 2012;Lambrechts et al. 2014), and create a new environment where the back-reaction can affect the gas and dust dynamics (Kanagawa et al. 2018).
On smaller scales, the dust back-reaction triggers the streaming instability, locally enhancing the concentration of dust particles until the solids become gravitationally unstable (Youdin & Goodman 2005), and close to the midplane, the friction between layers of gas and dust results in a Kelvin-Helmholtz instability between the two components (Johansen et al. 2006).
Finally, one scenario that we did not cover in our parameter space is when the turbulence is so low (α ν = 0) that the disk advection is reversed all the way to the inner boundary, which could lead to further perturbations at the snow line location, although a proper treatment of the dust sublimation should be included to account for this scenario.
Among our results, we could not reproduce the accumulation of dust in the outer regions of the disk described by Gonzalez et al. (2017) as the dust particles drift towards the inner regions before creating any perturbation in the outer gas disk.We also find that by taking into account the growth limits, the backreaction is less efficient than previously thought (Kanagawa et al. 2017) as the fragmentation barrier prevents the particles to grow to sizes beyond St frag and limiting the effect of the back-reaction even if the gas surface density decreases.
We do not expect our results to be significantly affected by changes in the disk mass or the stellar mass.Since particles sizes around the snow line are limited by the fragmentation barrier, the changes in any of these two parameters will only affect the physical size of the particles, but not their Stokes number (Eq. ( 9)), which controls the dynamical contribution of the particles to the gas motion.The timescales and the snow line location would change accordingly, but the qualitative results presented in this work should hold true.

Layered accretion by dust settling
Because large particles settle towards the midplane, the backreaction push onto the gas can be stronger at the disk midplane than at the surface (Kanagawa et al. 2017), which can result in the upper layers flowing inward (unperturbed by the dust), while the inner layers flow outward (due to the dust back-reaction).Depending on the particle sizes, this might result in different accretion rates at different heights.
While our approach to treat the vertical structure correctly traces the net mass transport (Sect.2.2.1, Appendix B) it does not provide information about layered accretion flow.To check if there is a substantial inflow of material at the upper layers, we calculate the accretion rate at every height (using the vertical model from Takeuchi & Lin 2002, see Appendix B) and measure the total mass inflow and outflow separately (Fig. 13).We find that inside the snow line (r < 2.7 au), where the dust particles are small (St ∼ 10 −4 ) and well mixed with the gas, the backreaction damps the gas motion uniformly at all heights, and the total inflow is of 3.0 × 10 −9 M yr −1 , only 6.0 × 10 −10 M yr −1 higher than the net accretion rate onto the star.
In the regions beyond the snow line dust accumulation (r > 3.0 au), we find that the accretion rate is layered, with the disk midplane flowing outward, while the surface layers move inward.The material inflow in this case is comparable to that of a dust free disk ( Ṁ ∼ 10 −8 M yr −1 ), even if the net mass flux is positive.This is in agreement with the results of Kanagawa et al. (2017).For these regions, we find that the dust back-reaction can revert the gas flow up to 2h d , which for the large dust particles at 3 au (St ∼ 10 −2 ) corresponds to 0.6h g .Interestingly, at the snow line location where the dust particles accumulate (2.7 au < r < 3.0 au), the dust back-reaction is strong enough to perturb the gas surface density.The steeper negative surface density slope found at the snow line causes the viscous accretion to be reduced or reversed at all heights (Takeuchi & Lin 2002).The accretion inflow is then reduced to a value of Ṁ = 8.0 × 10 −10 M yr −1 for the simulation with ε 0 = 0.03 (in which the reduced inflow only occurs above 0.7h g ), and to Ṁ = 6 × 10 −12 M yr −1 for the simulation with ε 0 = 0.05, where the inflow only occurs 2.5h g above the midplane.
The steepening of the surface density slope at the water snow line was not observed in the previous results of Kanagawa et al. (2017) as they did not include the snow line or a dust growth and recondensation model.We find that this perturbation caused by the dust accumulation at the snow line is key to reducing or stopping the accretion inflow over a wide vertical range, which can be larger than the dust scale height itself.
Given that the disk mass inside the snow line is of 2.0 × 10 −3 M , the composition of the gas phase described in Sect.4.2 should be corrected for the material flowing from the outer disk into the inner regions.For the simulation with ε 0 = 0.03, the H 2 + He ratio should be higher by a 20%, considering that the inflow is reduced by over an order of magnitude at the snow line location (though not completely stopped).For the simulation with ε 0 = 0.05, all our results hold.

Observational implications
The perturbation caused by the dust back-reaction at the snow line is only effective if the viscous turbulence is low, the dustto-gas ratio is high, and it only acts at early times of the disk evolution, when dust is supplied towards the inner regions.Given these constraints, we want to find out which disk properties would fit in this parameter space and what signatures we can expect to find if the back-reaction is effectively perturbing the gas.

Ideal targets
Young Class 0 and Class I disks seem to have typical sizes around 100-200 au (Najita & Bergin 2018, Table 1), so solids can be delivered to the inner regions only until 0.5−1 Myr, before A149, page 11 of 18 A&A 635, A149 (2020) the disk is depleted of dust (unless a pressure bump prevents particles from moving towards the star).This means that older disks (t > 1 Myr) are unlikely to present any perturbation from the back-reaction push.Then, among young disks and assuming viscous accretion, only those with low accretion rates of: could be subject to the back-reaction damping, as a low viscous evolution (α ν 10 −3 ) is required for the dust to affect the gas.
In terms of the dimensionless accretion parameter introduced by Rosotti et al. (2017), defined as: a disk of age τ would require η 0.1 for the dust back-reaction to effectively perturb the gas.

On the gas orbital velocity
If the concentration of dust in any region is high, then the gas pressure support is reduced and the orbital velocity approaches to the Keplerian velocity v K (Eq.( 15)).
At the midplane, where large grains concentrate, the gas motion deviates from the Keplerian velocity by: where the √ St/α t factor measures the concentration of large particles at the midplane by settling (see Appendix B).If in our disk the initial pressure velocity around the snow line was v P ≈ 2 × 10 −3 v K , then the dust back-reaction and the accumulation of water vapor makes the gas orbit at velocities of ∆v g,θ ≈ 7 × 10 −4 v K .At 2.7 au, where the snow line is located in our simulations, this correspond to a difference from the Keplerian velocity of approximately ∆v g,θ ≈ 10 m s −1 .
We expect that in future observations, the deviations from the Keplerian velocity could be used to constrain the dust content.Teague et al. (2018) already showed that the deviations from the Keplerian velocity can be used to kinematically detect a planet, reaching a precision of 2 m s −1 .Better characterizations of the orbital velocity profiles in dust rings may then be used to differentiate between a planet perturbation and a dust back-reaction perturbation, based on the profile shape.
Unfortunately, the spatial resolution required to observe this variation is less than 10 mas, for a disk at a distance of 100 au and a snow line at 3 au from the star, and next-generation instruments would be required.The velocity deviation could be easier to detect for disks around Herbig stars where the snow line is located at larger radii.

Shadows casted by dust accumulation
A recent study by Ueda et al. (2019) has shown that dust can accumulate at the inner edge of a dead zone (a region with low ionization and low turbulence, Gammie 1996) and cast shadows that extend up to 10 au.We notice that our accumulation of dust at the snow line is similar to the dead zone scenario, in the sense that high dust-to-gas ratios are reached in a narrow region of the inner disk (Fig. 4).Therefore, we hypothesize that similar shadows could be found in the regions just outside the snow line if enough dust is present.Still, radiative transfer simulations would be needed to determine the minimum dust-to-gas ratio necessary to cast a shadow.

Effects of the snow line traffic jam
The fast drift of the icy particles particles and the traffic jam at the snow line results in the accumulation of both small silicate dust and water vapor inside the snow line, even if the effect of the dust back-reaction is ignored (see Figs. 9 and 10).We find that if the initial dust supply is large enough (high 0 and large disk size) then, during the early stages of the disk evolution (t 1 Myr), we can expect the material accreted into the star to be rich in silicates and refractory materials carried by the dust (see Fig. 4), rich in oxygen (which is carried by the water vapor), and relatively poor in hydrogen, helium, and other volatile elements mixed with the gas outside the water snow line (see Fig. 8), such as nitrogen and neon.The X-ray emission could provide estimates of the abundance ratios in the accreted material (Günther et al. 2006), although the coronal emission of neon in young stars could mask some of these abundances (H.M. Günther, priv. comm.).
The increased concentration of water vapor in the warm inner regions would also enhance the emission from the water rotational lines.These lines have been already detected in different disks (Carr & Najita 2008;Salyk et al. 2008) in the mid-IR with Spitzer IRS, and could be further observed in the future using Mid-Infrared Instrument at the James Webb Space Telescope (MIRI, Rieke et al. 2015).Additionally, the excess of water should lead to low C/O ratios inside the snow line for young protoplanetary disks (Öberg et al. 2011;Booth & Clarke 2018).

Summary
In this study, we include the effects of the dust back-reaction on the gas in a model of the water snow line, which is known to act as a concentration point for dust particles due to the change in the fragmentation velocity between silicates and ices, together with the recondensation of water vapor into the surface of icy particles (Dra ¸żkowska & Alibert 2017).Our model shows how the dust back-reaction can perturb the gas dynamics and disk evolution, although the parameter space required for this to happen is limited.
In the vicinity of the snow line, provided that the global dustto-gas ratio is high (ε 0 0.03) and the viscosity low (α ν 10 −3 ), the effects of the dust back-reaction are to: -revert the net gas flux outside the snow line; -reduce the gas inflow at the snow line by over an order of magnitude; -damp the gas accretion rate onto the star to a 30−50% of its initial value; -reduce the hydrogen-helium content in the inner regions and concentrate water vapor at the snow line; -concentrate solids at the snow line, reaching dust-to-gas ratios of 0.8.These effects build up as long as dust is supplied from the outer disk into the snow line, with the duration set by the growth and drift timescale of the outer regions.After the dust reservoir is exhausted, the back-reaction effects decay in the viscous timescale of the inner regions.For a disk with size r c = 100 au, we find that dust accumulates only during the first 0.4 Myr and that the perturbation onto the gas has disappeared by the age of 1 Myr.

A149, page 12 of 18
The high dust-to-gas ratios required to trigger the backreaction effects and the traffic jam at the snow line can result in an enhanced water content in the inner regions, as well as in the accretion onto the star to be enriched with refractory materials and oxygen and, perhaps, a shadow that would be cast outside the snow line location by the accumulation of dust particles.Other types of dust traps could present similar behaviors, although each case must be revisited individually to evaluate the magnitude of the perturbation of the back-reaction into the gas velocity.3)) obtained from the simulations (red -DustPy, blue -TwoPopPy), and the analytical limits given by α eq, min (green) and α eq, max (black).The value obtained from the simulations is in between the two limits, in agreement with the analytical model.Right: accretion rate measured from the simulations, and the steady state accretion rate for the different α eq limits.
-the disk is initialized with a fully grown particle distribution (so that the back-reaction effects are uniform through the disk); -the back-reaction coefficients (in this test case) are implemented assuming that the dust-to-gas ratio is vertically uniform.If the simulations are working properly, then the disk will remain in steady state, and the accretion rate will be constant in radius with a value given by the damped equivalent viscosity α eq (Eq.(A.6)).Since in this test case all the particles are small < α ν ) and the size distribution is constant with radius, the dust-to-gas ratio and the back-reaction effects should also remain approximately uniform in time.
As shown in Fig. A.1, after 0.1 Myr the disk surface density between 5 and 100 au remains close to the steady state profile, with a deviation of less than 0.1% relative to its initial value.Figure A.2 shows that the mass accretion rate of the gas in the simulations is Ṁ ≈ 5.6 × 10 −8 M yr −1 and constant through the disk, in agreement with a steady state solution.More importantly, the value of the accretion is constrained between the minimum and maximum values given by α eq, min and α eq, max and Eq.(A.6).In terms of the viscous accretion, the back-reaction effect in our setup is equivalent to reduce the viscous turbulence α ν to a value of α eq ≈ 0.57 α ν .
Both twopoppy and DustPy deliver similar results, with a relative difference of roughly 5% in the α eq and Ṁ values.From here we can conclude that the back-reaction effects observed in the two population model are expected to be in agreement with those from a proper particle distribution.

A.2. Where the viscous approximation breaks
While we can always write the gas velocity in the form of Eq. (A.4) using the α eq parameter (Eq.(A.3)), the global disk evolution will still differ from a regular viscous evolution (unless α eq ∝ α ν ) as the value of γ ν does not depend on the slope of α eq .
In particular, the back-reaction effects cannot be treated as a viscous process if St /( + 1) α ν (Dipierro et al. 2018).In this case, the back-reaction push becomes more important than the inward viscous transport and results in negative equivalent α eq values, meaning that mass will be transported against the pressure gradient.Also, in the outer regions of the disk where the surface density profile becomes steeper (as in the self-similar solution Lynden-Bell & Pringle 1974), the viscous evolution spreads the gas outwards (γ ν < 0, v ν > 0).In these regions the dust backreaction pushes the gas in the same direction as the viscous spreading (2Bv P > 0) and, therefore, contributes to evolve the outer disk faster than the inner disk.

Appendix B: Modeling the vertical structure
In Sect.2.2.1, we discussed how to obtain the gas radial velocity from the net mass flux.For our simulations we considered the effect of the vertical structure of the gas and the settling of the dust on the back-reaction coefficients, but ignored the vertical profile of the pressure velocity v P (z) and viscous velocity v ν (z).In this appendix, we show that our results hold if we assume a standard vertical profile for the viscous and pressure velocity, and why our simple approximation works in the same way.
B.1.Vertical profiles for v ν and v P Following the Takeuchi & Lin (2002) model (see also Kanagawa et al. 2017;Dipierro et al. 2018), the vertical velocity profiles of v ν , v P are: where in this case, p and q are the local exponents of the gas surface density and temperature profiles.We include this information in the gas and dust velocity derived from the mass fluxes (Eqs.( 22) and ( 28)) and check how our results are affected.Figure B.1 shows the dust-to-gas ratio profile, and the evolution of the gas accretion present only minor differences when considering the vertical structure for the viscous and pressure velocities (Eqs.(B.1) and (B.2)).For the simulations with ε ≥ 0.03, only at the snow line location is the dust-to-gas ratio more spread over the radii, reducing the maximum concentration by a factor of a few.Consequently, the accretion rate onto the star is approximately a 5% higher when using the Takeuchi & Lin (2002) prescription.We find that our approximation reproduces the radial velocity profile calculated with the Takeuchi & Lin (2002) model (Fig. B.2) well, except in a narrow region beyond the snow line where the change in the gas surface density slope creates a spike in the gas velocity.However, this variation does not alter the rest of the simulation fields and our results are maintained independently of the vertical structure prescription used for the viscous and pressure velocities.

B.2. Explaining the vertical approximation
There are two distinguishable regimes for the gas and dust interaction.The first regime is the one in which the particles are small (St < α ν ).In this case, the particles are well mixed with the gas and, therefore, the dust-to-gas ratio (and the back-reaction coefficients) is uniform in the vertical direction.In this regime, the back-reaction pushing term is also negligible, and the gas velocity is well approximated by the damped viscous velocity vg,r ≈ Āv ν .
The second regime is the case in which the particles are large (St α ν ).In this case, the dust settles towards the midplane and the back-reaction push becomes important.The dust particles near the midplane will move the gas against the local pressure gradient with a velocity of v g,r (z h d ) ≈ 2 B(z h d )v P .Meanwhile, the upper layers of the disk have a low dust concentration and allow the gas to move inward with the viscous velocity, therefore the gas velocity above the characteristic dust scale height can be approximated by v g,r (z h d ) ≈ v ν .Then the velocity derived from the net flux (Eq.( 22)) will yield a good approximation considering the upper layers flowing inward (dominated by the viscous flow) and the midplane layers flowing against the pressure gradient (dominated by the back-reaction push).
This approximation seems to be valid for most of the disk, except on a narrow region where slope of the gas density is reversed (at r ≈ 3 au).However, this seems to be more related to a resolution problem than a physical reason; the surface density slope was smoothed in this case to avoid numerical problems in this region.Because of this region, we prefer using an approximate solution over the Takeuchi & Lin (2002) prescription.Solving the integral we obtain: From Eq. (D.3) we can also infer that the dust drifting from a radius r 1 to r 0 loses angular momentum (and delivers it to the gas) at a rate of: Jdrift = v 0 r 0 Ṁd (r 1 ) where the accretion rate of dust at r 1 is: Ṁd (r 1 ) = 2πr 1 d (r 1 )v d,r (r 1 ), (D.7) where we will assume an uniform dust-to-gas ratio, using Σ d = ε Σ g .Then, considering only the drifting component of the dust velocity (see Eq. ( 3)) we obtain that: where have taken the limit of small particles (St 1), and expanded the expression for the pressure velocity given in Eq. (A.2) using h/r = c s /v K .Now we can write the gap opening timescale as: with C(x) a function of the ratio x = r 1 /r 0 : We apply this formula to our simulations, using a scale height of h/r = 0.05, γ P = 2.75, and study the time required to clear a gap between r 0 = 2.5 au (which is approximately the location of the snow line) and r 1 = 10 au, and find: where we have taken r 0 = 2.5 au and r 1 = 10 au.We note that this condition is similar to 2Bv P Av ν (Eq.( 14)), which indicates whether gas motion is locally dominated by the dust back-reaction.In conclusion, looking at the values of t clear and t ν , it is easy to see why some disks create a gap-like perturbation in Fig. C.1.The disks with α = 10 −2 have viscous timescales that are too short in comparison with the clearing timescale by an order of magnitude, while the disks with α = 10 −4 are easily dominated by the dust back-reaction, provided that dust is delivered for enough time to complete a clearing timescale.

Fig. 1 .
Fig. 1.Stokes number radial profile after 0.4 Myr.Inside the water snow line (located between 2.5 and 3.0 au) the dust can grow only up to St ∼ 10 −4 .Outside the snow line it can reach values of St ∼ 10 −2 −10 −1 .
Fig.6.Gas accretion rate over time, measured at 0.5 au.The accretion rate decreases over time, dropping to a 85% of the initial value for the "Low ε 0 " simulation and to a 30−45% for the higher ε 0 cases.

Fig. 8 .
Fig. 8. H 2 , He mass fraction profile after 0.4 Myr.The mass fraction of light gases is lower inside the snow line as the dust crossing the snow line delivers water vapor.As the global dust-to-gas ratio increases, the back-reaction push outside the snow line reduces the flux of H 2 , He into the inner regions.

Fig
Fig. 11.Surface density profiles of gas (red) and dust (blue) at different times (solid lines).The initial condition corresponds to the self-similar profile (dashed lines).Top: simulation initially behaves in the same way as the power law profile until 0.1 Myr.Mid: at 0.4 Myr the dust supply gets exhausted before the back-reaction push can further deplete the gaseous disk.Bottom: after 1 Myr, the gas profile looks very similar to its initial condition, but most of the dust has been accreted.

Fig. 13 .
Fig. 13.Top: mass flux for the simulation with ε 0 = 5%, in the radial and vertical direction, obtained using the Takeuchi & Lin (2002) vertical velocity profiles.The blue regions show the material outflow, and the red regions show the inflow.Bottom: accretion inflow (red), outflow (blue), and the total mass accretion (dotted) profiles.
Fig. A.2. Left: equivalent α ν value (Eq.(A.3)) obtained from the simulations (red -DustPy, blue -TwoPopPy), and the analytical limits given by α eq, min (green) and α eq, max (black).The value obtained from the simulations is in between the two limits, in agreement with the analytical model.Right: accretion rate measured from the simulations, and the steady state accretion rate for the different α eq limits.
Fig. B.1.Comparison of the main results (top: dust-to-gas ratio profiles, bottom: accretion rate evolution) using our simple approximation for the vertical structure (solid lines) vs. the Takeuchi & Lin (2002) vertical structure model (dotted lines).
v d,r (r 1 ) = 2St v P (r 1 ) = −St γ P need to compare it with the viscous timescale, which can be understood in this case as the time necessary to close the gap: the snow line location of r 0 = 2.5 au, gives us a time of: we can derive from Eqs. (D.1), (D.9) and (D.12) a condition on α, ε , and St to see if the dust back-reaction can clear a gap: