Self-induced dust traps: overcoming planet formation barriers

Abstract

Planet formation is thought to occur in discs around young stars by the aggregation of small dust grains into much larger objects. The growth from grains to pebbles and from planetesimals to planets is now fairly well understood. The intermediate stage has however been found to be hindered by the radial-drift and fragmentation barriers. We identify a powerful mechanism in which dust overcomes both barriers. Its key ingredients are (i) backreaction from the dust on to the gas, (ii) grain growth and fragmentation and (iii) large-scale gradients. The pile-up of growing and fragmenting grains modifies the gas structure on large scales and triggers the formation of pressure maxima, in which particles are trapped. We show that these self-induced dust traps are robust: they develop for a wide range of disc structures, fragmentation thresholds and initial dust-to-gas ratios. They are favoured locations for pebbles to grow into planetesimals, thus opening new paths towards the formation of planets.

hydrodynamics, methods: numerical, protoplanetary discs

1 INTRODUCTION

In the core-accretion paradigm, planets are thought to originate from solid cores which accrete the surrounding material of their protoplanetary discs (Lissauer 1993; Pollack et al. 1996; Alibert et al. 2005; Mordasini et al. 2012). This requires primordial dust grains to concentrate and grow to form the building blocks of planets (Dominik et al. 2007). However, there are two important barriers that need to be overcome for this to happen. Grains feel the aerodynamic drag of the gas in the disc, causing the grains to settle to the mid-plane and drift inwards as they lose angular momentum. The efficiency of this relative motion of dust with respect to gas is controlled by the Stokes number St, i.e. the ratio of the drag stopping time to the orbital time, which depends on grain size and local gas conditions. Small (St ≪ 1, strong gas–dust coupling) and large (St ≫ 1, weak coupling) grains drift slowly, whereas grains with St ∼ 1 (i.e. millimetre sized at a few tens of au in observed discs, see e.g. Laibe, Gonzalez & Maddison 2012; Dipierro et al. 2015) experience the fastest radial drift, which can lead to their accretion on to the star before growing into planet embryos (Adachi, Hayashi & Nakazawa 1976; Weidenschilling 1977). Efficient grain growth is needed for dust particles to overcome this so-called radial-drift barrier (Laibe, Gonzalez & Maddison 2014). Larger decimetre-sized grains, on the other hand, which have higher relative velocities (Weidenschilling & Cuzzi 1993), shatter rather than stick upon collision, leading to the ‘fragmentation barrier’ (Blum & Wurm 2008). This raises the question: how, in spite of these barriers, do planets manage to form?

Dust grains drift towards the maximum of the gas pressure (Weidenschilling 1977), which for a smooth disc is located at the inner disc edge, resulting in their radial motion towards the central star. Any local gas pressure maximum in the disc will attract and trap grains both in the vicinity and those drifting inwards from the outer regions, forming a so-called dust trap (Whipple 1972). Such a trap suppresses any further radial drift of dust grains. The confinement of solids at this location increases their density, which accelerates grain growth (see Section 2), and reduces their relative velocities, which in turn prevents grain fragmentation. Dust traps are therefore an efficient way to overcome both the radial-drift and fragmentation barriers. Several types of dust traps have been proposed, e.g. vortices (Barge & Sommeria 1995; Lyra & Mac Low 2012; Méheut et al. 2012; Regály et al. 2012; Zhu et al. 2014) or planet gap edges (Paardekooper & Mellema 2004, 2006; Rice et al. 2006; Fouchet et al. 2007; Fouchet, Gonzalez & Maddison 2010; Ayliffe et al. 2012; Pinilla, Benisty & Birnstiel 2012; Zhu et al. 2012, 2014). They all require special conditions in the gas pressure profile or an external factor such as the presence of an existing planet.

Due to the high numerical cost, important processes are often neglected when modelling these dust traps. For example, there exist detailed hydrodynamical studies of instabilities in a gas (modelled as a fluid) and dust (represented by Lagrangian particles) disc with self-gravity but in local simulations and without grain growth or fragmentation (Johansen et al. 2007) or sophisticated treatments of growth and fragmentation but in 1D simulations of the dust evolution on a static gas background (Birnstiel, Dullemond & Brauer 2010). More recently, different models of dust evolution have been performed on a 2D static (Dra̧żkowska, Windmark & Dullemond 2013) or 1+1D evolving gas disc (Dra̧żkowska, Alibert & Moore 2016). To date, no numerical study has simultaneously treated all relevant physical ingredients in global 3D simulations of protoplanetary discs. Our two-fluid (gas+dust) smoothed particle hydrodynamics (SPH) code allows global 3D simulations of discs where the coupled evolution of gas and dust is treated self-consistently. We include both the aerodynamic drag of gas on dust and its backreaction, i.e. the drag of dust on gas, as well as the growth and fragmentation of dust grains (see Section 2). Backreaction is usually ignored because it is negligible when the dust-to-gas ratio ε is low, which is the case most of the time in astrophysical objects. Its effects become important, however, in situations where ε locally approaches unity, like in planet formation. Taking into account these processes simultaneously, in this work we demonstrate that dust traps can develop spontaneously in protoplanetary discs under a variety of conditions.

We detail our SPH code and growth/fragmentation model in Section 2 and present the resulting gas and dust distributions in Section 3. We discuss the origin of the self-induced dust traps, whose formation mechanism is summarized in Fig. 1, as well as explain how this mechanism differs from the streaming instability in Section 4. We summarize our work and conclude in Section 5.

Figure 1.

Self-induced dust trap formation mechanism. (1) Without backreaction, growing grains drift faster and faster as they approach St = 1, and either re-fragment in the inner disc or are accreted on to the star (grey). Backreaction slows down the drift, leaving enough time for grains to grow to large sizes, decouple from the gas and overcome radial drift (brown). (2) Once dust has piled up in a radial concentration, backreaction drags the gas outwards (blue), forming a gas pressure maximum, i.e. a dust trap. (3) The self-induced dust trap stops inward-drifting grains and strengthens the dust concentration.

2 METHODS

2.1 Hydrodynamical simulations

We use our 3D, two-phase (gas+dust), SPH code to model protoplanetary discs (Barrière-Fouchet et al. 2005). We adopt a vertically isothermal equation of state which mimics the internal thermal structure of typical discs below a thin surface layer (Woitke et al. 2016). The code does not include self-gravity. This has a negligible effect on our results since the disc-to-star mass ratio is of the order of 1 per cent.¹ Our implementation of gas–dust coupling via aerodynamic drag is that of Monaghan (1997) and conserves the total linear and angular momentum rigorously (the reader is referred to Monaghan 1997 for general tests and to Maddison 1998 for additional tests, specific to discs). In particular, we include here the drag of dust on gas, which we call the ‘backreaction’, often neglected in the literature. This term becomes important when properly describing the evolution of the disc mid-plane, where grains have settled and concentrated to reach dust-to-gas ratios of order unity. In that respect, we obtain larger dust densities in the mid-plane by treating dust settling dynamically, compared to 2D codes where equilibrium models are used. We adopt the standard SPH artificial viscosity (Monaghan 1989) with |$\alpha _\mathrm{\scriptscriptstyle AV}=0.1$| and |$\beta _\mathrm{\scriptscriptstyle AV}=0.5$|⁠, corresponding to a uniform Shakura–Sunyaev (Shakura & Sunyaev 1973) α ∼ 10⁻², as discussed in Fouchet et al. (2007). The scheme has been shown to naturally reproduce the expected properties of Prandtl-like turbulence (Arena & Gonzalez 2013). In particular, the dust scaleheight in the stationary regime for grains of constant size reproduces that predicted by Dubrulle, Morfill & Sterzik (1995).

