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

TL;DR

This study uses numerical simulations to show that dust back-reaction at the snowline can significantly reduce gas accretion and cause dust accumulation, affecting planet formation in protoplanetary disks.

Contribution

It demonstrates the impact of dust back-reaction on gas dynamics at the snowline, highlighting conditions under which it can halt or reverse gas flow.

Findings

Dust back-reaction can reduce gas accretion to below 50%.

Dust accumulates at the snowline with high dust-to-gas ratios.

Effectiveness depends on initial dust ratio, turbulence, disk size, and age.

Abstract

Context. The water snowline divides dry and icy solid material in protoplanetary disks, and has been thought to significantly affect planet formation at all stages. If dry particles break up more easily than icy ones, then the snowline causes a traffic 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 snowline onto the gas dynamics. Methods. Using numerical simulations, we model the global radial evolution of an axisymmetric protoplanetary disk. Our model includes particle growth, evaporation and recondensation of water, and the back-reaction of dust onto the gas, taking into account the vertical distribution of dust particles. Results. We find that the dust back-reaction can stop and even reverse the net flux of gas outside the snowline, decreasing the gas accretion rate onto the star to under $50\%$ of its initial value. At the same time the dust accumulates at the snowline, reaching dust-to-gas ratios of $\epsilon \gtrsim 0.8$, and delivers large amounts of water vapor towards the inner disk, as the icy particles cross the snowline. However, the accumulation of dust at the snowline and the decrease in the gas accretion rate only take place if the global dust-to-gas ratio is high ($\varepsilon_0 \gtrsim 0.03$), if the viscous turbulence is low ($\alpha_\nu \lesssim 10^{-3} $), if the disk is large enough ($r_c \gtrsim 100\, \textrm{au}$), and only during the early phases of the disk evolution ($t \lesssim 1\, \textrm{Myr}$). Otherwise the dust back-reaction fails to perturb the gas motion.

Tables

Table 1. Table 1: Parameter space.

Simulation	${\epsilon }_0$
Low ${\epsilon }_0$	0.01
Mid ${\epsilon }_0$	0.03
High ${\epsilon }_0$	0.05

Equations

$$\frac{\partial }{\partial t}\left(r{\Sigma }_{\text{g}}\right)+\frac{\partial }{\partial r}(r{\Sigma }_{\text{g}}{v}_{\text{g,r}})=0,$$

$$\frac{\partial }{\partial t}\left(r{\Sigma }_{\text{d}}\right)+\frac{\partial }{\partial r}(r{\Sigma }_{\text{d}}{v}_{\text{d,r}})-\frac{\partial }{\partial r}\left(r{D}_{\text{d}}{\Sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{d}}}{{\Sigma }_{\text{g}}}\right)\right)=0,$$

$$v_{\text{d},r}=\frac{1}{1+{\mathrm{St}}^2}{v}_{\text{g},r}+\frac{2\mathrm{St}}{1+{\mathrm{St}}^2}\Delta v_{\text{g},\theta },$$

$$\Delta v_{\text{d},\theta }=\frac{1}{1+{\mathrm{St}}^2}\Delta v_{\text{g},\theta }-\frac{\mathrm{St}}{2(1+{\mathrm{St}}^2)}{v}_{\text{g},r},$$

$$\mathrm{St}=t_{\text{stop}}{\Omega }_K,$$

$$t_{\text{stop}}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_s}{{\rho }_{\text{g}}}\frac{a}{c_s},$$

$$c_s=\sqrt{\frac{k_{\text{B}}T}{\mu m_{\text{H}}}},$$

$$\mathrm{St}=\frac{\pi }{2}\frac{a{\rho }_s}{{\Sigma }_{\text{g}}}.$$

$${\mathrm{St}}_{\text{frag}}=\frac{1}{3}\frac{v_{\text{frag}}^2}{{\alpha }_t{c}_s^2},$$

$${\mathrm{St}}_{\text{drift}}={\left|\frac{\text{dln}P}{\text{dln}r}\right|}^{-1}\frac{v_K^2}{c_s^2}ϵ,$$

$$P=\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{h}_{\text{g}}}{c}_s^2,$$

$$D_{\text{d}}=\frac{\nu }{(1+ϵ)},$$

$$\nu ={\alpha }_{\nu }{c}_s^2{\Omega }_K^{-1},$$

$$v_{\text{g},r}=A{v}_{\nu }+2B{v}_P,$$

$$\Delta v_{\text{g},\theta }=-A{v}_P+\frac{1}{2}B{v}_{\nu }.$$

$$v_{\nu }=-\frac{3}{{\Sigma }_{\text{g}}\sqrt{r}}\frac{\partial }{\partial r}(\nu {\Sigma }_{\text{g}}\sqrt{r}).$$

$$v_P=-\frac{1}{2}{\left(\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{h}_{\text{g}}}{\Omega }_K\right)}^{-1}\frac{\partial P}{\partial r}.$$

$$A=\frac{X+1}{Y^2+(X+1{)}^2},$$

$$B=\frac{Y}{Y^2+(X+1{)}^2},$$

$$X=\sum _m\frac{1}{1+\mathrm{St}(m{)}^2}ϵ(m),$$

$$Y=\sum _m\frac{\mathrm{St}(m)}{1+\mathrm{St}(m{)}^2}ϵ(m),$$

$${\Sigma }_{\text{g}}{\bar{v}}_{\text{g},r}={\int }_{-\infty }^{+\infty }{\rho }_{\text{g}}(z)v_{\text{g},r}(z)\text{d}z,$$

$${\rho }_{\text{g}}(z)=\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{h}_{\text{g}}}\exp \left(-\frac{z^2}{2h_{\text{g}}^2}\right),$$

$${\rho }_{\text{d}}(z,m)=\frac{{\Sigma }_{\text{d}}(m)}{\sqrt{2\pi }{h}_{\text{d}}(m)}\exp \left(-\frac{z^2}{2h_{\text{d}}^2(m)}\right),$$

$$h_{\text{d}}(m)=h_{\text{g}}\cdot \min \left(1,\sqrt{\frac{{\alpha }_t}{\min (\mathrm{St},1/2)(1+{\mathrm{St}}^2)}}\right).$$

$$ϵ(z,m)=\frac{{\rho }_{\text{d}}(z,m)}{{\rho }_{\text{g}}(z)},$$

$${\bar{v}}_{\text{g},r}=\frac{1}{{\Sigma }_{\text{g}}}{\int }_{-\infty }^{+\infty }{\rho }_{\text{g}}(z)\left(A(z)v_{\nu }+2B(z)v_P\right)\text{d}z=\bar{A}{v}_{\nu }+2\bar{B}{v}_P,$$

$${\Sigma }_{\text{d}}(m){\bar{v}}_{\text{d},r}(m)={\int }_{-\infty }^{+\infty }{\rho }_{\text{d}}(z,m)v_{\text{d},r}(z,m)\text{d}z.$$

$$A_{\text{single}}=\frac{ϵ+1+{\mathrm{St}}^2}{(ϵ+1{)}^2+{\mathrm{St}}^2},$$

