© M. Gárate et al. 2020

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 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 back-reaction 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 back-reaction 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 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): $$\frac{\partial }{\partial t}\left(r{\Sigma }_{\text{g}}\right)+\frac{\partial }{\partial r}(r{\Sigma }_{\text{g}}{v}_{\text{g},\text{r}})=0,$$(1) $$\frac{\partial }{\partial t}\left(r{\Sigma }_{\text{d}}\right)+\frac{\partial }{\partial r}(r{\Sigma }_{\text{d}}{v}_{\text{d},\text{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,$$(2)

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.

2.1 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): $$v_{\text{d},r}=\frac{1}{1+{\text{St}}^2}{v}_{\text{g},r}+\frac{2\text{St}}{1+{\text{St}}^2}\Delta v_{\text{g},\theta },$$(3) $$\Delta v_{\text{d},\theta }=\frac{1}{1+{\text{St}}^2}\Delta v_{\text{g},\theta }-\frac{\text{St}}{2(1+{\text{St}}^2)}{v}_{\text{g},r},$$(4)

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: $$\text{St}=t_{\text{stop}}{\Omega }_{\text{K}},$$(5)

where Ω_K is the Keplerian angular velocity and t_stop is: $$t_{\text{stop}}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_{\text{s}}}{{\rho }_{\text{g}}}\frac{a}{c_{\text{s}}},$$(6)

with a the particle size, ρ_s the material density of the solids, ρ_g the gas density. The isothermal sound speed c_s is: $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu m_{\text{H}}}},$$(7)

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: $$\text{St}=\frac{\pi }{2}\frac{a{\rho }_{\text{s}}}{{\Sigma }_{\text{g}}}.$$(8)

The size of the particles is not static in time (Birnstiel et al. 2010, 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: $${\text{St}}_{\text{frag}}=\frac{1}{3}\frac{v_{\text{frag}}^2}{{\alpha }_t{c}_{\text{s}}^2},$$(9)

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: $${\text{St}}_{\text{drift}}={\left|\frac{\text{dln\hspace{0.17em}}P}{\text{dln\hspace{0.17em}}r}\right|}^{-1}\frac{v_{\text{K}}^2}{c_{\text{s}}^2}\varepsilon ,$$(10)

with v_K the Keplerian velocity, and P the isothermal gas pressure at the midplane: $$P=\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{h}_{\text{g}}}{c}_{\text{s}}^2,$$(11)

with the gas scale height h_g = c_s∕Ω_K.

Additionally, we assume that dust diffuses with $$D_{\text{d}}=\frac{\nu }{(1+\varepsilon )},$$(12)

with ɛ = Σ_d∕Σ_g the vertically integrated dust-to-gas ratio, and ν the turbulent viscosity of the gas (Shakura & Sunyaev 1973): $$\nu ={\alpha }_{\nu }{c}_{\text{s}}^2{\Omega }_{\text{K}}^{-1},$$(13)

controlled by the viscous turbulence parameter α_ν.

We notice that, as in Carrera et al. (2017), our model consi- ders two different turbulence parameters: α_t for the dustturbulence (that controls the dust fragmentation, Eq. (9)) and α_ν for the viscous turbulence (that controls the gas viscosity, Eq. (13)).

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

2.2 Gas dynamics

The gas velocities, considering the dust back-reaction onto the gas, have the following form: $$v_{\text{g},r}=A{v}_{\nu }+2B{v}_{\text{P}},$$(14) $$\Delta v_{\text{g},\theta }=-A{v}_{\text{P}}+\frac{1}{2}B{v}_{\nu }.$$(15)

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): $$v_{\nu }=-\frac{3}{{\Sigma }_{\text{g}}\sqrt{r}}\frac{\partial }{\partial r}(\nu {\Sigma }_{\text{g}}\sqrt{r}).$$(16)

Similarly, if there is no dust, the gas orbits at sub-Keplerian speeds due to the pressure support, this pressure velocity is given by: $$v_{\text{P}}=-\frac{1}{2}{\left(\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{h}_{\text{g}}}{\Omega }_{\text{K}}\right)}^{-1}\frac{\partial P}{\partial r}.$$(17)