We simulate discs of mass M_disc = 0.01 M_⊙ orbiting a 1 M_⊙ star. Our initial gas disc has a power-law surface density Σ_g = Σ₀(r/r₀)^−p and temperature T = T₀(r/r₀)^−q. The disc extends from r_in to r_out and is allowed to spread. Particles are removed from the simulation if they escape to radii larger than r_esc or if they drift inside of r_in, where they are assumed to be subsequently accreted by the star. The sound speed varies as c_s ∝ r^−q/2 and the disc aspect ratio as |$H/r\propto r^\frac{1-q}{2}$|⁠. We consider two disc models: a flat disc with p = 0 and q = 1, commonly used in planet-gap opening studies (Paardekooper & Mellema 2004, 2006; Fouchet et al. 2007, 2010; Gonzalez et al. 2015a), and a steep disc with p = 1 and q = 1/2, representing an average observed disc from Williams & Best (2014). The disc radial extension sets the value of Σ₀. The flat disc has a constant H/r = 0.05, which sets the temperature scale, whereas the steep disc has T₀ = 200 K at r₀ = 1 au. The relevant model parameters are listed in Table 1.

Table 1.

Parameters of our flat and steep disc models. Values with subscript 0 are given for r₀ = 1 au.

Model	p	q	rin	rout	resc	Σ0	T0	H0/R0
			(au)	(au)	(au)	(kg m−2)	(K)
Flat	0	1	4	120	160	19.67	623	0.05
Steep	1	1/2	10	300	400	487.74	200	0.0283

We start our simulations with 200 000 gas particles distributed according to the surface density profile with initial Keplerian velocities, and let the gas disc relax for eight orbits at r = 40 au for the flat disc and 24 orbits at r = 100 au for the steep disc. We then overlay the gas particles with an equal number of dust particles, with the same positions and velocities, to reproduce an initially uniform dust-to-gas ratio ε₀ and let the disc evolve. The gas particles have equal masses, as do the dust particles, but the two phases differ in mass by a factor ε₀. Solids have a typical bulk density ρ_s = 1000 kg m⁻³, an initial size s₀ = 10 μm and are allowed to grow and fragment. The results are insensitive to the value of s₀, since small grains initially grow very efficiently (Laibe et al. 2008). Fragmentation of dust grains occurs when dust relative velocities V_rel reach a threshold V_frag which is independent of the grain size. We explore values of V_frag of 10, 15, 20 and 25 m s⁻¹. These values of the fragmentation threshold span the expected range for materials from compact icy grains (10 m s⁻¹, Blum & Wurm 2008) to porous aggregates (several tens of m s⁻¹, Yamamoto, Kadono & Wada 2014). Since our modelled discs are entirely outside the water ice line, we do not take lower values, typical of compact silicate grains (1 m s⁻¹, Blum & Wurm 2008), into account. We consider initial dust-to-gas ratios ε₀ of 1 per cent, typical of the interstellar medium and commonly used for protoplanetary discs, for both the flat and steep discs, and higher values of 3 and 5 per cent in the steep disc model, as indicated by recent observations of a sample of Class II discs (Williams & Best 2014). The parameters used in the different models are provided in Table 2. Simulations are performed up to 200 000 yr for the flat disc and 400 000 yr for the steep disc, to allow the disc to form a stationary dust trap.

Table 2.

Simulation suite.

Model	ε0	Vfrag (m s−1)	Backreaction
Flat	1 per cent	10, 15, 20, 25	Yes
	1 per cent	15	No
Steep	1 per cent	10, 15, 20, 25	Yes
	1 per cent	15	No
	3 per cent	10	Yes
	5 per cent	10, 15, 20, 25	Yes

2.2 Drag in the Epstein regime

The gas densities of the discs we model are sufficiently low that the mean free path of gas molecules is larger than the grain sizes (Laibe et al. 2012). Hence, grains are in the Epstein aerodynamic drag regime and the expression of the Stokes number is

\begin{equation} \mathrm{St}\equiv \frac{\Omega _\mathrm{K}\rho _\mathrm{s}s}{\rho _\mathrm{g}c_\mathrm{s}}, \end{equation}

(1)

where Ω_K and ρ_g denote the Keplerian angular velocity and the gas density, respectively. Dust grains experiencing the fastest radial drift have St = 1, corresponding to a size

\begin{equation} s_\mathrm{drift} \equiv \frac{\rho _\mathrm{g}c_\mathrm{s}}{\rho _\mathrm{s}\Omega _\mathrm{K}}. \end{equation}

(2)

If grains cannot grow fast enough to overcome St = 1 before drifting out of the disc, the dust population is drift-limited: s_drift is the maximum size grains can reach at a given location in the disc.

2.3 Growth model

Grain growth is implemented as described in Laibe et al. (2008). The dust size distribution is assumed to be locally monodisperse, i.e. narrowly peaked around a local mean value, and two colliding grains with relative velocities below the fragmentation threshold will stick. The mass m of the resulting larger grain doubles in a mean collision time τ_coll, i.e. dm/dt = m/τ_coll, which translates to

\begin{equation} \frac{\mathrm{d}s}{\mathrm{d}t}=\frac{\rho _\mathrm{d}}{\rho _\mathrm{s}}\,V_\mathrm{rel}=\epsilon \,\frac{\rho _\mathrm{g}}{\rho _\mathrm{s}}\,V_\mathrm{rel}, \end{equation}

(3)

where ρ_d is the volume density of the dust phase and ε ≡ ρ_d/ρ_g is the dust-to-gas ratio (Stepinski & Valageas 1997). The relative velocity between grains V_rel originates from the gas turbulent motion, which is transferred to dust via aerodynamic drag. We adopt the subgrid model of turbulence developed by Stepinski & Valageas (1997), which gives

\begin{equation} V_\mathrm{rel}=\sqrt{2}\,V_\mathrm{t}\frac{\sqrt{\mathrm{Sc}-1}}{\mathrm{Sc}}. \end{equation}

(4)

|$V_\mathrm{t}\equiv \sqrt{\tilde{\alpha }}c_{\rm s}$| denotes the turbulent velocity, where |$\tilde{\alpha }\equiv 2^{1/2}\,\mathrm{Ro}\,\alpha$| is an effective turbulent diffusivity which depends on Ro, the Rossby number for turbulent motions (taken to be uniform and equal to 3). Sc denotes the Schmidt number of the grains, given by

\begin{equation} \mathrm{Sc} \equiv (1+\mathrm{St})\sqrt{1+\frac{\Delta \boldsymbol {v}^2}{V_\mathrm{t}^2}}, \end{equation}

(5)

where |$\Delta \boldsymbol {v}=\boldsymbol {v}_\mathrm{d}-\boldsymbol {v}_\mathrm{g}$| is the differential velocity between dust and gas. The time-scale for growth typically scales as (Ω_K ε)⁻¹ ∝ r^3/2 ε⁻¹ (Laibe et al. 2008): dust grains grow more rapidly in the inner disc and when the dust-to-gas ratio is large.

Our simulations show that |$\left| \Delta \boldsymbol {v} \right| \ll V_\mathrm{t}$|⁠, i.e. Sc is dominated by the effect of gas–dust coupling determined by the value of St. This gives a useful approximation for the variations of the relative velocities:

\begin{equation} V_\mathrm{rel}\simeq \sqrt{2\tilde{\alpha }}\,\frac{\sqrt{\mathrm{St}}}{1+\mathrm{St}}\,c_\mathrm{s}. \end{equation}

(6)

V_rel remains smaller than

\begin{equation} V_\mathrm{rel}^\mathrm{max}\equiv \sqrt{\frac{\tilde{\alpha }}{2}}\,c_\mathrm{s}, \end{equation}

(7)

which occurs when St = 1. In our simulations, |$V_\mathrm{rel}^\mathrm{max}\simeq 0.15\,c_\mathrm{s}$|⁠.

Our dust SPH particles have equal fixed masses and represent dust grains of the same size: with growth, they represent fewer physical grains of larger size (and mass), but still numerous enough to maintain the validity of the numerical scheme (Laibe et al. 2008). The assumption that growth occurs only between grains of the same size has little effect on the global size distribution. Indeed, a given volume contains SPH particles representing grains of different sizes. Our growth model results in a spread of (larger) sizes in that volume. If allowed, growth between grains of different sizes would produce sizes within that same range and only slightly depopulate both ends of the size distribution in that volume.

The spread in the relative velocities between gas and dust is treated in a Lagrangian way by our SPH code. This results in a spread in V_rel at a given location in the disc, similarly to the velocity distributions used by Windmark et al. (2012) and Garaud et al. (2013). The subgrid description of turbulence is a source of uncertainty of the model. Several expressions for V_rel have been proposed in the literature, e.g. Ormel & Cuzzi (2007) and Pan & Padoan (2010) – see Laibe (2014) for a discussion. We have run test simulations with several of them, producing qualitatively similar dust behaviours and, importantly, they all lead to a dust pile-up.