$$B_{\text{single}}=\frac{ϵ\mathrm{St}}{(ϵ+1{)}^2+{\mathrm{St}}^2}.$$

Full text

*Aims. *We aim to evaluate the effect of high dust concentrations around the snowline onto the gas dynamics.

*Methods. *Using numerical simulations, we model the global radial evolution of an axisymmetric protoplanetary disk. Our model includes particle growth, evaporation and recondensation of water, and the back-reaction of dust onto the gas, taking into account the vertical distribution of dust particles.

*Results. *We find that the dust back-reaction can stop and even reverse the net flux of gas outside the snowline, decreasing the gas accretion rate onto the star to under $50\%$ of its initial value. At the same time the dust accumulates at the snowline, reaching dust-to-gas ratios of $\epsilon\gtrsim 0.8$, and delivers large amounts of water vapor towards the inner disk, as the icy particles cross the snowline. However, the accumulation of dust at the snowline and the decrease in the gas accretion rate only take place if the global dust-to-gas ratio is high ($\varepsilon_{0}\gtrsim 0.03$), if the viscous turbulence is low ($\alpha_{\nu}\lesssim${10}^{-3}\text{,}), if the disk is large enough ($r_{c}\gtrsim$100\text{\,}\mathrm{a}\mathrm{u}), and only during the early phases of the disk evolution ($t\lesssim$1\text{,}\mathrm{M}\mathrm{y}\mathrm{r}$$). Otherwise the dust back-reaction fails to perturb the gas motion.

1 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 sub-keplerian speed due to its own pressure support.