The back-reaction coefficients are defined as follows: $$A=\frac{X+1}{Y^2+{(X+1)}^2},$$(18) $$B=\frac{Y}{Y^2+{(X+1)}^2},$$(19)

where X and Y are the following sums defined by (Tanaka et al. 2005; Okuzumi et al. 2012; Dipierro et al. 2018): $$X={\sum }_m\frac{1}{1+\text{St}{(m)}^2}\varepsilon (m),$$(20) $$Y={\sum }_m\frac{\text{St}(m)}{1+\text{St}{(m)}^2}\varepsilon (m),$$(21)

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³ρ_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.

2.2.1 Effect of the vertical structure on the net mass flux

The corrected gas radial velocity ${\overline{v}}_{\text{g},\text{r}}$ can be obtained from the net mass flux, which is defined as: $${\Sigma }_{\text{g}}{\overline{v}}_{\text{g},r}={\displaystyle {\int }_{-\infty }^{+\infty }{\rho }_{\text{g}}}(z)v_{\text{g},r}(z)\text{d}z,$$(22)

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: $${\rho }_{\text{g}}(z)=\frac{{\Sigma }_{\text{g}}}{\sqrt{2\pi }{h}_{\text{g}}}\exp \left(-\frac{z^2}{2h_{\text{g}}^2}\right),$$(23) $${\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),$$(24)

where the vertical scale height of the dust species with mass m is defined in Birnstiel et al. (2010) as: $$h_{\text{d}}(m)=h_{\text{g}}\cdot \min \left(1,\sqrt{\frac{{\alpha }_t}{\min (\text{St},1/2)(1+{\text{St}}^2)}}\right).$$(25)

From these profile we can obtain the dust-to-gas ratio of every particle species at height z with: $$\varepsilon (z,m)=\frac{{\rho }_{\text{d}}(z,m)}{{\rho }_{\text{g}}(z)},$$(26)

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: $${\overline{v}}_{\text{g},r}=\frac{1}{{\Sigma }_{\text{g}}}{\displaystyle {\int }_{-\infty }^{+\infty }{\rho }_{\text{g}}}(z)\left(A(z)v_{\nu }+2B(z)v_{\text{P}}\right)\text{d}z=\overline{A}{v}_{\nu }+2\overline{B}{v}_{\text{P}},$$(27)

where Ā and $\overline{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: $${\Sigma }_{\text{d}}(m){\overline{v}}_{\text{d},r}(m)={\displaystyle {\int }_{-\infty }^{+\infty }{\rho }_{\text{d}}}(z,m)v_{\text{d},r}(z,m)\text{d}z.$$(28)

The gas and dust velocities derived from the net mass flux (${\overline{v}}_{\text{g},r}$ and ${\overline{v}}_{\text{d},r}$) are used to transport the gas and dust in the advection-diffusion Eqs. (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: (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): $$A_{\text{single}}=\frac{\varepsilon +1+{\text{St}}^2}{{(\varepsilon +1)}^2+{\text{St}}^2},$$(29) $$B_{\text{single}}=\frac{\varepsilon \text{St}}{{(\varepsilon +1)}^2+{\text{St}}^2}.$$(30)

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

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 Awe 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 α_ν 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)).

2.3 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: $$\mu =({\Sigma }_{{\text{H}}_2}+{\Sigma }_{\text{vap}}){\left(\frac{{\Sigma }_{{\text{H}}_2}}{{\mu }_{{\text{H}}_2}}+\frac{{\Sigma }_{\text{vap}}}{{\mu }_{\text{vap}}}\right)}^{-1},$$(31)