2.4 Fragmentation model

When two particles collide with high relative kinetic energy, they break into a collection of fragments of various sizes. In a steady state, the size distribution of fragments is commonly represented by a power law. We implement fragmentation by comparing the relative velocity between grains V_rel to the fragmentation threshold V_frag. If V_rel > V_frag, grains shatter, which we model by decreasing the size of the representative SPH particles following

\begin{equation} \frac{\mathrm{d}s}{\mathrm{d}t}=-\epsilon \,\frac{\rho _\mathrm{g}}{\rho _\mathrm{s}}\,V_\mathrm{rel}. \end{equation}

(8)

For grains of a given Stokes number, fragmentation occurs preferentially in the inner disc regions, where they have higher relative velocities (see equation 6, in which c_s decreases as a function of r). Appendix A gives the initial location of the limit between growth and fragmentation as a function of grain size and distance to the star in our disc models. To ensure consistency with the growth model, as well as to properly conserve physical quantities and avoid a prohibitively large increase in the number of SPH particles, we keep a locally monodisperse size distribution. This means we follow only the largest fragment of the distribution. Equation (8) is equivalent to dm/dt = −m/τ_coll and implies that the largest fragment has lost most of the original grain's mass, i.e. that fragmentation leads to small sizes. Our model is thus conservative: by producing fragments that are much smaller than their parent grains, it delays their subsequent evolution to larger sizes and makes planetesimal formation more difficult. For numerical efficiency, we limit the minimum fragmentation size to s₀ = 10 μm. This choice has little effect on the overall size distribution. Indeed, very small grains always have small V_rel – see equation (6) – hence are unaffected by fragmentation. They grow very efficiently and quickly forget their initial size (Laibe et al. 2008).

For a range of velocities below V_frag, colliding grains are expected to bounce instead of sticking or shattering (Zsom et al. 2010), which we do not take into account. This would lower the global growth efficiency, and can be mimicked to first order by reducing V_frag (however, this would still produce small fragments, which are not expected in bouncing collisions). We chose to keep only one free parameter and not to include bouncing in this work, leaving an exploration of its effects for a future study. Our treatment is simplified compared to that of Birnstiel et al. (2010) or Garaud et al. (2013), in particular we do not keep track of the smallest fragments who maintain a population of small grains in the disc that are responsible for the infrared emission, nor do we consider collisions between unequal-sized particles. Yet, the resulting global size distribution we achieve when backreaction is turned off is very similar to that obtained by Brauer, Dullemond & Henning (2008) and Birnstiel et al. (2010), where backreaction was neglected (see Section 4).

3 RESULTS

The evolution of the radial grain size distribution in the flat model with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent is displayed in Fig. 2. Its full evolution is shown in Video 1 (available online). At early times, growing grains fragment almost immediately in most of the disc, keeping sizes smaller than 100 μm. Only grains in the very outer disc have relative velocities low enough to avoid fragmentation and grow slowly. After ∼80 000 yr, they have reached the size where St = 1 (∼1 mm for these disc conditions), drift inwards and continue to grow. At ∼100 000 yr, the largest grains begin decoupling from the gas, slow their radial drift and pile up between ∼70 and 100 au. Locally, ε becomes large, which causes the backreaction to strongly reduce the dust drift velocity, and that is the case for all values of St, including the fastest drift regime when St ∼ 1 (Nakagawa, Sekiya & Hayashi 1986). This amplifies this pile-up, helping grains escape the radial-drift barrier. To see a similar behaviour without backreaction, Birnstiel et al. (2010) had to artificially decrease the drift efficiency by 80 per cent. The critical point is that backreaction also adds an outwards term to the gas radial velocity, usually negligible compared to the inwards accretion flow, except again when ε is large and St ∼ 1 (see Section 4.3 and Appendix C for more details). Near the radial dust concentration, gas is locally dragged outwards, accumulates and creates a pressure maximum. This results in a self-induced dust trap, with no need for any special condition. After 200 000 yr, most of the dust grains are concentrated in a dense ring extending from ∼40 to ∼65 au, have grown to cm sizes and are well decoupled from the gas, with low relative velocities. For higher values of V_frag, the dust trap forms closer to the star (see Appendix B for details).

Figure 2.

Radial grain size distribution in the flat disc for V_frag = 15 m s⁻¹ and ε₀ = 1 per cent. Grains grow, decouple from the gas and pile-up, resulting in a self-induced dust trap. The blue/red colour code represents the Stokes number St (left), and the brown/green colour code represents the ratio V_rel/V_frag (right). Four snapshots at 50 000, 100 000, 150 000 and 200 000 yr are shown, from top to bottom. Video 1 (available online) shows the full evolution.

We performed the same simulation without the backreaction in order to emphasize the critical role it plays. Fig. 3 shows that no grain has been able to escape the radial-drift barrier (see Video 1, available online, for the full evolution), similarly to what is found in other studies in which backreaction is neglected, e.g. Brauer et al. (2008), Birnstiel et al. (2010) and Windmark et al. (2012). Once grains have reached the size for which St = 1, they drift inwards to regions where their relative velocities are high and exceed V_frag. The grains then fragment and replenish the small grains reservoir. The largest fragments keep drifting and are accreted on to the star. No pile-up forms in this case.

Figure 3.

Top: same as Fig. 2 after 200 000 yr. Bottom: same simulation with backreaction turned off. In this case, when grains grow and reach St ∼1, they drift rapidly towards the central star while fragmenting and replenishing the small grains reservoir. Video 1 (available online) shows the full evolution.

Self-induced dust traps also develop in the steep disc for a range of conditions. For ε₀ = 1 per cent and V_frag = 10 m s⁻¹, a trap forms at ∼30 au (top panel of Fig. 4). Behind it, the region of large, decoupled grains extends out to ∼100 au because of the different disc structure. For higher values of V_frag, fragmentation is less of a hindrance and grains can grow efficiently at smaller radial distances where the growth time-scale, varying as r^3/2 (Laibe et al. 2008), is shortest. This behaviour is similar to the case where fragmentation is not included (equivalent to V_frag = +∞) and grains pile up very efficiently in the inner disc, as predicted theoretically (Laibe 2014) and seen in simulations (Gonzalez et al. 2015a), at the same location for all values of V_frag ≥ 15 m s⁻¹ (bottom three panels of Fig. 4). The gas then accumulates and the self-induced dust trap forms. Note that this radial concentration is not caused by the inner boundary: grains initially in the inner disc drift inwards and leave the disc before dust piles up at 15–20 au (see Video 2, available online, which shows the evolution of radial grain size distribution in the steep disc with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent, with and without backreaction). For higher values of ε₀, grains grow faster (equation 3), decouple sooner during their inward drift, and therefore pile up at larger radii – see Fig. 5. Results for the full range of values of ε₀ and V_frag are detailed in Appendix B.

Figure 4.

Radial grain size distribution in the steep disc after 400 000 yr, for ε₀ = 1 per cent and different values of the fragmentation velocity. From top to bottom: V_frag = 10, 15, 20 and 25 m s⁻¹. For V_frag = 10 m s⁻¹, a self-induced dust trap forms at r ∼ 30 au, at the fragmentation front of the particles where V_rel(St = 1) = V_frag. For V_frag ≥ 15 m s⁻¹, the trap forms at r ∼ 15–20 au, where growth induces a strong pile-up. A log scale has been used to emphasize these different locations.

Figure 5.

Radial grain size distribution in the steep disc after 400 000 yr, for V_frag = 10 m s⁻¹ and different values of the initial dust-to-gas ratio. From top to bottom: ε₀ = 1, 3 and 5 per cent. Larger grain concentrations give rise to more effective dust trapping.

4 DISCUSSION

The formation mechanism of the self-induced dust trap is summarized in the sketch presented in Fig. 1, and detailed in the following subsections.

4.1 Without backreaction: radial-drift and fragmentation barriers