On the other side, dust particles couple to the gas motion according to their size (Nakagawa et al., 1986; Takeuchi & Lin, 2002), 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, which causes them to drift towards the gas pressure maximum (Whipple, 1972; Weidenschilling, 1977), which in a typical disk is 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), outer edges of gaps carved by planets (Dipierro & Laibe, 2017; Kanagawa et al., 2018), snowlines (Brauer et al., 2008b; 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 snowline 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 snowline can grow to larger sizes (Gundlach et al., 2011) and drift faster to the inner regions. After crossing the snowline, 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 snowline on the gas dynamics, by considering the effect of the dust back-reaction onto the gas. We want to find under which conditions the dust can slow down or revert the gas accretion rate, and test if further structures can appear beyond the snowline.

We use 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 use the model of Birnstiel et al. (2012), that includes the size evolution of solids by using representative species, and implement the modifications introduced by Dra̧żkowska & Alibert (2017), that model the evaporation and recondensation of water at the snowline.

The paper is structured as follows. In section 2 we describe the gas and dust velocities considering the back-reaction, and present our model for the snowline. In section 3 we present the setup of our simulations and list the parameter space explored. In section 4 we show the conditions in which the accumulation of dust at the snowline results in strong back-reaction effects able to damp the accretion of gas to the inner regions. In section 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 section 6. We include a further study of the back-reaction equations and a semi-analytical test for the interested reader in Appendix A.

2 Gas and Dust evolution

The evolution of gas and dust can be described with the advection-diffusion equations as in Birnstiel et al. (2010):

[TABLE]

[TABLE]

where $r$ is the radial distance to the star, $\Sigma$ is the surface density, $v_{\textrm{r}}$ is the radial velocity, and $D_{\textrm{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 feels the stellar gravity, the pressure force, the viscous force, and the drag from multiple dust species, while each dust species only feels the stellar gravity and the drag force from the gas.

2.1 Dust Dynamics

Solving the momentum conservation equations, the radial and azimuthal velocities of the dust are as in Weidenschilling (1977); Nakagawa et al. (1986); Takeuchi & Lin (2002):

[TABLE]

[TABLE]

where for convenience the dust azimuthal velocity is written relative to the keplerian velocity $v_{K}$ as $\Delta v_{\textrm{d},\theta}=v_{\textrm{d},\theta}-v_{K}$. The same convention is used for the gas azimuthal velocity $\Delta v_{\textrm{g},\theta}$.

The Stokes number $\mathrm{St}$ is the dimensionless stopping time that measures the level of coupling of a dust species to the gas motion and is defined as:

[TABLE]

where $\Omega_{K}$ is the keplerian angular velocity and $t_{\textrm{stop}}$ is:

[TABLE]

with $a$ the particle size, $\rho_{s}$ the material density of the solids, $\rho_{\textrm{g}}$ the gas density. The isothermal sound speed $c_{s}$ is:

[TABLE]

where $k_{\textrm{B}}$ is the Boltzmann constant, $T$ the gas temperature, $m_{\textrm{H}}$ the hydrogen mass, and $\mu$ the mean molecular weight.

From Equation 3 and 4 it can be inferred that small particles ($\mathrm{St}\ll 1$) move along with the gas, while large particles ($\mathrm{St}\gg 1$) are decoupled from it. Particles with $\mathrm{St}\sim 1$ feel the head-wind 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:

[TABLE]

The size of the particles is not static in time (Birnstiel et al., 2010, 2012), dust grows until it reaches the fragmentation barrier, where the particles are destroyed by high velocity collisions among themselves (Brauer et al., 2008a), 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:

[TABLE]

where $v_{\textrm{frag}}$ is the fragmentation velocity which depends on the dust composition, and $\alpha_{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:

[TABLE]

with $v_{K}$ the Keplerian velocity, and $P$ the isothermal gas pressure at the midplane:

[TABLE]

with the gas scale height $h_{\textrm{g}}=c_{s}/\Omega_{K}$.

Additionally, we assume that dust diffuses with

[TABLE]

with $\epsilon=\Sigma_{\textrm{d}}/\Sigma_{\textrm{g}}$ the vertically integrated dust-to-gas ratio, and $\nu$ the turbulent viscosity of the gas (Shakura & Sunyaev, 1973):

[TABLE]

controlled by the viscous turbulence parameter $\alpha_{\nu}$.

Notice that as in Carrera et al. (2017), our model considers two different turbulence parameters: $\alpha_{t}$ for the dust turbulence (that controls the dust fragmentation, Equation 9), and $\alpha_{\nu}$ for the viscous turbulence (that controls the gas viscosity, Equation 13).

The $(1+\epsilon)^{-1}$ factor in Equation 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+\mathrm{St}^{2})^{-1}$ factor from Youdin & Lithwick (2007) since the particle sizes in our simulations remain small ($\mathrm{St}^{2}\ll 1$).

2.2 Gas Dynamics

The gas velocities, considering the dust back-reaction onto the gas, have the following form:

[TABLE]

[TABLE]

The gas velocity depends on the viscous velocity $v_{\nu}$, 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):

[TABLE]

Similarly, if there is no dust, the gas orbits at sub-keplerian speeds due to the pressure support, this pressure velocity is given by:

[TABLE]

The back-reaction coefficients are defined as follows:

[TABLE]

[TABLE]

where $X$ and $Y$ are the following sums defined by (Tanaka et al., 2005; Okuzumi et al., 2012; Dipierro et al., 2018):

[TABLE]

[TABLE]

where $\epsilon(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\pi a^{3}\rho_{s}/3$ and Equation 8.

In principle, Equation 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 section 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 section 2.2.2, and provide an approximated expression valid for the case with a single dust species.

2.2.1 Effect of the vertical structure on the net mass flux

The corrected gas radial velocity $\bar{v}_{\textrm{g,r}}$ can be obtained from the net mass flux, which is defined as:

[TABLE]

where $z$ is the distance to the midplane.

The vertical profile of the radial velocity $v_{\textrm{g},r}(z)$ depends on: the vertical density distributions of gas and dust ($\rho_{\textrm{g}}(z)$, $\rho_{\textrm{d}}(m,z)$), and the vertical profiles of the viscous and pressure velocities ($v_{\nu}(z)$, $v_{P}(z)$).

Assuming that the gas and dust are in vertical hydrostatic equilibrium, their respective density profiles are:

[TABLE]

[TABLE]

where the vertical scale height of the dust species with mass $m$ is defined in Birnstiel et al. (2010) as:

[TABLE]

From these profile we can obtain the dust-to-gas ratio of every particle species at height $z$ with:

[TABLE]

and plug it into Equation 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 (Equation 14) to every height $v_{\textrm{g},r}(z)=A(z)v_{\nu}(z)+2B(z)v_{P}(z)$.

Now, the final step would be to define the vertical profiles of $v_{\nu}(z)$ and $v_{P}(z)$, however we find that assuming $v_{\nu}(z)=v_{\nu}$ (as given in Equation 16), and $v_{P}(z)=v_{P}$ (as given in Equation 17) is a good approximation. The net flux is then calculated with Equation 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:

[TABLE]

where $\bar{A}$ and $\bar{B}$ are the back-reaction coefficients corrected for the vertical structure (i.e., derived from the mass flux in Equation 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$:

[TABLE]

The gas and dust velocities derived from the net mass flux ($\bar{v}_{\textrm{g},r}$ and $\bar{v}_{\textrm{d},r}$) are used to transport the gas and dust in the advection-diffusion equations 1 and 2.

2.2.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 a “damping” factor (coefficient $A$), that slows the radial viscous evolution and reduces the pressure support, and 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):

[TABLE]

[TABLE]

From here we can see that both coefficients have values between 0 and 1, and that if the particles are small ($\mathrm{St}^{2}\ll\epsilon$), then $A\approx(\epsilon+1)^{-1}$ and $B\approx\mathrm{St}\,\epsilon\,(\epsilon+1)^{-2}$.

In the case where the gas velocity $v_{\textrm{g,r}}$ is dominated by the viscous term such that $Av_{\nu}>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 semi-analytical test comparing the evolution of a simulation with back-reaction and dust growth to the standard viscous evolution of a disk with a modified $\alpha_{\nu}$ parameter.

An equivalent expression for the gas and dust velocities, including the contribution of the back-reaction coefficients (Equation 18, 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 Equation 14 and 15 are equivalent to their Eqs. 11 and 12, while the integrals $X$ and $Y$ are equivalent to their $\lambda_{i}$ (Eq. 17).

2.3 Evaporation and recondensation at the snowline

To include the snowline 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 $\Sigma_{\textrm{g}}$, and is advected according to Equation 1. The water vapor, with surface density $\Sigma_{\textrm{vap}}$, is advected with the same velocity as the gas, but also diffuses according to the concentration gradient. The mean molecular weight of the gas phase is then:

[TABLE]

where $\mu_{\textrm{H}_{2}}=2.3$ and $\mu_{\textrm{vap}}=18$ are respectively the mean molecular weights of the hydrogen-helium mixture and the water vapor, and $\Sigma_{\textrm{H}_{2}}=\Sigma_{\textrm{g}}-\Sigma_{\textrm{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 $\Sigma_{\textrm{d}}$, evolved according to Equation 2, and have a material density of:

[TABLE]

where $\rho_{\textrm{sil}}$ = $3\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$ and $\rho_{\textrm{ice}}$ = $1\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-3}$ are the densities of the silicates and ices, respectively, and $\Sigma_{\textrm{sil}}=\Sigma_{\textrm{d}}-\Sigma_{\textrm{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_{\textrm{frag}}=$10\text{,}\mathrm{m},\mathrm{s}^{-1} (Wada et al., [2011](#bib.bib60); Gundlach et al., [2011](#bib.bib23); Gundlach & Blum, [2015](#bib.bib22)) if there is more than $1\%$ of ice in the mixture, and the fragmentation velocity of silicates $v_{\textrm{frag}}=$1\text{\,}\mathrm{m}\,\mathrm{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:

[TABLE]

with $P_{\textrm{eq,0}}=$1.14\text{\times}{10}^{13}\text{,}\mathrm{g},\mathrm{c}\mathrm{m}^{-1}\mathrm{s}^{-2} and $A=$6062\text{\,}\mathrm{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:

[TABLE]

When the water vapor pressure is below this threshold ($P_{\textrm{vap}}<P_{\textrm{eq}}$) the ice evaporates into vapor as follows:

[TABLE]

and vice-versa, if the vapor pressure is higher then it recondenses into ice with:

[TABLE]

where the factor next to $\pm(P_{\textrm{vap}}-P_{\textrm{eq}})$ transforms the pressure difference at the midplane into surface density.

As shown by Birnstiel et al. (2010); Dra̧żkowska & Alibert (2017), at the snowline a traffic jam of dust is created because of the difference in the fragmentation velocities of silicates and ices. Recondensation also contributes to enhance the amount of solids when the vapor diffuses and freezes back beyond the snowline (Stammler et al., 2017).

3 Simulation Setup

We use the code twopoppy (Birnstiel et al., 2012) to study the global evolution of a protoplantary disk for $0.4\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ around a solar mass star, advecting the gas and the dust according to the back-reaction velocities described in section 2.1 and 2.2, with the snowline model of Dra̧żkowska & Alibert (2017) summarized above in section 2.3.

3.1 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 (Equation 9 and 10), with a factor correction: $\textrm{St}_{\textrm{max}}=\min(0.37\cdot\textrm{St}_{\textrm{frag}},0.55\cdot\textrm{St}_{\textrm{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_{\textrm{m}}=0.97$ for the drift limited case, and $f_{\textrm{m}}=0.75$ for the fragmentation limited case.

3.2 Disk Initial conditions

The gas surface density and temperature profile are defined by the following power laws:

[TABLE]

[TABLE]

with $r_{0}=$1\text{,}\mathrm{a}\mathrm{u}, $\Sigma_{0}=$1000\text{\,}\mathrm{g}\,\mathrm{c}\mathrm{m}^{-2}, $T_{0}=$300\text{,}\mathrm{K}$$, $p=1$ and $q=1/2$.

The disk surface density initially extends until $r=$300\text{,}\mathrm{a}\mathrm{u}$$. 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 snowline in section 4.4.

We start the simulations with an uniform dust-to-gas ratio $\varepsilon_{0}$ such that $\Sigma_{\textrm{d}}=\varepsilon_{0}\Sigma_{\textrm{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 $\alpha_{t}=${10}^{-3}\text{,}, and an initial size of $a_{0}=$1\text{\,}\mu\mathrm{m}.

3.3 Grid and Boundary Conditions

The region of interest in our simulation extends from $0.1\text{\,}$ to $300\text{\,}\mathrm{a}\mathrm{u}$, with $n_{r}=482$ logarithmically spaced radial cells.

To avoid possible effects of the boundary conditions in our region of interest, we add 20 additional grid cells in the inner region between $0.05\text{\,}$ and $0.1\text{\,}\mathrm{a}\mathrm{u}$, and 58 additional grid cells in the outer region between $300\text{\,}$ and $600\text{\,}\mathrm{a}\mathrm{u}$. In total, our simulation consist on $560$ grid cells from $0.05\text{\,}$ to $600\text{\,}\mathrm{a}\mathrm{u}$.

The additional cells at the inner region where 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 $\Sigma_{\textrm{g,d}}\cdot 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 (Equation 22 and 28) we construct a local vertical grid at every radius with $n_{z}=300$ points, logarithmically spaced between {10}^{-5}\text{\,}$\,h_{\textrm{g}}$ and 10\text{,}$\,h_{\textrm{g}}$.

3.4 Parameter Space

The two most important parameters that control the strength of the back-reaction are the global dust-to-gas ratio $\varepsilon_{0}$, and the gas viscous turbulence $\alpha_{\nu}$.

We will focus our study in three simulations with “Low”, “Mid”, and “High” global dust-to-gas ratios, with the respective values for $\varepsilon_{0}$ summarized in Table 1.

For the sake of clarity, through the paper we will use a single value for the viscous turbulence, with $\alpha_{\nu}=${10}^{-3}\text{,}$$. This turbulence is low enough for the back-reaction effects to start affecting the gas dynamics (i.e., the term $2\bar{B}v_{P}$ becomes comparable to $\bar{A}v_{\nu}$ in the gas velocity, Equation 14).

For completeness, in Appendix C we further extend our parameter space111The simulation data files, including the extended parameter space, are available in Zenodo: doi.org/10.5281/zenodo.3552597 to include different values for the viscous turbulence $\alpha_{\nu}$, though for simplicity we keep the dust turbulence constant, with $\alpha_{t}=${10}^{-3}\text{,}$$.

4 Dust accumulation and gas depletion at the snowline

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 snowline, the Stokes number changes by 2 orders of magnitude (Figure 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 snowline location, caused by the higher concentration of water vapor which increases the fragmentation limit (by increasing the the mean molecular weight, and decreasing the sound speed, see Equation 7, 9 and 31).

In the “Low $\varepsilon_{0}$” simulation (Figure 2, top panel), the change in particle size alone causes a traffic jam at the snowline 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 snowline 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 $\varepsilon_{0}$” simulation (Figure 2, middle panel). First we notice an increment in the gas density profile at the snowline 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 snowline in this case, as the higher dust-to-gas ratio damps more efficiently the viscous velocity ($|\bar{A}v_{\nu}|<|v_{\nu}|$), 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 snowline, resulting in higher dust concentrations.

We also observe a small decrease in the gas surface density outside the snowline, 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 $\varepsilon_{0}$” simulation (Figure 2, bottom panel). As in the “Mid $\varepsilon_{0}$” simulation, the solids also accumulate at the snowline location, but now the icy dust particles outside the snowline exert a stronger push onto the gas, and reverse the gas accretion of the outer regions. This results in a depletion of gas outside the snowline (between $r>$2.5\text{,}\mathrm{a}\mathrm{u}$$), reaching a minimum density of $\sim 50\%$ of its initial value.

Furthermore, the drop in gas density outside the snowline 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\text{\,}\mathrm{a}\mathrm{u}$, the gas is depleted to $\sim 65\%$ of its initial value. Only the additional water vapor supplied by the dust crossing the snowline prevents a further depletion of gas. The evolution of this simulation is illustrated in Figure 3, where we can see the initial traffic jam caused by the change in particle size ($t=$0.01\text{,}\mathrm{M}\mathrm{y}\mathrm{r}), followed by a further concentration of solids once the vapor accumulates in snowline ($t=$0.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}), and finally the depletion of gas outside the snowline, accompanied by the extended accumulation of dust ($t=$0.4\text{,}\mathrm{M}\mathrm{y}\mathrm{r}$$).

From Figure 4 we see that the dust-to-gas ratios can reach extremely high values depending on the simulation parameters. The “Low $\varepsilon_{0}$” simulation reaches a concentration of $\epsilon\approx 0.1$ in the inner regions (where the particles are small), because of the traffic jam, but no further accumulation occurs outside the snowline.

In the “Mid $\varepsilon_{0}$” case, the dust-to-gas ratio reaches a high value of $\epsilon\approx 0.85$ at the snowline, and $\epsilon\approx 0.4$ at $1\text{\,}\mathrm{a}\mathrm{u}$. The dust is more concentrated towards the snowline in this case because the back-reaction slows down the viscous diffusion (Equation 14), yet as time passes the dust should spread more evenly towards the inner regions.

The most extreme case is the “High $\varepsilon_{0}$” simulation, where the dust accumulates both inside and outside the snowline. 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 $\epsilon\approx 0.5-1.0$. Outside the snowline the dust back-reaction depletes the gas and reduces the pressure gradient, creating another concentration point between 2.5\text{\,}$-$4\text{\,}\mathrm{a}\mathrm{u} where the dust-to-gas ratio reaches values of $\epsilon\approx 1.0-2.0$. The recondensation of vapor also contributes to enhance the concentration of solids outside the snowline (Dra̧żkowska & Alibert, 2017; Stammler et al., 2017).

4.1 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 (Equation 14 to 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 $\varepsilon_{0}$” simulation the dust-to-gas ratio is higher in the inner regions (where grain sizes are small), and lower at the outer 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_{\textrm{g,r}}\approx 0.85\,v_{\nu}$ inside the snowline and $v_{\textrm{g,r}}\approx 0.80\,v_{\nu}$ outside the snowline, 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 $\dot{M}=$8\text{\times}{10}^{-9}\text{,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$$ to $6.8\text{\times}{10}^{-9}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$ (Figure 6). Once the dust supply is depleted, the accretion rate should return to its steady state value.

In the “High $\varepsilon_{0}$” simulation, where the dust concentrations are high inside and outside the snowline, we can see the full effects of dust back-reaction. In the inner regions ($r<$2.5\text{,}\mathrm{a}\mathrm{u}) the particles are small ($\mathrm{St}\sim${10}^{-4}\text{\,}), so the gas velocity is dominated by the term $\bar{A}v_{\nu}$, which corresponds to the viscous velocity damped by a factor of $\bar{A}\approx(1+\epsilon)^{-1}$. In the outer region ($r>$2.5\text{,}\mathrm{a}\mathrm{u}) where the particles are large ($\mathrm{St}\gtrsim${10}^{-2}\text{\,}), the velocity is dominated by the pressure velocity term $2\bar{B}v_{P}$, which moves the gas outward, against the pressure gradient (Equation 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 snowline in terms of mass transport, the accretion rate into the star is considerably reduced. As solid particles accumulate around the snowline, and the inner regions become more and more depleted of gas, the accretion rate reaches a value as low as $\dot{M}=$2.5\text{\times}{10}^{-9}\text{,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$$. 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 snowline (Ciesla & Cuzzi, 2006).

Meanwhile, the mass outside the snowline is transported outwards at a rate of $\sim${10}^{-9}\text{,}$-${10}^{-8}\text{,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$$. 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 $\varepsilon_{0}$” simulation is consistently in between the “Low $\varepsilon_{0}$” and “High $\varepsilon_{0}$” cases, with that the gas flux is practically frozen ($v_{\textrm{g,r}}\approx 0$) in the outer regions ($r>$4\text{,}\mathrm{a}\mathrm{u}). All simulations show that the back-reaction push is particularly strong in a narrow region outside the snowline (between $r\approx$2.5\text{\,}$-$4\text{\,}\mathrm{a}\mathrm{u}), where the concentration of icy particles increases because of the recondensation of water vapor.

In subsection 5.3 we comment on the effects of dust settling on the accretion rate at different heights.

4.2 Depletion of $\textrm{H}_{2}$ and He inside the snowline.

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 snowline, causing the inner regions to become relatively depleted of gas.

In particular, the dust back-reaction reduces the supply of the $\textrm{H}_{2},\textrm{He}$ to the inner regions, as outside the snowline this is the dominant gas component.

At the same time, the icy grains cross the snowline and deliver water vapor to the inner regions. Therefore, the gas will present a lower $\textrm{H}_{2},\textrm{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 $\epsilon_{0}$, while the dust back-reaction affects how it is distributed.

Figure 8 shows that even in the “Low $\varepsilon_{0}$” case, the mass fraction of $\textrm{H}_{2},\textrm{He}$ is reduced to a $90\%$.

For the “Mid $\varepsilon_{0}$” and “High $\varepsilon_{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 snowline, with a $\textrm{H}_{2},\textrm{He}$ mass fraction between $40\%-65\%$. The depletion is more concentrated towards the snowline because the damping term of the gas velocity ($\bar{A}\,v_{\nu}$) slows down the viscous diffusion of water vapor.

After the dust supply is exhausted, the region inside the snowline will be gradually refilled with gas from the outer regions in the viscous timescale ($t_{\nu}\approx$0.5\text{,}\mathrm{M}\mathrm{y}\mathrm{r}$$ at $4\text{\,}\mathrm{a}\mathrm{u}$), and the $\textrm{H}_{2},\textrm{He}$ mixture will be replenished to become the dominant component once more.

4.3 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 Figure 9 we turn off the back-reaction effects ($v_{\textrm{g,r}}=v_{\nu}$, $\Delta v_{\textrm{g},\theta}=-v_{P}$), and ignore the collective effect of dust on its diffusivity ($D_{\textrm{d}}=\nu$). The simulation with $\varepsilon_{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 $\varepsilon_{0}\geq 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 snowline. 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 snowline, due to the perturbed pressure gradient and the slower dust motion. For the simulations with $\varepsilon_{0}\geq 0.03$, the dust back-reaction increases the dust-to-gas ratio by over an order of magnitude outside the water snowline. This agrees with previous results of Dra̧żkowska & Alibert (2017); Hyodo et al. (2019) where the dust back-reaction was incorporated as the collective drift of the dust species.

4.4 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, but also on how long the back-reaction is effectively acting.

In the “High $\varepsilon_{0}$” case, the dust first creates a small depletion into the gas outside the snowline, 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 snowline, while in reality the supply has a limit given by the disk size.

We made a test simulation with $\varepsilon_{0}=0.05$ as in the “High $\varepsilon_{0}$” case, but this time starting with a self-similar profile (Lynden-Bell & Pringle, 1974), following:

[TABLE]

with a cut-off radius of $r_{c}=$100\text{,}\mathrm{a}\mathrm{u}$$.

From Figure 11 we can see the evolution of this simulation until $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. Though we still observe that dust accumulates at the snowline, reaching dust-to-gas ratios between $\epsilon=0.7-0.8$, and that the back-reaction push still creates a small dip in the gas surface density outside the snowline, 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 snowline is observed.

The effect that still remains present is the decrease of the accretion rate (Figure 12). As long as dust is delivered at the snowline, the accretion rate of gas is damped, and the mass fraction of the $\textrm{H}_{2},\textrm{He}$ mixture is decreased in the inner regions.

We find that between 0.4\text{\,}$-$0.5\text{\,}\mathrm{M}\mathrm{y}\mathrm{r} the dust concentration reaches its maximum value at the snowline (roughly the time required for the dust in the outer regions to grow and drift through the disk), and the accretion rate reaches its minimum of $3.0\text{\times}{10}^{-9}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$, where only $60\%$ of the accretion flow corresponds to $\textrm{H}_{2},\textrm{He}$.

After $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ the dust is completely depleted, the disk surface density roughly recovers the self similar profile and the accretion rate rises back again.

5 Discussion

5.1 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 snowline (Figure 4), damp the gas accretion rate (Figure 6), and deplete the inner regions from hydrogen and helium (Figure 8).

All of these effects can be traced back to the push exerted by the dust back-reaction onto the gas (Equation 14), that reduces the pressure gradient (which enhances dust accumulation), and slows down the flux of material from outside the snowline 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 ($\bar{A}v_{\nu}\sim 2\bar{B}v_{P}$, Equation 14), which occurs roughly when the particles have large Stokes number and high dust-to-gas ratios such that $\mathrm{St}\,\epsilon/(\epsilon+1)\sim\alpha$ (Kanagawa et al., 2017; Dipierro et al., 2018).

In an inviscid disk ($\alpha_{\nu}\approx 0$), the gas velocity is dominated by the term $2\bar{B}v_{P}$, and the gas moves against the pressure gradient (Tanaka et al., 2005). On the other side, if the disk is highly turbulent ($\alpha_{\nu}\gg\mathrm{St\epsilon}$), then the gas evolves with a damped viscous velocity $\bar{A}v_{\nu}$. 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 $\varepsilon_{0}\gtrsim 0.03$, and a low viscous turbulence of $\alpha_{\nu}\lesssim${10}^{-3}\text{,}$$ (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 (Figure 11). In particular, for a disk with cut-off radius of $r_{c}=$100\text{,}\mathrm{a}\mathrm{u}$$ the dust drifts from the outer regions to the snowline in $0.4\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$. Afterwards, the back-reaction effects decay in a viscous timescale of the inner regions (roughly another $0.5\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$).

Moreover, part of the dust accumulated at the snowline 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 snowline 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 form in the first place.

5.2 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 back-reaction 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 snowline (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 ($\alpha_{\nu}=0$) that the disk advection is reversed all the way to the inner boundary, which could lead to further perturbations at the snowline location, though 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 back-reaction is less efficient than previously thought (Kanagawa et al., 2017), as the fragmentation barrier prevents the particles to grow to sizes beyond $\textrm{St}_{\textrm{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 snowline 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 (Equation 9) which controls the dynamical contribution of the particles to the gas motion. The timescales and the snowline location would change accordingly, but the qualitative results presented in this work should hold true.

5.3 Layered accretion by dust settling

Because large particles settle towards the midplane, the back-reaction 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 traces correctly the net mass transport (Section 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 (Figure 13).

We find that inside the snowline ($r<$2.7\text{,}\mathrm{a}\mathrm{u}), where the dust particles are small ($\mathrm{St}\sim${10}^{-4}\text{\,}) and well mixed with the gas, the back-reaction damps the gas motion uniformly at all heights, and the total inflow is of $3.0\text{\times}{10}^{-9}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$, only $6.0\text{\times}{10}^{-10}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$ higher than the net accretion rate onto the star.

In the regions beyond the snowline dust accumulation ($r>$3.0\text{,}\mathrm{a}\mathrm{u}), 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 ($\dot{M}\sim${10}^{-8}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}), 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_{\textrm{d}}$, which for the large dust particles at $3\text{\,}\mathrm{a}\mathrm{u}$ ($\mathrm{St}\sim${10}^{-2}\text{,}$$) corresponds to $0.6h_{\textrm{g}}$.

Interestingly, at the snowline location where the dust particles accumulate (2.7\text{\,}\mathrm{a}\mathrm{u}$<r<$3.0\text{\,}\mathrm{a}\mathrm{u}), the dust back-reaction is strong enough to perturb the gas surface density. The steeper negative surface density slope found at the snowline 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 $\dot{M}=$8.0\text{\times}{10}^{-10}\text{,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r} for the simulation with $\varepsilon_{0}=0.03$ (in which the reduced inflow only occurs above $0.7h_{\textrm{g}}$), and to $\dot{M}=$6\text{\times}{10}^{-12}\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r} for the simulation with $\varepsilon_{0}=0.05$, where the inflow only occurs $2.5h_{\textrm{g}}$ above the midplane.

The steepening of the surface density slope at the water snowline was not observed in the previous results of Kanagawa et al. (2017) as they did not include the snowline or a dust growth and recondensation model. We find that this perturbation caused by the dust accumulation at snowline is key to reduce or stop 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 snowline is of $2.0\text{\times}{10}^{-3}\text{\,}\mathrm{M}_{\odot}$, the composition of the gas phase described in Section 4.2 should be corrected for the material flowing from the outer disk into the inner regions. For the simulation with $\varepsilon_{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 snowline location (though not completely stopped). For the simulation with $\varepsilon_{0}=0.05$ all our results hold.

5.4 Observational Implications

The perturbation caused by the dust back-reaction at the snowline is only effective if the viscous turbulence is low, if the dust-to-gas ratio is high, and only acts at early times of the disk evolution, while dust is supplied towards the inner regions. Given these constraints, we want to find 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.

5.4.1 Ideal targets

Young Class 0 and Class I disks seem to have typical sizes around $100\text{\,}$ - $200\text{\,}\mathrm{a}\mathrm{u}$ (Najita & Bergin, 2018, Table 1), so solids can be delivered to the inner regions only until 0.5\text{\,}$-$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}, before the disk is depleted of dust (unless a pressure bump prevents particles from moving towards the star). This means that older disks ($t>$1\text{,}\mathrm{M}\mathrm{y}\mathrm{r}$$) 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:

[TABLE]

could be subject to the back-reaction damping, as a low viscous evolution ($\alpha_{\nu}\lesssim${10}^{-3}\text{,}$$) is required for the dust to affect the gas. In terms of the dimensionless accretion parameter introduced by Rosotti et al. (2017), defined as:

[TABLE]

a disk of age $\tau$ would require $\eta\lesssim 0.1$ for the dust back-reaction to effectively perturb the gas.

5.4.2 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}$ (Equation 15).

At the midplane, where large grains concentrate, the gas motion deviates from the keplerian velocity by:

[TABLE]

where the $\sqrt{\mathrm{St}/\alpha_{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 snowline was $v_{P}\approx$2\text{\times}{10}^{-3}\text{,}$v_{K}$, then the dust back-reaction and the accumulation of water vapor makes the gas orbit at velocities of $\Delta v_{\textrm{g},\theta}\approx$7\text{\times}{10}^{-4}\text{,}$v_{K}$. At $2.7\text{\,}\mathrm{a}\mathrm{u}$, where the snowline is located in our simulations, this correspond to a difference from the keplerian velocity of approximately $\Delta v_{\textrm{g},\theta}\approx$10\text{,}\mathrm{m}\mathrm{/}\mathrm{s}$$.

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\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$. 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\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$, for a disk at a distance of $100\text{\,}\mathrm{a}\mathrm{u}$ and a snowline at $3\text{\,}\mathrm{a}\mathrm{u}$ 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 snowline is located at larger radii.

5.4.3 Shadows casted by dust accumulation

A recent study of Ueda et al. (2019) showed 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\text{\,}\mathrm{a}\mathrm{u}$.

We notice that our accumulation of dust at the snowline 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 (Figure 4). Therefore, we hypothesize that similar shadows could be found in the regions just outside the snowline 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.

5.4.4 Effects of the snowline traffic jam

The fast drift of the icy particles particles and the traffic jam at the snowline results in the accumulation of both small silicate dust and water vapor inside the snowline, even if the effect of the dust back-reaction is ignored (see Figure 9 and 10).

We find that if the initial dust supply is large enough (high $\epsilon_{0}$ and large disk size), then during the early stages of the disk evolution ($t\lesssim$1\text{,}\mathrm{M}\mathrm{y}\mathrm{r}$$) we can expect the material accreted into the star to be rich in silicates and refractory materials carried by the dust (see Figure 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 snowline (see Figure 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), though the coronal emission of neon in young stars could mask some of these abundances (H. M. Günther, private communication).

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 snowline for young protoplanetary disks (Öberg et al., 2011; Booth & Clarke, 2018).

6 Summary

In this study we included the effects of the dust back-reaction on the gas in a model of the water snowline, which is known to act as a concentration point for dust particles due to the change in the fragmentation velocity between silicates and ices, and 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, though the parameter space required for this to happen is limited.

In the vicinity of the snowline, provided that the global dust-to-gas ratio is high ($\varepsilon_{0}\gtrsim 0.03$) and the viscosity low ($\alpha_{\nu}\lesssim${10}^{-3}\text{,}$$), the effects of the dust back-reaction are:

• •

Revert the net gas flux outside the snowline.

• •

Reduce the gas inflow at the snowline 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 snowline.

• •

Concentrate solids at the snowline reaching dust-to-gas ratios of $\epsilon\gtrsim 0.8$.

These effects build up as long as dust is supplied from the outer disk into the snowline, 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\text{,}\mathrm{a}\mathrm{u}$$, we find that dust accumulates only during the first $0.4\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$, and that the perturbation onto the gas has disappeared by the age of $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$.

The high dust-to-gas ratios required to trigger the back-reaction effects, and the traffic jam at the snowline, can result in an enhanced water content in the inner regions, in the accretion onto the star to be enriched with refractory materials and oxygen, and perhaps a shadow to be casted outside the snowline location by the accumulation of dust particles.

Other types of dust traps could present similar behaviors, though each case must be revisited individually to evaluate the magnitude of the perturbation of the back-reaction into the gas velocity.

Acknowledgements.

We would like to thank S. Facchini, H. M. Günther, and G. Rosotti for their comments and discussions during the early stages of this manuscript. We also want to thank the anonymous referee for his/her comments, that improved the extent and clarity of this work. The authors acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 714769, by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Ref no. FOR 2634/1, and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311.

Appendix A Semi-analytical test for back-reaction simulations.

In this section we intend to rewrite the radial velocity of the gas (Equation 14) in a similar way to the standard viscous velocity of Lynden-Bell & Pringle (1974). The viscous velocity and the pressure velocity (Equation 16 and 17) can be rewritten in the following form:

[TABLE]

[TABLE]

with $\gamma_{\nu}=\textrm{d}\mathrm{ln}(\nu\,\Sigma_{\textrm{g}}\,\sqrt{r})/\textrm{d}\mathrm{ln}\,r$ and $\gamma_{P}=\textrm{d}\mathrm{ln}P/\textrm{d}\mathrm{ln}\,r$.

Using these expressions, we can rewrite the gas radial velocity (Equation 14) as the viscous velocity in Equation 43, with the following $\alpha_{\nu}$-equivalent parameter:

[TABLE]

[TABLE]

This means that we can understand the evolution of a gas disk considering dust back-reaction, as a viscous evolution with a modified $\alpha_{\nu}$ value (we discuss the limits of this interpretation in Appendix A.2). From this point we can make further simplifications to develop a semi-analytical test for a back-reaction simulation using a standard viscous evolution model.

Our first simplification is that the surface density and temperature follow a power law profile with $\Sigma\propto r^{-p}$ and $T\propto r^{-q}$, which sets the factor $\gamma_{P}/(3\gamma_{\nu})$ involving the density and temperature gradients to:

[TABLE]

In particular, if the disk is in steady state with $p=1$ and $q=1/2$, then $\gamma_{P}/(3\gamma_{\nu})=-11/6$, and the accretion rate is:

[TABLE]

Now, assuming that the distribution of dust particles has sizes between $0<\mathrm{St}<\mathrm{St}_{\textrm{max}}$, and that $\mathrm{St}^{2}\ll\epsilon$, we can constrain the value of $\alpha_{\textrm{eq}}$ using the single particle approximation for the coefficients $A_{\textrm{single}}$ and $B_{\textrm{single}}$ (Equation 29 and 30). Then the minimum value that $\alpha_{\textrm{eq}}$ can take, given by the largest size particles, is:

[TABLE]

and the maximum value that $\alpha_{\textrm{eq}}$ can take, given by the smallest particles with $\mathrm{St}\approx 0$, is:

[TABLE]

A.1 Setting up a test simulation

From the equivalent viscosity equation (Equation 45) we can set a test to ensure that the back-reaction effects in a numerical simulation are acting according to the theoretical model.

We prepare a test for the code twopoppy that was used throughout the paper (Birnstiel et al. 2012), and also for the code DustPy (Stammler and Birnstiel, in prep.), that solves the Smoluchowski equation for particle growth by sticking and fragmentation of multiple dust species as in Birnstiel et al. (2010), along with the advection-diffusion equations (Equation 1 and 2).

The test disk has the following set-up:

• •

The surface density and temperature have steady state power law profiles with $p=1$ and $q=1/2$.

• •

To enhance the back-reaction damping and obtain obvious deviations from the regular dust-free evolution we set an unrealistic disk with $\epsilon=0.5$.

• •

The fragmentation velocity follows $v_{\textrm{frag}}\propto r^{-q}$ so that the maximum particle size (Equation 9) has a constant value of $\mathrm{St}_{\textrm{max}}=$5\text{\times}{10}^{-3}\text{,}$$.

• •

The viscous turbulence is set to $\alpha_{\nu}=${10}^{-2}\text{,}$$, so that the back-reaction is not strong enough to reverse the accretion of gas.

• •

The dust diffusion is turned off, so that the dust is only advected through the velocity $v_{\textrm{d,r}}$ (Equation 3).

• •

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 $\alpha_{\textrm{eq}}$ (Equation 48). Since in this test case all the particles are small ($\mathrm{St}<\alpha_{\nu}$) 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 Figure 14, after $0.1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ the disk surface density between 5\text{\,}$-$100\text{\,}\mathrm{a}\mathrm{u} remains close to the steady state profile, with a deviation of less than $0.1\%$ relative to its initial value.

Figure 15 shows that the mass accretion rate of the gas in the simulations is $\dot{M}\approx$5.6\text{\times}{10}^{-8}\text{,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{y}\mathrm{r}$$ 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 $\alpha_{\textrm{eq, min}}$ and $\alpha_{\textrm{eq, max}}$ and Equation 48.

In terms of the viscous accretion, the back-reaction effect in our setup is equivalent to reduce the viscous turbulence $\alpha_{\nu}$ to a value of $\alpha_{\textrm{eq}}\approx 0.57\,\alpha_{\nu}$.

Both twopoppy and DustPy deliver similar results, with a relative difference of roughly $5\%$ in the $\alpha_{\textrm{eq}}$ and $\dot{M}$ 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 Equation 46 using the $\alpha_{\textrm{eq}}$ parameter (Equation 45), the global disk evolution will still differ from a regular viscous evolution (unless $\alpha_{\textrm{eq}}\propto\alpha_{\nu}$), as the value of $\gamma_{\nu}$ does not depend on the slope of $\alpha_{\textrm{eq}}$.

In particular, the back-reaction effects cannot be treated as a viscous process if $\mathrm{St}\,\epsilon/(\epsilon+1)\gtrsim\alpha_{\nu}$ (Dipierro et al. 2018). In this case the back-reaction push becomes more important than the inward viscous transport, and results in negative equivalent $\alpha_{\textrm{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 ($\gamma_{\nu}<0$, $v_{\nu}>0$). In these regions the dust back-reaction 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 section 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_{\nu}(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_{\nu}$ 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_{\nu}$, $v_{P}$ are:

[TABLE]

[TABLE]

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 (Equation 22 and 28) and check how our results are affected.

Figure 16 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 (Equation 51 and 52). For the simulations with $\varepsilon\geq 0.03$, only at the snowline location the dust-to-gas ratio is more spread over 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 well the radial velocity profile calculated with the Takeuchi & Lin (2002) model (Figure 17), except in a narrow region beyond the snowline 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 for gas and dust interaction.

The first regime is the one in which the particles are small ($\mathrm{St}<\alpha_{\nu}$). 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 $\bar{v}_{\textrm{g},r}\approx\bar{A}v_{\nu}$.

The second regime is the case in which the particles are large ($\mathrm{St}\gtrsim\alpha_{\nu}$). 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_{\textrm{g},r}(z\lesssim h_{\textrm{d}})\approx 2\bar{B}(z\lesssim h_{\textrm{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_{\textrm{g},r}(z\gtrsim h_{\textrm{d}})\approx v_{\nu}$. Then, the velocity derived from the net flux (Equation 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\approx$3\text{,}\mathrm{a}\mathrm{u}$$), however this seem 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.

Appendix C Parameter space exploration

Here we extend the parameter space described in Table 1 to show the effect of different turbulent viscosity parameters $\alpha_{\nu}$ in the global disk evolution (Figure 18), along with the dust-to-gas ratio profile, the $\textrm{H}_{2},\textrm{He}$ mass fraction profiles, and the gas accretion rate evolution (Figure 19).

We can summarize the plots in a few remarks:

• •

Lower $\alpha_{\nu}$ values leads to higher dust-to-gas ratios, as the dust is more concentrated towards the snowline and spreads more slowly to the inner boundary.

• •

Higher $\varepsilon_{0}$ values lead to a stronger perturbation onto the gas surface density, but only if the turbulent viscosity is low enough ($\alpha_{\nu}\leq${10}^{-3}\text{,}$$).

• •

For an initial $\varepsilon_{0}\geq 0.03$ and $\alpha_{\nu}\leq${10}^{-3}\text{,}$$, the dust accumulates both inside and outside the snowline, always reaching dust-to-gas ratios above $\epsilon\geq 0.8$, even if the turbulence is high.

• •

For $\varepsilon_{0}\geq 0.03$ and $\alpha_{\nu}\leq${10}^{-3}\text{,}$$, the $\textrm{H}_{2},\textrm{He}$ mass fraction is reduced to values between $0.2-0.65$ inside the snowline, though the distribution of water vapor depends on the turbulent viscosity $\alpha_{\nu}$.

• •

For the case with $\varepsilon_{0}=0.01$ and $\alpha_{\nu}=${10}^{-4}\text{,}$$, the dust concentration is enhanced in a narrow region inside the snowline because of the low viscosity.

• •

The accretion rate can be reduced down to $30\%$ of its initial value depending, on the simulation parameters.

• •

In the high turbulence case ($\alpha_{\nu}=${10}^{-2}\text{,}) the dust concentration inside the snowline is reduced to $\epsilon\approx 0.01-0.1$, though it can reach to $\epsilon=0.02-0.2$ outside the snowline ($r\approx$3\text{\,}$-$4\text{\,}\mathrm{a}\mathrm{u}) due to the recondensation of water vapor.

From these plots we can extract that a global dust-to-gas ratio $\varepsilon_{0}\geq 0.03$ and a viscous turbulence $\alpha_{\nu}\leq${10}^{-3}\text{,}$$ are required for the dust back-reaction to perturb the gas surface density, deplete the the inner regions from hydrogen-helium ($\Sigma_{\textrm{H}_{2}}/\Sigma_{\textrm{g}}\lesssim 0.5$), reach high concentrations of dust inside and outside the snowline ($\epsilon\gtrsim 0.5$), and finally, to damp the accretion rate ($\dot{M}/\dot{M}_{0}\approx 0.3-0.5$).

Appendix D Criterion for gap opening through dust back-reaction.

In this section we will derive a criterion to determine if the dust back-reaction can clear a gap in the gas, based on the transport of angular momentum and the global disk properties.

The condition to clear a gap, is that the clearing timescale $t_{\textrm{clear}}$ must be shorter than the viscous timescale $t_{\nu}$:

[TABLE]

Now we proceed to derive $t_{\textrm{clear}}$ from the exchange of angular momentum between the dust and gas.

The angular momentum of a parcel of material with mass $m$, at a radius $r$, and with orbital velocity $v_{K}=\sqrt{GM/r}$ is:

[TABLE]

The angular momentum required to transport a ring of material from a radius $r$ to $r_{1}$ is:

[TABLE]

where $v_{0}$ is the keplerian velocity at radius $r_{0}$, and the mass of a gas ring is $\textrm{d}m=2\pi r\Sigma_{\textrm{g}}\textrm{d}r$.

Assuming that the gas surface density is $\Sigma_{\textrm{g}}=\Sigma_{0}(r_{0}/r)$, then the total angular momentum required to clear a gap between $r_{0}$ and $r_{1}$ is given by:

[TABLE]

Solving the integral we obtain:

[TABLE]

From Equation 55 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:

[TABLE]

where the accretion rate of dust at $r_{1}$ is:

[TABLE]

where we will assume an uniform dust-to-gas ratio, using $\Sigma_{\textrm{d}}=\varepsilon\,\Sigma_{\textrm{g}}$.

Then, considering only the drifting component of the dust velocity (see Equation 3) we obtain that:

[TABLE]

where have taken the limit of small particles ($\textrm{St}\ll 1$), and expanded the expression for the pressure velocity given in Equation 44 using $h/r=c_{s}/v_{K}$.

Now we can write the gap opening timescale as:

[TABLE]

with $C(x)$ a function of the ratio $x=r_{1}/r_{0}$:

[TABLE]

We apply this formula to our simulations, using a scale height of $h/r=0.05$, $\gamma_{P}=2.75$, and study the time required to clear a gap between $r_{0}=$2.5\text{,}\mathrm{a}\mathrm{u} (which is approximately the location of the snowline) and $r_{1}=$10\text{\,}\mathrm{a}\mathrm{u}, and find:

[TABLE]

Now, we only need to compare it with the viscous timescale, which can be understood in this case as the time necessary to close the gap:

[TABLE]

which at the snowline location of $r_{0}=$2.5\text{,}\mathrm{a}\mathrm{u}$$ give us a time of:

[TABLE]

To conclude we can derive from Equation 53, (61) and (64) a condition on $\alpha$, $\varepsilon$ and St to see if the dust back-reaction can clear a gap:

[TABLE]

where we have taken $r_{0}=$2.5\text{,}\mathrm{a}\mathrm{u} and $r_{1}=$10\text{\,}\mathrm{a}\mathrm{u}. Notice that this condition is similar to $2Bv_{P}\gtrsim Av_{\nu}$ (Equation 14), which indicates if gas motion is locally dominated by the dust back-reaction.

To conclude, looking at the values of $t_{\textrm{clear}}$ and $t_{\nu}$, it is easy to see why some disks create a gap-like perturbation in Figure 18. The disks with $\alpha=${10}^{-2}\text{,} have viscous timescales that are too short in comparison with the clearing timescale by an order of magnitude, while the disks with $\alpha=${10}^{-4}\text{\,} are easily dominated by the dust back-reaction, provided that dust is delivered for enough time to complete a clearing timescale.