where ${\mu }_{{\text{H}}_2}=2.3$ and μ_vap = 18 are respectively the mean molecular weights of the hydrogen-helium mixture and the water vapor, and ${\Sigma }_{{\text{H}}_2}={\Sigma }_{\text{g}}-{\Sigma }_{\text{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: $${\rho }_{\text{s}}=({\Sigma }_{\text{sil}}+{\Sigma }_{\text{ice}}){\left(\frac{{\Sigma }_{\text{sil}}}{{\rho }_{\text{sil}}}+\frac{{\Sigma }_{\text{ice}}}{{\rho }_{\text{ice}}}\right)}^{-1},$$(32)

where ρ_sil = 3 g cm⁻³ and ρ_ice = 1 g cm⁻³ 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⁻¹ (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⁻¹ (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_{\text{eq}}=P_{\text{eq},\text{0}}\exp (-A/T),$$(33)

with P_eq,0 = 1.14 × 10¹³g cm⁻¹s⁻² 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: $$P_{\text{vap}}=\frac{{\Sigma }_{\text{vap}}}{\sqrt{2\pi }{h}_{\text{g}}}\frac{k_{\text{B}}T}{{\mu }_{\text{vap}}{m}_{\text{H}}}.$$(34)

When the water vapor pressure is below this threshold (P_vap < P_eq) the ice evaporates into vapor as follows: $$\Delta {\Sigma }_{\text{vap}}=\min \left(\sqrt{2\pi }{h}_{\text{g}}\frac{{\mu }_{\text{vap}}{m}_{\text{H}}}{k_{\text{B}}T}(P_{\text{eq}}-P_{\text{vap}}),{\Sigma }_{\text{ice}}\right),$$(35)

and vice-versa, if the vapor pressure is higher then it recondenses into ice with: $$\Delta {\Sigma }_{\text{ice}}=\min \left(\sqrt{2\pi }{h}_{\text{g}}\frac{{\mu }_{\text{vap}}{m}_{\text{H}}}{k_{\text{B}}T}(P_{\text{vap}}-P_{\text{eq}}),{\Sigma }_{\text{vap}}\right),$$(36)

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).

3 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.

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 (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.

3.2 Disk initial conditions

The gas surface density and temperature profile are defined by the following power laws: $${\Sigma }_{\text{g}}(r)={\Sigma }_0{\left(\frac{r}{r_0}\right)}^{-p},$$(37) $$T(r)=T_0{\left(\frac{r}{r_0}\right)}^{-q},$$(38)

with r₀ = 1 au, Σ₀ = 1000 g cm⁻², T₀ = 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 ε₀ such that Σ_d = ε₀Σ_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⁻³, and an initial size of a₀ = 1 μm.

3.3 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 accountdust settling (Eqs. (22) and (28)), we construct a local vertical grid at every radius with n_z = 300 points, logarithmically spaced between 10⁻⁵ and 10 h_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 ε₀, 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 ε₀ summarized in Table 1.

For the sake of clarity, through the paper, we use a single value for the viscous turbulence, with α_ν = 10⁻³. This turbulence is low enough for the back-reaction effects to start affecting the gas dynamics (i.e., the term $2\overline{B}{v}_{\text{P}}$ becomes comparable to Āv_ν in the gas velocity, Eq. (14)).

For completeness, in Appendix C we further extend our parameter space¹ to include different values for the viscous turbulence α_ν, although for simplicity, we keep the dust turbulence constant, with α_t = 10⁻³.

4 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 ε₀” 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 ε₀” 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 ε₀ ” simulation (Fig. 2, bottom panel). As in the “Mid ε₀” 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 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 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 ε₀ ” 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 ε₀” 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, as time passes, the dust should spread more evenly towards the inner regions.

The most extreme case is the “High ε₀” 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).

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 upto St ~ 10−4. Outside the snow line it can reach values of St ~ 10−2−10−1.

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.

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.

Fig. 4 Dust-to-gas ratio radial profile for the three simulations after 0.4 Myr. The simulations with a high global dust-to-gas ratio (ε0 ≥ 0.03), shown an enhanced dust accumulation outside the snow line, reaching ɛ ≈ 0.8−2.0.

Fig. 5 Radial gas velocities after 0.4 Myr (solid lines), and initial viscous velocity (dashed line). Outside the snow line, the dust back-reaction can stop, and even reverse, the gas flux for the simulations with ε0 ≥ 0.03.

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 (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 ε₀ ” 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_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⁻⁹M_⊙ yr⁻¹ to 6.8 × 10⁻⁹M_⊙ yr⁻¹ (Fig. 6). Once the dust supply is depleted, the accretion rate should return to its steady state value.

In the “High ε₀” 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⁻⁴), so the gas velocity is dominated by the term Āv_ν, which corresponds to the viscous velocity damped by a factor of $\overline{A}\approx {(1+\varepsilon )}^{-1}$. In the outer region (r > 2.5 au) where the particles are large (St ≳×10⁻²), the velocity is dominated by the pressure velocity term $2\overline{B}{v}_{\text{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⁻⁹ M_⊙ yr⁻¹. 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⁻⁹−×10⁻⁸ 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 ε₀” simulation is consistently in between the “Low ε₀” and “High ε₀” cases, with the gas flux practically frozen (v_g,r ≈ 0) in the outer regions (r > 4 au). 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), wherethe 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.

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. 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\overline{B}{v}_{\text{P}}$ (blue) (see Eq. (14)). In the inner regions, the pushing term $2\overline{B}{v}_{\text{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\overline{B}{v}_{\text{P}}$ overcomes the viscous evolution, and pushes gas against the pressure gradient.

4.2 Depletion of H₂ 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₂, 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₂, 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 ɛ₀, while the dust back-reaction affects how it is distributed.

Figure 8 shows that even in the “Low ε₀ ” case, the mass fraction of H₂, He is reduced to a 90%.

For the “Mid ε₀” and “High ε₀” 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₂, 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₂, He mixture willbe 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 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.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.03, we observe thatwithout 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.03, the dust back-reaction 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.

Fig. 8 H2, 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 H2, He into the inner regions.

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

In the “High ε₀” 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.05 as in the “High ε₀” case, but this time starting with a self-similar profile (Lynden-Bell & Pringle 1974), following: $${\Sigma }_{\text{g}}(r)={\Sigma }_0{\left(\frac{r}{r_0}\right)}^{-p}\exp (-r/r_{\text{c}}),$$(39)

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 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 dustis delivered at the snow line, the accretion rate of gas is damped, and the mass fraction of the H₂, 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 the accretion rate reaches its minimum of 3.0 × 10⁻⁹M_⊙ yr⁻¹, where only 60% of the accretion flow corresponds to H₂, 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.

Fig. 9 Comparison of the surface density profiles when the back-reaction 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.

Fig. 10 Comparison of the dust-to-gas ratio profiles when the back-reaction 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.

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 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 ($\overline{A}{v}_{\nu }~2\overline{B}{v}_{\text{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\overline{B}{v}_{\text{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.03 and a low viscous turbulence of α_ν ≲ 10⁻³ (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).

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.

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. 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 H2, 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.

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 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 back-reaction 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.

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 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⁻⁴) 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 × 10⁻⁹M_⊙ yr⁻¹, only 6.0 × 10⁻¹⁰M_⊙ yr⁻¹ 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⁻⁸M_⊙ yr⁻¹), 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⁻²) 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⁻¹⁰M_⊙ yr⁻¹ for the simulation with ε₀ = 0.03 (in which the reduced inflow only occurs above 0.7h_g), and to Ṁ = 6 × 10⁻¹²M_⊙ yr⁻¹ for the simulation with ε₀ = 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⁻³ 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.03, the H₂ + 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.05, all our results hold.

Fig. 13 Top: mass flux for the simulation with ε0 = 5%, in the radialand 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.

5.4 Observational implications

The perturbation caused by the dust back-reaction at the snow line is only effective if the viscous turbulence is low, the dust-to-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.

5.4.1 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 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: $$\text{Dot}M\lesssim \times {10}^{-8}{M}_{\odot }{\text{yr}}^{-1}\left(\frac{M_{\text{disk}}}{0.1}{M}_{\odot }\right){\left(\frac{r_{\text{c}}}{100\text{au}}\right)}^{-1}{\left(\frac{M_{\text{star}}}{M_{\odot }}\right)}^{-1/2}\left(\frac{T_0}{300 ,\text{K}}\right),$$(40)

could be subject to the back-reaction damping, as a low viscous evolution (α_ν ≲ 10⁻³) is required for the dust to affect the gas. In terms of the dimensionless accretion parameter introduced by Rosotti et al. (2017), defined as: $$\eta =\frac{\tau \text{ Dot}M}{M_{\text{disk}}},$$(41)

a disk of age τ would require η ≲ 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 (Eq. (15)).

At the midplane, where large grains concentrate, the gas motion deviates from the Keplerian velocity by: $$\Delta v_{\text{g},\theta }\approx -\frac{v_{\text{P}}}{1+\max (1,\sqrt{\text{St}/{\alpha }_t})\cdot \varepsilon },$$(42)

where the $\sqrt{\text{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 snow line was v_P ≈ 2 × 10⁻³v_K, then the dust back-reaction and the accumulation of water vapor makes the gas orbit at velocities of Δv_g,θ ≈ 7 × 10⁻⁴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⁻¹.

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⁻¹. 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.

5.4.3 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.

5.4.4 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 ɛ₀ 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).

6 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 dust-to-gas ratio is high (ε₀ ≳ 0.03) and the viscosity low (α_ν ≲ 10⁻³), 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.

The high dust-to-gas ratios required to trigger the back-reaction 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.

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

Fig. A.1 Left: initial and final surface density of the test simulations. If the disk evolution with back-reaction is equivalent to a regular viscous evolution, the steady state should be maintained through the simulation. Right: surface density residuals relative to the initial state. After 0.1 Myr of evolution, the simulations deviate by less a value of 0.05% from the steady state profile.

In this section, we intend to rewrite the radial velocity of the gas (Eq. (14)) in a similar way to the standard viscous velocity of Lynden-Bell & Pringle (1974). The viscous velocity and the pressure velocity (Eqs. (16) and (17)) can be rewritten in the following form: $$v_{\nu }=-3{\alpha }_{\nu }\frac{c_{\text{s}}^2}{v_{\text{K}}}{\gamma }_{\nu },$$(A.1) $$v_{\text{P}}=-\frac{1}{2}\frac{c_{\text{s}}^2}{v_{\text{K}}}{\gamma }_P,$$(A.2)

with ${\gamma }_{\nu }=\text{dln}(\nu {\Sigma }_{\text{g}}\sqrt{r})/\text{dln}r$ and γ_P = dlnP∕dln r.

Using these expressions, we can rewrite the gas radial velocity (Eq. (14)) as the viscous velocity in Eq. (A.1), with the following α_ν-equivalent parameter: $${\alpha }_{\text{eq}}=A{\alpha }_{\nu }+\frac{{\gamma }_P}{3{\gamma }_{\nu }}B,$$(A.3) $$v_{\text{g}}=-3{\alpha }_{\text{eq}}\frac{c_{\text{s}}^2}{v_{\text{K}}}{\gamma }_{\nu }.$$(A.4)

This means that we can understand the evolution of a gas disk by considering dust back-reaction as a viscous evolution with a modified α_ν 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 Σ ∝r^−p and T ∝r^−q, which sets the factor γ_P∕(3γ_ν) involving the density and temperature gradients to: $$\frac{{\gamma }_P}{3{\gamma }_{\nu }}=-\frac{2p+q+3}{6(2-p-q)}.$$(A.5)

In particular, if the disk is in steady state with p = 1 and q = 1∕2, then γ_P∕(3γ_ν) = −11∕6, and the accretion rate is: $$\text{Dot}M=3\pi {\alpha }_{\text{eq}}\frac{c_{\text{s}}^2}{{\Omega }_{\text{K}}}{\Sigma }_{\text{g}}.$$(A.6)

Now assuming that the distribution of dust particles has sizes between 0 < St < St_max and that St² ≪ ɛ, we can constrain the value of α_eq using the single particle approximation for the coefficients A_single and B_single (Eqs. (29) and (30)). Then the minimum value that α_eq can take, as given by the largest size particles, is: $${\alpha }_{\text{eq},\text{min}}\approx \frac{{\alpha }_{\nu }}{\varepsilon +1}-\frac{11}{6}\frac{\varepsilon {\text{St}}_{\text{max}}}{{(\varepsilon +1)}^2},$$(A.7)

and the maximum value that α_eq can take, given by the smallest particles with St ≈ 0, is: $${\alpha }_{\text{eq},\text{max}}\approx \frac{{\alpha }_{\nu }}{\varepsilon +1}.$$(A.8)

A.1 Setting up a test simulation

From the equivalent viscosity equation (Eq. (A.3)), 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 is used throughout the paper (Birnstiel et al. 2012) and also for the code DustPy (Stammler and Birnstiel, in prep.) which 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 (Eqs. (1) and (2)).

The test disk has the following set-up:

• 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 ɛ = 0.5;

• the fragmentation velocity follows v_frag ∝ r^−q, so that the maximum particle size (Eq. (9)) has a constant value of St_max = 5 × 10⁻³;

• the viscous turbulence is set to α_ν = ×10⁻², 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_d,r (Eq. (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 α_eq (Eq. (A.6)). Since in this test case all the particles are small (St   < α_ν) 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.

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.

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⁻⁸M_⊙ yr⁻¹ 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 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 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: $$v_{\nu }(z)=\frac{\nu }{2r}\left(6p+q-3+(5q-9){\left(\frac{z}{h_{\text{g}}}\right)}^2\right),$$(B.1) $$v_{\text{P}}(z)=v_{\text{K}}{\left(\frac{h_{\text{g}}}{r}\right)}^2\left(p+\frac{q+3}{2}+\frac{q-3}{2}{\left(\frac{z}{h_{\text{g}}}\right)}^2\right),$$(B.2)

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.

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).

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.

Fig. B.2 Gas radial velocity profiles (derived from the net mass flux, Eq. (22)) using our simple vertical structure approximation vs. Takeuchi & Lin (2002) model (Eqs. (B.1) and (B.2)).

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 ${\overline{v}}_{\text{g},r}\approx \overline{A}{v}_{\nu }$.

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_{\text{g},r}(z\lesssim h_{\text{d}})\approx 2\overline{B}(z\lesssim h_{\text{d}})v_{\text{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.

Appendix C Parameter space exploration

Here we extend the parameter space described in Table 1 to show the effect of different turbulent viscosity parameters α_ν in the global disk evolution (Fig. C.1), along with the dust-to-gas ratio profile, the H₂, He mass fraction profiles, and the gas accretion rate evolution (Fig. C.2).

Fig. C.1 Surface density of gas (red) and dust (blue) after 0.4 Myr (solid lines) for different values of ε0 and αν. Initial conditions marked with dashed lines. The snow line location is marked with dotted lines.

Fig. C.2 Top row: dust-to-gas ratio radial profile. Middle row: H2, He mass fraction radial profile. Bottom row: accretion rate time evolution (divided by the initial steady state accretion Ṁ0). All for different values of ε0 and αν. The value of 0.5 is marked with a dotted line in every plot.

We can summarize the plots in a few remarks:

• lower α_ν values lead to higher dust-to-gas ratios as the dust is more concentrated towards the snow line and spreads more slowly to the inner boundary;

• higher ε₀ values lead to a stronger perturbation onto the gas surface density but only if the turbulent viscosity is low enough (α_ν ≤ 10⁻³);

• for an initial ε₀ ≥ 0.03 and α_ν ≤ 10⁻³, the dust accumulates both inside and outside the snow line, always reaching dust-to-gas ratios above ɛ ≥ 0.8, even if the turbulence is high;

• for ε₀ ≥ 0.03 and α_ν ≤ 10⁻³, the H₂,  He mass fraction is reduced to values between 0.2 and 0.65 inside the snow line, although the distribution of water vapor depends on the turbulent viscosity α_ν;

• for the case with ε₀ = 0.01 and α_ν = ×10⁻⁴, the dust concentration is enhanced in a narrow region inside the snow line 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 (α_ν = ×10⁻²) the dust concentration inside the snow line is reduced to ɛ ≈ 0.01−0.1, although it can reach ɛ = 0.02−0.2 outside the snow line (r ≈ 3−4 au) due to the recondensation of water vapor.

From these plots we can extract that a global dust-to-gas ratio ε₀ ≥ 0.03 and a viscous turbulence α_ν ≤ 10⁻³ are required for the dust back-reaction to perturb the gas surface density, deplete the inner regions from hydrogen-helium (${\Sigma }_{{\text{H}}_2}/{\Sigma }_{\text{g}}\lesssim 0.5$), reach high concentrations of dust inside and outside the snow line (ɛ ≳ 0.5), and, finally, to damp the accretion rate (Ṁ∕Ṁ₀ ≈ 0.3−0.5).

Appendix D Criterion for gap opening through dust back-reaction

In this section, we derive a criterion to determine whether 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_clear must be shorter than the viscous timescale t_ν: $$t_{\text{clear}}<t_{\nu }.$$(D.1)

Now we proceed to derive t_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_{\text{K}}=\sqrt{GM/r}$ is: $$J=v_{\text{K}}rm.$$(D.2)

The angular momentum required to transport a ring of material from a radius r to r₁ is: $$\text{d}J=\left(\sqrt{\frac{r_1}{r_0}}-\sqrt{\frac{r}{r_0}}\right)v_0{r}_0\text{d}m,$$(D.3)

where v₀ is the Keplerian velocity at radius r₀, and the mass of a gas ring is dm = 2πrΣ_gdr.

Assuming that the gas surface density is Σ_g = Σ₀(r₀∕r), then the total angular momentum required to clear a gap between r₀ and r₁ is given by: $$J_{\text{clear}}={\displaystyle \int \text{d}J}=2\pi v_0{r}_0^2{\Sigma }_0{\displaystyle {\int }_{r_0}^{r_1}\sqrt{\frac{r_1}{r_0}}}-\sqrt{\frac{r}{r_0}}\text{d}r.$$(D.4)

Solving the integral we obtain: $$J_{\text{clear}}=2\pi v_0{r}_0^3{\Sigma }_0\cdot \left(\frac{1}{3}{\left(\frac{r_1}{r_0}\right)}^{3/2}-{\left(\frac{r_1}{r_0}\right)}^{1/2}+\frac{2}{3}\right).$$(D.5)

From Eq. (D.3) we can also infer that the dust drifting from a radius r₁ to r₀ loses angular momentum (and delivers it to the gas) at a rate of: $$\text{Dot}{J}_{\text{drift}}=v_0{r}_0\text{ Dot}{M}_{\text{d}}(r_1)\left({\left(\frac{r_1}{r_0}\right)}^{1/2}-1\right),$$(D.6)

where the accretion rate of dust at r₁ is: $$\text{Dot}{M}_{\text{d}}(r_1)=2\pi r_1{\Sigma }_{\text{d}}(r_1)v_{\text{d},\text{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: $$v_{\text{d},\text{r}}(r_1)=2\text{St}{v}_{\text{P}}(r_1)=-\text{St}{\gamma }_{\text{P}}{\left(\frac{h}{r}\right)}^2{v}_0{\left(\frac{r_1}{r_0}\right)}^{-1/2},$$(D.8)

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: $$t_{\text{clear}}=\frac{J_{\text{clear}}}{\left|{\dot{J}}_{\text{drift}}\right|}=C\left(\frac{r_1}{r_0}\right){\left(\frac{h}{r}\right)}^{-2}{\epsilon }^{-1}{\text{St}}^{-1}{v}_0^{-1}{r}_0,$$(D.9)

with C(x) a function of the ratio x = r₁∕r₀: $$C(x)=\frac{x^2-3x+2\sqrt{x}}{3{\gamma }_{\text{P}}(\text{}\sqrt{x}-1)}.$$(D.10)

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₀ = 2.5 au (which is approximately the location of the snow line) and r₁ = 10 au, and find: $$t_{\text{clear}}\approx {\left(\frac{\epsilon }{0.01}\right)}^{-1}{\left(\frac{\text{St}}{0.03}\right)}^{-1}0.84\text{Myr}.$$(D.11)

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: $$t_{\nu }=r_0^2/\nu =\frac{r_0}{v_0\alpha }{\left(\frac{h}{r}\right)}^{-2},$$(D.12)

which, at the snow line location of r₀ = 2.5 au, gives us a time of: $$t_{\nu }={\left(\frac{\alpha }{{10}^{-3}}\right)}^{-1}0.25\text{Myr}.$$(D.13)

To conclude, 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: $$\frac{\epsilon \text{St}}{\alpha }\gtrsim C(r_1/r_0=4)\approx 1,$$(D.14)

where we have taken r₀ = 2.5 au and r₁ = 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⁻² have viscous timescales that are too short in comparison with the clearing timescale by an order of magnitude, while the disks with α = ×10⁻⁴ are easily dominated by the dust back-reaction, provided that dust is delivered for enough time to complete a clearing timescale.

References

All Tables

All Figures

	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 upto St ~ 10−4. Outside the snow line it can reach values of St ~ 10−2−10−1.
In the text

	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.
In the text

	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 the text

	Fig. 4 Dust-to-gas ratio radial profile for the three simulations after 0.4 Myr. The simulations with a high global dust-to-gas ratio (ε0 ≥ 0.03), shown an enhanced dust accumulation outside the snow line, reaching ɛ ≈ 0.8−2.0.
In the text

	Fig. 5 Radial gas velocities after 0.4 Myr (solid lines), and initial viscous velocity (dashed line). Outside the snow line, the dust back-reaction can stop, and even reverse, the gas flux for the simulations with ε0 ≥ 0.03.
In the text

	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.
In the text

Fig. 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\overline{B}{v}_{\text{P}}$ (blue) (see Eq. (14)). In the inner regions, the pushing term $2\overline{B}{v}_{\text{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\overline{B}{v}_{\text{P}}$ overcomes the viscous evolution, and pushes gas against the pressure gradient.

In the text

	Fig. 8 H2, 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 H2, He into the inner regions.
In the text

	Fig. 9 Comparison of the surface density profiles when the back-reaction 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.
In the text

	Fig. 10 Comparison of the dust-to-gas ratio profiles when the back-reaction 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.
In the text

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.

In the text

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 H2, 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.

In the text

	Fig. 13 Top: mass flux for the simulation with ε0 = 5%, in the radialand 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.
In the text

	Fig. A.1 Left: initial and final surface density of the test simulations. If the disk evolution with back-reaction is equivalent to a regular viscous evolution, the steady state should be maintained through the simulation. Right: surface density residuals relative to the initial state. After 0.1 Myr of evolution, the simulations deviate by less a value of 0.05% from the steady state profile.
In the text

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.

In the text

	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).
In the text

	Fig. B.2 Gas radial velocity profiles (derived from the net mass flux, Eq. (22)) using our simple vertical structure approximation vs. Takeuchi & Lin (2002) model (Eqs. (B.1) and (B.2)).
In the text

	Fig. C.1 Surface density of gas (red) and dust (blue) after 0.4 Myr (solid lines) for different values of ε0 and αν. Initial conditions marked with dashed lines. The snow line location is marked with dotted lines.
In the text

	Fig. C.2 Top row: dust-to-gas ratio radial profile. Middle row: H2, He mass fraction radial profile. Bottom row: accretion rate time evolution (divided by the initial steady state accretion Ṁ0). All for different values of ε0 and αν. The value of 0.5 is marked with a dotted line in every plot.
In the text