The first simulations of the dust evolution including fragmentation were performed by Brauer et al. (2008) and Birnstiel et al. (2010) in 1D along the radial direction on a fixed, vertically and azimuthally averaged gas disc structure, and neglecting backreaction. Their results are qualitatively reproduced in our 3D simulations with backreaction turned off, in which V_frag = 15 m s⁻¹ and ε₀ = 1 per cent (Fig. 3 and bottom panels of Video 1 for the flat disc and of Video 2 for the steep disc). In the inner disc, growing grains re-fragment rapidly and repopulate the reservoir of small grains, which drift very slowly. Grains growing in the outer disc reach sizes close to s_drift and are dragged towards the interior where they fragment to small sizes or are lost from the disc (see the bottom-right panels of Videos 1 and 2). Grains are not able to grow above s_drift, and hence they do not break the radial-drift barrier. This was the main result of previous studies neglecting backreaction (Brauer et al. 2008; Birnstiel et al. 2010). The dust-to-gas ratio shows an inward moving front, illustrating the dust drift (shown for the flat disc in the bottom panel of Fig. 6). In the mid-plane, ε is larger than its vertical average, due to dust settling, and reaches about 10 per cent. This 10-fold increase is consistent with the thickness of the dust layer predicted by Dubrulle et al. (1995) for α = 10⁻². Fig. 7 shows the final radial grain size distribution in the flat disc, colour-coded by altitude in the disc. It shows that the bulk of the grain population is close to the mid-plane and is drift-limited (the grain sizes stay below s_drift). In addition, Fig. 7 shows that the grains closer to the disc surface have sizes below the fragmentation size |$s_\mathrm{frag}^-$| (see Appendix A for its derivation), i.e. that they are fragmentation limited, something 1D simulations cannot see. Such simulations would also only notice a moderate increase in the vertically averaged dust-to-gas ratio (Fig. 6).

Figure 6.

Radial profiles of the dust-to-gas ratio ε in the flat disc with an initial value ε₀ = 1 per cent (marked by the horizontal dotted line), for V_frag = 15 m s⁻¹. Its mid-plane and vertically integrated values are shown at 100 000, 150 000 and 200 000 yr for the simulations with (top) and without (bottom) backreaction. The dust-to-gas ratio is enhanced in the mid-plane because of vertical dust settling. All radial profiles shown in this figure as well as in Figs 7, 8, 9, 10, B1, B2 and B3 are azimuthally averaged.

Figure 7.

Radial grain size distribution in the flat disc with V_frag = 15 m s⁻¹, ε₀ = 1 per cent without backreaction after 200 000 yr, colour-coded by altitude z. Radial profiles of the sizes s_drift (dashed line) and |$s_\mathrm{frag}^\pm$| (dash–dotted lines) are overplotted.

4.2 Backreaction effect no. 1: radial concentration of grains

When backreaction is considered, the dust drift velocity is reduced by a factor (1 + ε)², for all values of St ≲ 1 (Nakagawa et al. 1986). Negligible at first, this effect becomes important when dust settles and ε increases in the mid-plane. This results in a radically different dust evolution compared to simulations without backreaction (Fig. 2 and top panels of Video 1 for the flat disc and of Video 2 for the steep disc, for V_frag = 15 m s⁻¹). Initially, fragmentation dominates the dust evolution in the inner disc, similarly to the case without backreaction. However, grains exterior to r_frag grow unhindered by fragmentation and backreaction acts to slow down their drift. This helps them to grow above |$s_\mathrm{frag}^+$| (the top-right panels of Videos 1 and 2 show that they have V_rel < V_frag) before entering the inner fragmentation region (see Appendix A). They progressively decouple from the gas and slow down. In the flat disc, after ∼100 000 yr, dust has piled up and formed a radial concentration of marginally coupled dust grains at ∼90 au, as can be seen in the radial profile of ε in the mid-plane in Fig. 6. (Note that the large values of ε interior to 20–30 au are due to a decrease of the gas density rather than an increase of the dust density.) The dust concentration drifts slowly and keeps decelerating, helped by the increasing importance of backreaction when ε approaches unity. The drift stalls, with the inner edge of the dust concentration around 40 au.

4.3 Backreaction effect no. 2: dust concentrations modify the gas structure

The first stage of the self-induced dust trapping mechanism consists in the formation of the strongly peaked radial concentration of dust grains with St ∼ 1, as just described. We now explain the second stage, i.e. how dust backreaction modifies the gas structure around the dust concentration. In Appendix C, we give the analytical expression of the gas radial velocity as the sum of the viscous inwards flow and the outwards motion due to dust drag. In most situations, the viscous flow dominates, but in dust concentrations ε becomes large and the outwards drag locally becomes the dominant term. To demonstrate the effect on the gas structure, we numerically integrate the evolution equation of the gas surface density in the presence of a sharp dust concentration. We observe an outwards mass transfer and a gas accumulation at the outer edge of the dust concentration (see Fig. C2), in agreement with the analytical prediction. The gas evolution equation (equation C11) stresses the role played by the following three key ingredients: (i) backreaction, whose effect on the gas phase is represented by the f_drag term (defined in Appendix C), (ii) growth and fragmentation, assisted by backreaction, which produce local dust concentrations with ε ≃ St ≃ 1, hence large values of f_drag and (iii) global dust dynamics, giving rise to relevant large-scale gradients of f_drag. This mechanism cannot be simulated using local approximations.

Quantitatively, the relative importance between backreaction and viscosity is measured by the parameter x_br, a function of ε and St defined in equation (C8). As an example, we discuss the role played by backreaction in the case of the flat disc with V_frag = 15 m s⁻¹ in Fig. 8. Before 100 000 yr, backreaction is negligible (x_br ≪ 1 since ε and St are small). However, around 100 000 yr, the motion of the gas is driven by backreaction since both ε and St approach unity and x_br ≫ 1. The resulting outwards flow starts to modify the gas pressure profile, causing a slight bump at ∼90 au shown in Fig. 9. At later stages, the disc region affected by backreaction has drifted inwards, together with the radial concentration of slowly decoupling grains (Fig. 8, x_br decreases but remains well above unity). Fig. 9 shows the formation of a pressure maximum, which has become stationary by 200 000 yr, when the grains in the radial concentration have fully decoupled (see Fig. 2), and centred around 40 au. At this stage, the viscous and drag terms balance each other (Fig. 8 shows that x_br ∼ 1 at the location of the pressure maximum).

Figure 8.

Evolution of ε, St and x_br in the mid-plane of the flat disc with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent.

Figure 9.

Evolution of the gas pressure (top) and pressure gradient (bottom) in the mid-plane of the flat disc with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent.

The motion of the gas phase is illustrated by the profile of its radial density flux, plotted in Fig. 10. At 100 000 yr, it shows outwards mass transport exterior to ∼30 au and vanishing at ∼90 au where gas thus accumulates in a nascent pressure bump (Fig. 9). Outwards motion is seen again further out as the disc viscously spreads. At 125 000 and 150 000 yr, the flux profiles exhibit adjacent regions of outwards and inwards mass transport, i.e. converging motion towards 65 and 40 au, respectively, strengthening the pressure maximum seen in Fig. 9. At the end of the simulation, no gas motion is seen when the pressure maximum has fully formed and is stationary. As expected, no such behaviour is observed in the simulation without backreaction, which shows only the usual accretion in the inner disc and viscous spreading in the outer disc.

Figure 10.

Evolution of the radial density flux of the gas in the mid-plane of the flat disc with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent for the simulations with and without backreaction.

A dust trap has thus formed at 40 au, coincident with the inner edge of the dust accumulation seen in Figs 2 and 6 at the end of the simulation. Note that even though the gas accumulates at the outer edge of the initial (and narrow) dust concentration (Appendix C), inward drifting grains are subsequently stopped by the gas pressure maximum and thus pile up exterior to it, forming a wider dust accumulation.

4.4 Self-induced dust traps

The combined evolution of the gas and dust phases detailed above leads to the formation of a gas pressure maximum, which acts as a dust trap. This trap does not result from a pre-existing discontinuity in the disc, such as a dead zone edge or a snow line or from the opening a planet gap, but forms spontaneously – a self-induced dust trap. It then acts as any other dust trap and stops the inward motion of grains drifting from the outer disc, which accumulates them further. Trapped grains have sufficiently low collision velocities to stick and continue growing. At the end of our simulations, the mass of trapped dust ranges from a few to a few tens of Earth masses (with higher values for larger V_frag and larger ε₀). Enough solid mass has accumulated to allow the formation of planetesimals. Modelling their actual formation would require self-gravity, which is not included in our simulations, and is out of the scope of this paper.

Backreaction plays a notable role at every step of the evolution. First, it enhances the radial concentration of dust grains by slowing down their drift, which increases the dust-to-gas ratio ε and accelerates grain growth, which in turn slows the particles even more. Secondly, it reshapes the gas structure by transferring gas from the inner to the outer regions. A pressure maximum forms, preventing further drift of the particles towards the star.

In summary, without backreaction the dust population stays drift- or fragmentation-limited, as commonly seen by other authors. It is only when the self-consistent evolution of both gas and dust is fully taken into account, with backreaction included, that a different behaviour emerges: drifting, growing and fragmenting grains are able to break the radial-drift and fragmentation barriers and form a stationary radial concentration.

This self-induced trapping mechanism is robust. Whatever the disc conditions, fragmentation threshold, or initial dust-to-gas ratio, dust traps form, but at different locations and different evolutionary times. The larger the dust mass, the stronger the effect. In particular, higher dust-to-gas ratios cause grains to pile up at larger radii, as well as grow to larger sizes. If they end up forming planetesimals, this could explain the presence of planets orbiting several tens of au from their stars: a value of ε₀ slightly larger than average is enough for pebbles to accumulate in the outer disc.

Note that, since grains in our self-induced dust traps have reached mm to cm sizes, they would appear as bright rings in dust continuum observations in the mm range (see Gonzalez et al. 2015b for ALMA – the Atacama Large Millimeter/submillimeter Array – simulated images of such a dust concentration). At later times, if planetesimals form at the trap location, a deficit of emission would be seen as a concentric gap.

4.5 Comparison with the streaming instability

We aim to compare the formation of the self-induced dust traps identified above (defined by the equations given in Appendix C) to the streaming instability, discovered and defined by the set of equations given in Youdin & Goodman (2005). Both mechanisms are particular cases of two-stream instabilities in discs, which also encompass the generation of non-axisymmetric structures by dust backreaction (Lyra & Kuchner 2013). Both of these instabilities require a dust-and-gas mixture in rotation, and an additional background pressure gradient in the gas.

The streaming instability develops in the (r, z) plane of a disc whose background pressure gradient is enforced to remain constant. The transient linear growth stage has been studied analytically by Youdin & Goodman (2005), Youdin & Johansen (2007) and Jacquet, Balbus & Latter (2011) and shows efficient local enhancement of the dust density. The non-linear regime has been investigated numerically (Johansen et al. 2007; Bai & Stone 2010; Yang & Johansen 2014) and leads to strong dust clumping, assisted by the slowdown of radial drift in overdense regions. The streaming instability is a proposed mechanism to overcome the barriers of planet formation and leads to the formation of planetesimals once self-gravity becomes important (Johansen et al. 2007). The self-induced dust traps we describe in this work also lead to strong concentrations of dust grains, which no longer drift and can grow unhindered by fragmentation. They provide therefore a solution to the barriers of planet formation as well. The two mechanisms share similar physical ingredients, but we stress their differences below since they have important consequences for planet formation.

The streaming instability can grow even in the limit of a rigorously incompressible gas (Youdin & Goodman 2005), i.e. gas compression is not required for the local dust-to-gas ratio to be amplified. The perturbation must develop in the radial but also in the vertical direction to grow (i.e. the flow is not unstable if the wavenumber k_z = 0, e.g. see equations 29 and 32 of Jacquet et al. 2011 in the terminal velocity approximation). Self-induced dust traps require large-scale background gradients for the dust-to-gas ratio and the Stokes number, assisted by an efficient compression of the gas in the mid-plane to develop, i.e. the formation of the pressure maximum is triggered by the term on the right-hand side of equation (C11), which comes from the compressibility term |$\Sigma _\mathrm{g} {\nabla \cdot \boldsymbol {v}} \ne 0$|⁠. The instability develops along the sole radial direction. Global simulations of evolved dust distributions with backreaction are therefore needed to observe the formation of self-induced dust traps.

The time-scale for the streaming instability to develop under the terminal velocity approximation, when St ≪ 1 is |$\mathcal {O}\left((r/H)^{2}t_\mathrm{K}\right)$|⁠, see Youdin & Goodman (2005). In simulations, it is found to be ∼10–40 t_K (Johansen & Youdin 2007; Bai & Stone 2010; Yang & Johansen 2014). This is somewhat longer than the typical time required for the self-induced dust trap to form, τ ∼ 5 t_K (see equation C12). This suggests that where large-scale radial dust concentrations build up, self-induced dust traps should form, and that small-scale streaming instability may develop over longer time-scales. Concurrent appearance of both cannot be ruled out. The streaming instability and self-induced dust traps both lead to the formation of dust clumps, but along different paths. The streaming instability starts to concentrate solids as early as in its transient linear growth stage. Self-induced dust traps require first a local enhancement of the dust-to-gas ratio by settling or radial pile-up. Then, a pressure maximum forms, without affecting the local dust density. Dust is then collected by this trap over longer time-scales, as particles drifting from the outer disc regions accumulate in the pressure maximum.

Finally, the streaming instability requires low viscosities for the perturbation to grow without being damped by viscous dissipation. Equation 33 of Youdin & Goodman (2005) gives a maximum value of

\begin{equation} \alpha \lesssim \frac{g}{\Omega _\mathrm{K}} \left(\frac{H}{r} \frac{2\pi }{K}\right)^2, \end{equation}

(9)

where g denotes the linear growth rate of the instability and K the normalized wavenumber. In the favourable case where St ∼ 1 and ε ∼ 0.1–1, Youdin & Johansen (2007) find g/Ω_K ∼ 0.1 and K ∼ 1, which gives, with H/r ∼ 0.05, a maximum value α ∼ 10⁻², comparable to the values in our discs. These conditions are however found only transitorily in our simulations and the streaming instability would not be operative in our large viscosity discs when St and ε have different values. For example, for St ∼ 0.1 and ε ∼ 0.1–1, the streaming instability develops only when α ≲ 2–8 × 10⁻⁵. Self-induced dust traps, on the other hand, have no problem forming with large viscosities: backreaction dominates the gas motion for a wide range of St and ε values when α = 10⁻² – see Fig. C1 – and for a reduced range for α as high as 10⁻¹ (equation C8). We note that the simulations presented here lack the necessary resolution to pick up the small-scale fastest growing mode of the streaming instability. This is of no consequence in this work, since the viscosity is large enough for the streaming instability to operate only marginally.

Thus, although the streaming instability and self-induced dust traps involve some common physical processes, they are two distinct mechanisms with different outcomes (see also Section 4.6). We recommend to refer to these two mechanisms by two distinct names.

4.6 Where do planetesimals form?

The expected location and number of planetesimals are tightly related to the underlying formation mechanism. Self-induced dust traps favour an effective concentration of solids, but this concerns only a small number of rings located at a few tens of au from the central star, at most. If planetesimals form inside the pressure bump by gravity in a few orbits, a significant population of observable grains should remain preserved outside of the trap. If planetesimals form in low-viscosity discs by the streaming instability, there is no physical reason to expect the formation to be local: the conditions for the instability to develop should be fulfilled in extended regions of the disc, if not everywhere. Moreover, numerical simulations show that once activated, it is a powerful mechanism to convert grains into planetesimals (Johansen et al. 2007; Bai & Stone 2010). Hence, planetesimal formation by the streaming instability should be global and total, preserving only a tiny population of grains in regions where it develops. Thus, self-induced dust traps and the streaming instability provide two potential paths towards planetesimal formation. Statistics over future spatially resolved observations of young discs might help to determine the relevance of each of them.

5 CONCLUSION

The core-accretion paradigm of planet formation requires the concentration and growth of primordial grains, a mechanism assumed to be hampered by the radial-drift and fragmentation barriers. In this work, we show that both of these barriers to planet formation can be overcome by a self-induced dust trap, which requires no special conditions or new physics. The key ingredients are the self-consistent inclusion of the backreaction of the dust on to the gas due to aerodynamic drag, which becomes important when the dust-to-gas ratio locally approaches unity, along with growth and fragmentation of dust grains. In global disc simulations of dust dynamics, growing grains that have settled to the disc mid-plane progressively decouple from the gas and slow their radial migration. The dust piles up and these radial dust concentrations modify the gas structure on large scales due to the backreaction, producing a gas pressure maximum that acts as a dust trap. The low collision velocities in the dust trap further enhances grain growth, rapidly leading to pebble-sized grains. The subsequent growth of these pebbles to planetesimals can progress via either collisions or gravitational collapse as demonstrated by Johansen et al. (2014) and Levison, Kretke & Duncan (2015).

We demonstrate that this process is extremely robust and that self-induced dust traps form in different disc structures, with different fragmentation thresholds, and for a variety of initial dust-to-gas ratios. Changing these parameters result in self-induced dust traps at different locations in the disc, and at different evolutionary times. While seemingly counter-intuitive, fragmentation is a vital ingredient for planet formation as it helps to form dust traps at large distances from the star. Indeed, fragmentation only allows grains to grow exterior to a certain radial distance and when grains decouple from the gas and start piling up, they do so near that radius. Stronger fragmentation, with a lower fragmentation threshold, implies that this radius lies farther away from the star. This would suggest that most discs thus retain and concentrate their grains at specific locations in time-scales compatible with recent observations of structures in young stellar objects (ALMA Partnership et al. 2015; Andrews et al. 2016).

The critical role of the backreaction is demonstrated by comparing with simulations in which it is not included. As seen in previous work, without backreaction a gas pressure maximum does not form and hence grains continue drifting inwards towards the star. Grains with relative large velocities cannot overcome the fragmentation barrier, fragmenting upon collision and replenishing the small grain population, while the largest grains cannot overcome the radial-drift barrier and are rapidly accreted into the central star.

The self-induced dust trap presented here is not the same mechanism as the streaming instability. The differences between both have consequences for planetesimal formation: (i) self-induced dust traps develop on somewhat shorter time-scales than the streaming instability, (ii) they can develop even in highly viscous discs (α ≃ 10⁻²–10⁻¹) and (iii) they promote planetesimal formation in a small number of rings while the streaming instability favours large regions of the disc. It is likely, however, that both play a part in planet formation.

In this work, we have presented a powerful new mechanism to overcome the main bottleneck of planet formation – the growth from micrometre-sized grains to centimetre-sized and decimetre-sized pebbles. This vital phase of grain growth is mostly hindered by the radial-drift barrier and the fragmentation barrier. Once grains have reached pebble sizes, there are several possible pathways for their subsequent growth to planetesimals. Thus, self-induced dust traps bring a critical missing piece of the planet formation puzzle and provide a favourable environment for the growth towards planetesimals.

Acknowledgments

We thank the anonymous referee for thorough comments that helped us to clarify and improve the manuscript. This research was partially supported by the Programme National de Physique Stellaire and the Programme National de Planétologie of CNRS/INSU, France. JFG thanks the LABEX (LABoratoire d'EXcellence) Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the programme ‘Investissements d'Avenir’ (ANR-11-IDEX-0007) of the French government operated by the ANR (Agence Nationale de la Recherche). GL is grateful for funding from the European Research Council for the FP7 ERC advanced grant project ECOGAL. STM acknowledges partial support from PALSE (Programme Avenir Lyon Saint-Etienne). Simulations were run at the Common Computing Facility (CCF) of LABEX LIO. Figs 2, 3, 4, 5 and 7 were made with SPLASH (Price 2007).

1

Only our flat disc model (see Table 1) is marginally unstable to gravitational instability near its outer rim, which would help planet formation there. However, we are studying a different mechanism and neglecting self-gravity allows us to separate the effects.

2

We have run test simulations with several of them, producing qualitatively similar dust behaviours and, importantly, they all lead to a dust pile-up.

APPENDIX A: Growth versus fragmentation

At a given location in the disc, the grain size separating the ranges of growing and fragmenting grains is estimated by writing V_rel = V_frag in combination with equation (6). This gives a second order equation for St:

\begin{equation} \mathrm{St}^2+2\left(1-\frac{\tilde{\alpha }\,c_\mathrm{s}^2}{V_\mathrm{frag}^2}\right)\mathrm{St}+1=0. \end{equation}

(A1)

When |$V_\mathrm{frag}>V_\mathrm{rel}^\mathrm{max}$|⁠, equation (A1) has no solution: V_rel < V_frag for all sizes and grains always grow. From equation (7), this occurs exterior to a radius

\begin{equation} r_\mathrm{frag}=r_0\left(\frac{\tilde{\alpha }\,c_\mathrm{s0}^2}{2V_\mathrm{frag}^2}\right)^{{1}/{q}}, \end{equation}

(A2)

where the disc is colder and the turbulence less effective, resulting in low relative velocities between dust grains. Otherwise, the solutions to equation (A1) can be translated into limiting sizes for fragmentation

\begin{equation} s_\mathrm{frag}^\pm =s_\mathrm{drift}\left[\frac{\tilde{\alpha }\,c_\mathrm{s}^2}{V_\mathrm{frag}^2}-1\pm \frac{\sqrt{\tilde{\alpha }}\,c_\mathrm{s}}{V_\mathrm{frag}}\sqrt{\frac{\tilde{\alpha }\,c_\mathrm{s}^2}{V_\mathrm{frag}^2}-2}\ \right]. \end{equation}

(A3)

Grains grow if |$s<s_\mathrm{frag}^-$| or |$s>s_\mathrm{frag}^+$| and fragment if |$s_\mathrm{frag}^-<s<s_\mathrm{frag}^+$|⁠.

The radial profiles of both limiting sizes |$s_\mathrm{frag}^-$| and |$s_\mathrm{frag}^+$|⁠, and of the size of fastest drift s_drift, are plotted in Fig. A1 for the initial power-law setup of the flat and steep discs. At early times, grains are small and grow. In the inner regions of both discs, grains hit the fragmentation barrier at sizes |$s_\mathrm{frag}^-$| which are well below s_drift, thus keeping St ≪ 1 and stay small and well coupled to the gas. They drift very little and remain in the disc. Exterior to r_frag, grains grow unhampered by fragmentation and start drifting rapidly as their size approaches s_drift. If they drift faster than they grow, they will reach the fragmentation barrier, shatter to small sizes and stop drifting. Otherwise, they stay above |$s_\mathrm{frag}^+$|⁠, keep growing and progressively reach large St values, decouple from the gas and drift more and more slowly. Note that the profiles of s_drift and |$s_\mathrm{frag}^\pm$| evolve with time, following the evolution of the gas density (see equation 2). In particular, ρ_g decreases in the outer disc due to viscous spreading, which shifts the profiles to lower sizes and allows more grains to have |$s>s_\mathrm{frag}^+$|⁠.

Figure A1.

Initial radial profiles of the sizes s_drift and |$s_\mathrm{frag}^\pm$| in the power-law setup of the flat (top) and steep (bottom) discs. Between the branches |$s_\mathrm{frag}^-$| and |$s_\mathrm{frag}^+$|⁠, V_rel > V_frag and grains fragment, while they grow below |$s_\mathrm{frag}^-$| and above |$s_\mathrm{frag}^+$|⁠. Exterior to a radius r_frag (which depends on V_frag), grains always grow.

APPENDIX B: INITIAL DUST-TO-GAS RATIO AND FRAGMENTATION THRESHOLD

The dust growth rate ds/dt is proportional to the dust-to-gas ratio ε (equation 3) and the drift velocity dr/dt is a decreasing function of ε (Nakagawa et al. 1986). Therefore, |$\mathrm{d}s/\mathrm{d}r=(\mathrm{d}s/\mathrm{d}t)/(\mathrm{d}r/\mathrm{d}t)$| is an increasing function of ε. In discs with larger ε₀, grains grow faster, and when they decouple from the gas and stop drifting, they are farther away from the star (in the r–s plane, the slope of a grain's trajectory is steeper, similarly to what is seen in the top panel of Fig. 1 in the r–St plane, comparing the cases without and with backreaction). The self-induced dust traps form sooner and at larger radial distances, as can be seen in Fig. 5, presenting the radial grain size distribution for the steep disc with V_frag = 10 m s⁻¹ and ε₀ = 1, 3 and 5 per cent (Video 3, available online, shows its evolution) and in Fig. B1, plotting the corresponding gas pressure profiles. After the same evolutionary time, grains have grown to larger sizes.

Figure B1.

Radial profiles of the gas pressure in the mid-plane of the steep disc with V_frag = 10 m s⁻¹ after 400 000 yr for ε₀ = 1, 3 and 5 per cent.

The radius exterior to which fragmentation is not effective varies as |$r_\mathrm{frag}\propto V_\mathrm{frag}^{-2/q}$| (equation A2). For larger fragmentation thresholds, grains can grow unhindered by fragmentation closer to the star, and therefore more rapidly due to the radial dependence of the growth time-scale. Those grains are the ones that eventually overcome the radial-drift and fragmentation barriers and pile up. The self-induced dust trap thus forms closer to the star and sooner. This is illustrated in Fig. B2, showing the locations of the gas pressure maxima for the flat disc with ε₀ = 1 per cent and values of V_frag from 10 to 25 m s⁻¹. In the steep disc, the same trends are observed (see Fig. B3 and Video 4, available online) with some additional peculiarities. For ε₀ = 1 per cent, when V_frag ≥ 15 m s⁻¹, it is sufficiently large that grains manage to decouple very efficiently from the gas as they grow before reaching r_frag. Their subsequent evolution is thus independent of V_frag and is the same as for growing grains when fragmentation is not taken into account, i.e. V_frag = +∞ (Laibe et al. 2008; Laibe 2014). After a phase of rapid drift, they decouple from the gas and pile up in the inner disc. The self-induced dust trap thus forms at the same position for all values of V_frag ≥ 15 m s⁻¹ (Fig. B3, top). For ε₀ = 5 per cent, grain growth is so fast that a first dust trap forms very early and at large distances, between 100 and 200 au. The gas compression at the location of the trap lowers the gas surface density, and therefore |$s_\mathrm{frag}^+$|⁠, around it, which helps grains slightly interior to it to escape fragmentation. They drift inwards and grow, forming a second dust trap in the inner disc as discussed previously, at shorter radii for larger V_frag (Fig. B3, bottom). The formation of these two self-induced dust traps can only be seen in global simulations.

Figure B2.

Radial profiles of the gas pressure (top) and pressure gradient (bottom) in the mid-plane of the flat disc with ε₀ = 1 per cent after 200 000 yr for V_frag = 10, 15, 20 and 25 m s⁻¹.

Figure B3.

Radial profiles of the gas pressure in the mid-plane of the steep disc after 400 000 yr for V_frag = 10, 15, 20 and 25 m s⁻¹. Top: ε₀ = 1 per cent. Bottom: ε₀ = 5 per cent.

APPENDIX C: FORMING THE TRAP

C1 Equations of motion

After a few stopping times, the radial velocity of the viscous gas phase of a dusty disc is given by

\begin{equation} \begin{array}{r@{\, }c@{\, }l}v_\mathrm{g} & = & -f_\mathrm{drag} \, \displaystyle \frac{1}{\Sigma _\mathrm{g} \Omega } \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(c_\mathrm{s}^2 \Sigma _\mathrm{g} \right) + \frac{\displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left( \Sigma _\mathrm{g} \nu r^3 \frac{\mathrm{\partial} \Omega }{\mathrm{\partial} r} \right)}{r \Sigma _\mathrm{g} \displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} r} \left( r^2 \Omega \right)}, \\ & \equiv & v_\mathrm{g,drag} + v_\mathrm{g,visc}. \end{array} \end{equation}

(C1)

f_drag is a dimensionless function of the Stokes number St and the dust-to-gas ratio ε given by

\begin{equation} f_\mathrm{drag} \equiv \frac{\epsilon }{(1 + \epsilon )^2 \mathrm{St}^{-1} + \mathrm{St}}. \end{equation}

(C2)

The first term of the right-hand side of equation (C1) is the backreacting counterpart of the usual drift term obtained for a sub-Keplerian gaseous disc (Nakagawa et al. 1986). The minus sign indicates that dust backreaction makes the gas drift outwards by angular momentum conservation. The second term is the usual inwards viscous flow. In the zero-dust-mass limit ε = 0, equation (C1) reduces to the seminal diffusion equation of a viscous disc (Lynden-Bell & Pringle 1974). From equation (C1), the conservation of mass reads

\begin{eqnarray} \frac{\mathrm{\partial} \Sigma _\mathrm{g}}{\mathrm{\partial} t} \!=\! \!-\! \frac{1}{r} \frac{\mathrm{\partial} }{\mathrm{\partial} r}\!\left( \frac{\displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left( \Sigma _\mathrm{g} \nu r^3 \frac{\mathrm{\partial} \Omega }{\mathrm{\partial} r} \right)}{\displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} r} \!\left( r^2 \Omega \right)} \right) \!+\! \frac{1}{r} \frac{\mathrm{\partial} }{\mathrm{\partial} r} \!\left( f_\mathrm{drag} \frac{r}{\Omega } \frac{\mathrm{\partial} \left(c_\mathrm{s}^{2} \Sigma _\mathrm{g} \right)}{\mathrm{\partial} r} \right). \nonumber\\ \end{eqnarray}

(C3)

The frequency Ω is the local Keplerian frequency Ω_K corrected by perturbative terms according to

\begin{equation} \Omega = \Omega _\mathrm{K} + \frac{1}{2 v_\mathrm{K} \Sigma _\mathrm{g}}\frac{\mathrm{\partial} \left(c_\mathrm{s}^{2} \Sigma _\mathrm{g} \right)}{\mathrm{\partial} r}. \end{equation}

(C4)

Looking at orders of magnitude, we denote l the typical length over which surface density gradients develop in the disc, such that

\begin{eqnarray} \frac{1}{v_\mathrm{K} \Sigma _\mathrm{g}}\frac{\mathrm{\partial} \left(c_\mathrm{s}^{2} \Sigma _\mathrm{g} \right)}{\mathrm{\partial} r} = \mathcal {O}\left(\frac{H^{2}}{rl} \right) \Omega _\mathrm{K}, \end{eqnarray}

(C5)

\begin{eqnarray} \frac{1}{\Sigma _\mathrm{g} \Omega } \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(c_\mathrm{s}^2 \Sigma _\mathrm{g} \right) = \mathcal {O}\left(\frac{H^{2}}{rl} \right) v_\mathrm{K}, \end{eqnarray}

(C6)

\begin{eqnarray} \frac{\displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left( \Sigma _\mathrm{g} \nu r^3 \frac{\mathrm{\partial} \Omega }{\mathrm{\partial} r} \right)}{r \Sigma _\mathrm{g} \displaystyle \frac{\mathrm{\partial} }{\mathrm{\partial} r} \left( r^2 \Omega \right)} = \mathcal {O}\left(\alpha \frac{H^{2}}{\min \left(r,l \right)^{2}} \right) v_\mathrm{K} . \end{eqnarray}

(C7)

When l ≳ r, equation (C7) is dominated by the global shear of the nearly Keplerian rotation profile of the disc (the first term of the right-hand side of equation C4). However, for l ≪ r the viscous velocity is dominated by a local shear induced by the correction to the orbital frequency (the second term in equation C4, see Kanagawa et al. 2015). It should be noted that an additional viscous correction has been neglected in equation (C4). This approximation is valid as long as l/r ≥ α, a condition satisfied in practice.

The relative importance of backreaction is measured via the parameter x_br defined as

\begin{equation} x_\mathrm{br} \equiv \frac{\left| v_\mathrm{g,drag} \right| }{ \left| v_\mathrm{g,visc} \right| } \simeq \frac{1}{\alpha } \frac{(\min \left(l,r \right))^2}{lr} f_\mathrm{drag}. \end{equation}

(C8)

To estimate x_br, one should distinguish two cases, depending whether a local pressure maximum has developed in the disc or not. This is similar to gap formation in dusty discs (Dipierro et al. 2016).

C2 Formation of the pressure maximum

We consider first a disc in which no pressure maximum has formed yet. This implies l ≃ r. Equation (C3) reduces therefore to

\begin{eqnarray} \frac{\mathrm{\partial} \Sigma _\mathrm{g}}{\mathrm{\partial} t} \!=\! \frac{3}{r} \frac{\mathrm{\partial} }{\mathrm{\partial} r}\left( \sqrt{r} \frac{\mathrm{\partial} }{\mathrm{\partial} r} \!\left( \sqrt{r}\nu \Sigma _\mathrm{g} \right) \right) \!+\! \frac{1}{r} \frac{\mathrm{\partial} }{\mathrm{\partial} r} \left( f_\mathrm{drag} \frac{r}{\Omega _\mathrm{K}} \frac{\mathrm{\partial} \left(c_\mathrm{s}^{2} \Sigma _\mathrm{g} \right)}{\mathrm{\partial} r} \right), \nonumber\\ \end{eqnarray}

(C9)

and an estimate of the parameter x_br is

\begin{equation} x_\mathrm{br} \simeq \frac{f_\mathrm{drag}}{\alpha }. \end{equation}

(C10)

Using equation (C2), the values of x_br as a function of St and ε for our value of α = 10⁻² are shown in Fig. C1.

Figure C1.

Parameter x_br, quantifying the importance of backreaction on the gas motion, as a function of St, for ε = 0.01, 0.1 and 1, and α = 10⁻².

Initially, ε ≃ 10⁻² and St ≪ 1 in the entire disc: f_drag ∼ ε St and x_br ≪ 1. As expected, viscosity governs the gas evolution. Once dust has piled-up such that locally, ε ∼ 1 and St ∼ 1. f_drag is of order unity, implying a large x_br ∼ 1/α: v_g,visc becomes negligible, and backreaction takes over. This simple result is counter to most of the literature to date: gas dynamics can be dominated by dust, rather than by its viscosity. This is a local effect, since f_drag is a strongly peaked function around the dust concentration.

To detail how backreaction modifies the gas structure at the location of a radial dust concentration, we integrate equation (C3) numerically between r/r_c = 0.1 and r/r_c = 100, using a typical temperature profile of T ∝ r^−0.5 and a peaked profile of |$x_\mathrm{br}=x_\mathrm{br}^0+x_\mathrm{br}^\mathrm{peak}\mathrm{e}^{-\frac{(r-r_\mathrm{c})^2}{2\Delta r_\mathrm{c}^2}}$|⁠, with |$x_\mathrm{br}^0 = 0.01$|⁠, |$x_\mathrm{br}^\mathrm{peak} = 20$| and Δr_c/r_c = 0.05. We choose a conservative value of 20 for |$x_\mathrm{br}^\mathrm{peak}$| and Δr_c is the typical width of the dust peak. We adopt for Σ_g the values of the self-similar solution at the boundaries (Lynden-Bell & Pringle 1974), and the disc is evolved on a short time-scale to study the effect of the dust radial concentration on the gas structure (this formalism is not appropriate for long-term evolution since the dust surface density varies). Fig. C2 shows the evolution of gas profiles obtained with |$x_\mathrm{br}^\mathrm{peak}=0$| and |$x_\mathrm{br}^\mathrm{peak}=20$|⁠. In the absence of a peaked dust distribution, |$x_\mathrm{br}^\mathrm{peak}=0$| and the gas evolution is viscous and is almost identical to the self-similar solution. With a peaked dust distribution, here |$x_\mathrm{br}^\mathrm{peak}=20$|⁠, mass is transferred from the inner edge of the pile-up, where a density minimum forms, to its outer edge, creating a density maximum. The disc is not dissipative enough to prevent the deformation of the gas profile, even with the large value of α ≃ 10⁻² used here. It may be noted that in the limiting case of a disc with very low dissipation (α < 5 × 10⁻³), the density maximum may induce the formation of a vortex by the Rossby wave instability (Zhu et al. 2014). For q = 0.5, we find that the gas profile becomes ‘bumpy’ for |$x_\mathrm{br}^\mathrm{peak}\gtrsim 9$|⁠. Qualitatively, the typical radial scale of the pile-up is much shorter than the radial extension of the disc. Thus, an approximate equation of the gas evolution is

\begin{equation} \frac{\mathrm{\partial} \Sigma _\mathrm{g}}{\mathrm{\partial} t} \simeq \frac{1}{\Omega _\mathrm{K}} \frac{\mathrm{\partial} }{\mathrm{\partial} r} \left( c_{\rm s}^{2} \Sigma _\mathrm{g} \right) \frac{\mathrm{\partial} f_{\rm drag}}{\mathrm{\partial} r}. \end{equation}

(C11)

From equation (C11), the typical formation time of the gas density maximum is

\begin{equation} \tau \sim \frac{r_\mathrm{c}\,\Delta r_\mathrm{c}}{H^2}\,t_\mathrm{K}, \end{equation}

(C12)

where H is the disc pressure scaleheight and t_K the Keplerian time-scale, and depends on the history of the dust evolution, according to the choice of disc and fragmentation parameters. For our conservative value of Δr_c/r_c = 0.05 and a typical H/r_c ∼ 0.1, τ ∼ 5 t_K. Equation (C11) shows that the gas density increases (respectively decreases) close to the inflection point of the function f_drag with a negative (respectively positive) slope.

Figure C2.

Gas surface density profiles after 0.1 viscous time at the dust concentration radius. The grey dotted line plots x_br, representing the dust concentration (vertical scale on the right). With |$x_\mathrm{br}^0=0.01$| and |$x_\mathrm{br}^\mathrm{peak}=0$| (orange short-dashed line), the profile is almost superimposed with the purely viscous self-similar solution (blue long-dashed line). In the profile for |$x_\mathrm{br}^0=0.01$| and |$x_\mathrm{br}^\mathrm{peak}=20$| (solid green line), a maximum forms close to the peak's inflection point with negative slope located at r/r_c = 1.065 (vertical red dot–dashed line).

Once a pressure maximum has formed in the gas, l ≪ r. Local shear generates strong viscosity that acts to spread out the local density peak (a mechanism also discussed in Yang & Menou 2010 for Rayleigh-unstable discs) and prevent a too large concentrations of gas. The pressure maximum is stable and the gas flow is almost stationary when viscosity counterbalances the backreaction. This stable situation corresponds to x_br ≃ 1. Near the pressure maximum, v_g ≃ 0. The inwards flow of particles coming from the disc outer regions sustains an effective backreaction over long times at the trap location.

APPENDIX D: CAPTIONS OF ONLINE VIDEOS

Video 1: evolution of the radial grain size distribution in the flat disc with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent up to 200 000 yr with (top) or without (bottom) backreaction of dust on gas. The colour represents the Stokes number St (left) or the ratio V_rel/V_frag (right). Initially, grain growth is only efficient in the very outer disc, where grains have the lowest V_rel, but is slow [its time-scale varies as r^3/2 (Laibe et al. 2008)]. With backreaction, the drift velocity of grains is slightly lower, delaying their drift to regions where growth is faster, they therefore grow more slowly than without backreaction. After ∼70 000 yr, their size distribution has caught up to the case without backreaction outside of ∼100 au. In the inner regions, the size distributions in both cases differ greatly because grains decouple from the gas or not when backreaction is included or not, respectively. With backreaction, most grains are piled-up in the self-induced dust trap and grow. Without backreaction, grains in the inner disc have high V_rel and fragment, replenishing the small grains reservoir. Grains thus stay well coupled to the gas and remain in the disc.

Video 2: same as Video 1 for the steep disc with V_frag = 15 m s⁻¹ and ε₀ = 1 per cent evolved up to 400 000 yr.

Video 3: evolution of the radial grain size distribution in the steep disc with V_frag = 10 m s⁻¹ and ε₀ = 1, 3 and 5 per cent (from top to bottom) up to 400 000 yr. The colour represents the Stokes number St (left) or the ratio V_rel/V_frag (right).

Video 4: evolution of the radial grain size distribution in the steep disc with ε₀ = 1 per cent (left) or 5 per cent (right) and V_frag = 10, 15, 20 and 25 m s⁻¹ (from top to bottom) up to 400 000 yr. The colour represents the Stokes number St.