© ESO 2020

1 Introduction

High spatial resolution dust continuum observations with the Atacama Large Millimeter Array (ALMA) have shown that concentric rings and gaps are common features in protoplanetary discs (see e.g., ALMA Partnership et al. 2015; Andrews et al. 2016; Isella et al. 2016; Pinte et al. 2016). Detections of rings and gaps are common in millimeter continuum emission, which traces the population of millimeter-sized pebbles (e.g., Clarke et al. 2018; Dipierro et al. 2018; Fedele et al. 2018; Long et al. 2018, 2019); however, the same patterns have also been observed in the distribution of micron-sized dust grains via scattered light observations (e.g., Ginski et al. 2016; van Boekel et al. 2017; Avenhaus et al. 2018), and in the gas surface density via observations of molecular emission (e.g., Isella et al. 2016; Fedele et al. 2017; Teague et al. 2017; Favre et al. 2019). Since these features are seen in both the solid and the gas component of the disc, they have been interpreted as a signature of planet-disc interactions (e.g., Pinilla et al. 2012; Dipierro et al. 2015; Favre et al. 2019; Fedele et al. 2018). Numerical simulations have long predicted that massive planets will open gap(s) in the disc, creating local pressure maxima at the gap edges where particles can become trapped (e.g., Papaloizou & Lin 1984; Paardekooper & Mellema 2004; Dong et al. 2015). This idea has gained relevance, as Dullemond et al. (2018) recently found evidence that at least some of the rings seen in observations are due to dust trapping in radial pressure bumps. Since the dust-to-gas ratios in these pressure bumps can become significantly higher than the global value, they are a likely site for planetesimal formation via the streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007; Bai & Stone 2010a; Yang & Johansen 2014). This is the key process that we investigate in this work.

In this paper, we explore observational consequences of the formation of gaps by embedded planets. Particularly, we investigate whether particle trapping at the edges of planetary gaps is efficient enough to trigger the streaming instability and result in the formation of planetesimals. To this end, we used a dust evolution model including both radial drift, stirring, and coagulation, and we performed first principle calculations in 1D over large spatial scales (1−500 au) and long disc evolution lifetimes (1 Myr). The main points that we want to address can be summarized as follows: (1) do planetesimals form at the edges of planetary gaps; (2) if so, how efficient is this process and how does the efficiency vary with different disc and planet parameters; (3) what does the distribution of dust and pebbles look like for the different simulations; and (4) how do these distributions compare with observations of protoplanetary discs?

We find that planetesimals do indeed form at the edges of planetary gaps. For millimeter-sized pebbles and planetary masses larger than the pebble isolation mass, essentially all pebbles trapped at the pressure bump are turned into planetesimals. In combination with fast radial drift, this results in the depletion of pebbles in the region’s interior to the outermost planet, leaving us with something that resembles a transition disc (Andrews et al. 2011). When the particle size is lowered to 100 μm, a larger dust-to-gas ratio is required to trigger the streaming instability. Because of this, the planetesimal formation efficiency drops, and at least one ring of pebbles remains visible in the otherwise empty inner disc. When the planetary mass is lowered to less than the pebble isolation mass, trapping at the gap edges becomes less efficient. Pebbles are also able to partially drift through the planetary gaps, resulting in a continuous replenishing of pebbles to the inner disc, and the pebble distribution appears as a series of weak gaps and rings at the locations of the planets.

In Sect. 2, we present our models for disc evolution, dust growth, and planetesimal formation. The numerical set up of the simulations is described in Sect. 3. In Sect. 4, we present the results for the nominal model and the parameter study. In the nominal simulation, the maximum pebble size reached by coagulation is around one millimeter; the results from simulations where the maximum grain size was artificially decreased to 100 μm are presented inSect. 5. In Sect. 6, we compare our results to observations, and in Sect. 7 we discuss the shortcomings of the model. The most important results and findings are summarized in Sect. 8. In Appendix A, we describe the method used to model particle collisions, and in Appendix B we show particle size distributions for some interesting simulations in the parameter study.

2 Theory

In our model, we used a 1D gas disc that is evolving viscously. We applied planetary torques to the disc in order to simulate gap opening by planets. For the evolution of dust particles, we used a model containing both particle growth, stirring, and radial drift. We further included a model for planetesimal formation where the conditions for forming planetesimals were derived from streaming instability simulations.

2.1 Disc model

The initial surface density profile of the disc was chosen to be that of a viscous accretion disc, $$\Sigma =\frac{{\dot{M}}_0}{3\pi \nu }\exp \left[-\frac{R}{R_{\text{out}}}\right],$$(1)

where Σ is the surface density of the gas, Ṁ₀ is the initial disc accretion rate, ν is the kinematic viscosity of the disc, R is the semimajor axis, and R_out is the location of the outer disc edge (e.g., Pringle 1981). The evolution of the surface density was solved using the standard 1D viscous evolution equation from Lin & Papaloizou (1986) for a disc that is being perturbed by a planet, $$\frac{\partial \Sigma }{\partial t}=\frac{1}{R}\frac{\partial }{\partial R}\left[3R^{1/2}\frac{\partial }{\partial R}(\nu \Sigma R^{1/2})-\frac{2\Lambda \Sigma R^{3/2}}{{(G{M}_{*})}^{1/2}}\right].$$(2)

In the above equation, t is the time, Λ is the torque density distribution, G is the gravitational constant, and M_* is the stellarmass. Equation (2) is essentially the continuity equation in cylindrical coordinates, $$\frac{\partial \Sigma }{\partial t}=\frac{1}{R}\frac{\partial }{\partial R}\left(\Sigma v_{\text{R}}R\right),$$(3)

where the radial velocity v_R has two components, which can be obtained from the comparison of the two equations. The kinematic viscosity was approximated using the alpha approach from Shakura & Sunyaev (1973), $$\nu ={\alpha }_{\text{visc}}\Omega H^2,$$(4)

where α_visc is a parameter related to the efficiency of viscous transport, $\Omega ={(G{M}_{*}/R^3)}^{1/2}$ is the Keplerian angular velocity, and H is the scale height of the disc. The scale height was calculated as $$H=\frac{c_{\text{s}}}{\Omega },$$(5)

where c_s is the sound-speed, $$c_{\text{s}}={\left(\frac{k_{\text{B}}T}{\mu m_{\text{H}}}\right)}^{1/2}.$$(6)

In the above equation, k_B is the Bolzmann constant, T is the temperature, m_H is the mass of the hydrogen atom, and μ is the mean molecular weight, which was set to be 2.34 for a solar-composition mixture of hydrogen and helium (Hayashi 1981). We used a fixed powerlaw structure for the temperature, $$T=T_{\text{const}}\times {(R/\text{AU})}^{-\zeta },$$(7)

with a radial temperature gradient of ζ = 3∕7 and a midplane temperature of T_const = 150 K at 1 au (Chiang & Goldreich 1997).

2.2 Planetary torque

The effect on the disc due to the planet is governed by the torque density distribution, Λ, which here is defined as the rate of angular momentum transfer from the planet to the disc per unit mass. For modeling of the torque density distribution, we follow D’Angelo & Lubow (2010), $$\Lambda =-F(x,\beta ,\zeta ){\Omega }_a^2{a}^2{q}^2{\left(\frac{a}{H_a}\right)}^4.$$(8)

In the above equation F is a dimensionless function, x = (R − a)∕H_a, β, and ζ are the negative radial gradients of surface density and temperature, respectively, and q is the planet-to-star mass ratio. The subscript a denotes the location of the planet. The analytic expression used for function F is $$\begin{array}{ll}F(x,\beta ,\zeta )=\hfill & \left\{p_1\exp \left[-\frac{{(x+p_2)}^2}{p_3^2}\right]+p_4\exp \left[-\frac{{(x+p_5)}^2}{p_6^2}\right]\right\}\hfill \\ \hfill & \times \tanh (p_7-p_{8}x),\hfill \end{array}$$(9)

where the parameters (p₁, …, p₈) are providedas a fit to actual simulations. Table 1 of D’Angelo & Lubow (2010) gives values for these parameters for a set of discrete values of β and ζ. As mentioned in the previous subsection, we used a fixed radial temperature gradient of ζ = 3∕7, and we chose to simplify the problem even further by also using a constant surface density gradient of β = 15∕14. For these values of β and ζ the parameters (p₁, …, p₈) take on the values listed in Table 1.

2.3 Pebble isolation mass

The pebble isolation mass (M_iso) is defined asthe mass when the planet can perturb the local pressure gradient in the midplane enough to make it zero outsidethe gap, thus creating a pressure bump. Pebbles drifting inward are trapped at the pressure bump, resulting in a locally enhanced solid density just outside the orbit of the planet (e.g., Lambrechts et al. 2014; Pinilla et al. 2016; Weber et al. 2018). The planetary masses used in our simulations are always in units of pebble isolation masses.

We used an analytical fitting formula for the pebble isolation mass, which was derived by Bitsch et al. (2018) using 3D hydrodynamical simulations of planet-disc interactions, $$M_{\text{iso}}=25M_{\oplus }\times f_{\text{fit}},$$(10)

and $$f_{\text{fit}}={\left[\frac{H/R}{0.05}\right]}^3\left[0.34{\left(\frac{\log ({\alpha }_3)}{\log ({\alpha }_{\text{visc}})}\right)}^4+0.66\right]\left[1-\frac{\frac{\partial \ln P}{\partial \ln R}+2.5}{6}\right].$$(11)

In the above equation α₃ = 0.001.

The pebble isolation mass is extremely dependent on the gap depth, which is known to vary between 1D and multidimensional simulations (e.g., Lin & Papaloizou 1986; Kanagawa et al. 2015; Hallam & Paardekooper 2017). In general, 1D simulations produce narrower and deeper gaps than their higher dimensional analogs. This means that the mass required to reach a radial pressure gradient of zero is smaller in 1D simulations than in 3D simulations. In other words, the pebble isolation masses derived by Bitsch et al. (2018) are significantly higher than the pebble isolation mass obtained in our simulations using the tabulated torques of D’Angelo & Lubow (2010).

We performed our own 1D simulations to calculate the pebble isolation mass, and found during comparison that the pebble isolation masses obtained in 1D and 3D simulations are related by a scalar factor, which only seems to depend on α_visc. Approximate values of this scalar, which we denote k_3D∕1D, can be found in Table 2 for a range of values of α_visc. As an example: if α_visc = 0.01, then the pebble isolation mass at 11.8 au is 59.6 M_⊕, according to Eq. (10) and by using the values quoted for T_const, β, and ζ. The planetary mass required to obtain a zero pressure gradient outside a planet orbiting at the same semimajor axis in our simulations is 59.6∕1.5 = 39.7 M_⊕. Figure 3 of Johansen et al. (2019) shows a similar systematic difference between the 1D and the 3D gap depth.

In our simulations, to avoid working with an artificially low pebble isolation mass due to the 1D approach, we modified the magnitude of the torque density distribution in the following way, $$\Lambda \to \frac{\Lambda }{k_{3\text{D}/1\text{D}}^2}.$$(12)

In other words, for α = 0.01 we simply divided Eq. (8) by 1.5², and then the pebble isolation masses obtained from Eq. (10) are correct even in 1D. The power of two is there because the torque density is proportional to the planetary mass square (Eq. (8)).

2.4 Dust evolution

We adopted the approach of Lagrangian super-particles for the solid component of the disc. Each super-particle represents multiple identical physical solid particles, and each super-particle i has its own position x_i and velocity v_i. The particlevelocity was taken as the sum of the drift velocity and the turbulent velocity, the algorithms for which are described below.

2.4.1 Drift velocity

The radial drift velocity of a dust particle in a disc that is accreting gas is $$v_{\text{D}}=-\frac{2{\tau }_{\text{s}}}{1+{\tau }_{\text{s}}^2}\left(\eta v_{\text{kep}}-\frac{1}{2{\tau }_{\text{s}}}{v}_{\text{R}}\right),$$(13)

where η is the difference between the azimuthal gas velocity (v_θ) and the Keplerian velocity (v_kep), v_R is the gas velocity in the radial direction, and τ_s is the Stokes number, which is sometimes referred to as the dimensionless stopping time (Nakagawa et al. 1986; Guillot et al. 2014). The η parameter is directly related to the pressure gradient of the disc as $$\eta =-\frac{1}{2}{\left(\frac{H}{R}\right)}^2\frac{\partial \ln P}{\partial \ln R}.$$(14)

The Stokes number for a particle in the Epstein regime is $${\tau }_{\text{s}}=\frac{s{\rho }_{•}}{H\rho };$$(15)

it is important to note that the factor $\sqrt{\pi /8}$ in Eq. (6) from Carrera et al. (2015) is not included in this work and would not have changed the results significantly. In the above equation s is the particle radius, ρ is the gas density in the midplane (related to the gas surface density through ρ = Σ∕(2πH)), and ρ_∙ is the solid density. We adopt the value ρ_∙ = 1000 kg m⁻³ throughout this work. This density could approximately represent the density of icy pebbles with a significant porosity.

2.4.2 Turbulent velocity

Solid particles in the protoplanetary disc experience a drag force due to the turbulent motion of gas in the disc. The resultant turbulent diffusion of the particles can be modeled as a damped random walk, which we implemented by using the algorithm from Ros et al. (2019). They calculated the turbulent diffusion coefficient (D) by applying a force acceleration (f) to the particles on a time-scale τ_for, and damping the turbulent velocity on the correlation time-scale τ_cor. The forcing time-step τ_for was set to equal the time-step of the simulation, and the correlation time is the approximate time-scale over which a particle maintains a coherent direction, which was calculated as the inverse of the Keplerian angular velocity (τ_cor = Ω⁻¹). As an addition to the algorithm from Ros et al. (2019), the falling of D with Stokes number was implemented here following Eq. (37) in Youdin & Lithwick (2007). We used a value of one for the dimensionless eddy time.

2.4.3 Particle collisions

Our model of particle collisions follows the results of Güttler et al. (2010). They combine laboratory collision experiments and theoretical models to show that collisions among dust particles in the disc lead to either sticking, bouncing, or fragmentation. The outcome is determined by the mass of the projectile (the lightest of the colliding particles) and the collision velocity, which is calculated as the sum of the relative speed from drift, Brownian motion, and turbulent motion. The result of the collision also varies depending on the mass ratio of the projectile to target particle and on whether the particles are porous or compact. In our simulations, we limited ourselves to porous particles, and drew the line between equal-size particles and nonequal-size particles at a mass ratio of 10 (using effectively only the two upper panels of Fig. 11 in Güttler et al. 2010). If the outcome is sticking, then the target mass was either doubled, or multiplied by (1 + M_projectile∕M_target) if the total mass in the projectile particle M_projectile was smaller than the total mass in the target particle M_target. For fragmentation, we set all of the target and projectile particles to the mass of the projectile. For a complete description of the collision algorithm, see Appendix A.

2.5 Planetesimal formation

Carrera et al. (2015) performed hydrodynamical simulations of particle-gas interactions to find out under which conditions solid particles in the disc come together in dense filaments that can collapse under self-gravity to form planetesimals. By doing so, they mapped out for which solid concentrations and particle Stokes number filaments emerge. This map was revised by Yang et al. (2017) who expanded on the investigation by using longer simulation times and significantly higher resolutions. The critical curves on the map for when the solid concentration is large enough to trigger particle clumping are given by Yang et al. (2017) as $$\log Z_{\text{c}}=0.3{(\log {\tau }_{\text{s}})}^2+0.59\log {\tau }_{\text{s}}-1.57({\tau }_{\text{s}}>0.1),$$(16)

and $$\log Z_{\text{c}}=0.1{(\log {\tau }_{\text{s}})}^2+0.20\log {\tau }_{\text{s}}-1.76({\tau }_{\text{s}}<0.1),$$(17)

where Z_c = Σ_solid∕Σ_total and the logarithm is with base 10. These equations were derived for a laminar disc model; however, unless the degree of turbulence is very high, they may also be valid for nonlaminar discs (Yang et al. 2018). For example, Yang et al. (2018) find a critical solid-to-gas ratio of 2% when using τ_s = 0.1 particles and a vertical turbulence strength of 10⁻³, driven by density waves excited by the magnetorotational instability in the turbulent surface layers.

2.5.1 Pressure dependence

The map from Yang et al. (2017) determines whether or not the streaming instability forms filaments based on the solid abundanceand the Stokes number; however, the degree of clumping is also strongly dependent on the radial pressure gradient (Bai & Stone 2010b; Abod et al. 2019; Auffinger & Laibe 2018). Simulations by Bai & Stone (2010b) show that the critical solid abundance required to trigger particle clumping via the streaming instability monotonically increases with the radial pressure gradient. In other words, planetesimal formation is most likely to occur inside pressure bumps where the pressure gradient is small, a result which has also been obtained from a linear analysis (Auffinger & Laibe 2018). From Bai & Stone (2010b), it appears that the critical solid abundance is roughly linearly dependent on the radial pressure gradient. In order to catch this dependency, we scaled the metallicity threshold from Yang et al. (2017) with the local radial pressure gradient divided by the background pressure gradient. Here we made the assumption that the background pressure gradient in our simulations is the same as in Yang et al. (2017). This is a simplification; however, a comparison of simulations using the different background pressure gradients show that the results are unaffected.

When the pressure gradient is exactly zero, there is no particle drift, meaning that the streaming instability is formally absent. However, there is another planetesimal formation mechanism that is still active, which we have not not included in our simulations, namely secular gravitational instability (e.g., Youdin 2011; Takahashi & Inutsuka 2014). Recently, Abod et al. (2019) performed simulations of planetesimal formation for various pressure gradient conditions (including a zero pressure gradient). They find that many results and conclusions obtained in their study of planetesimal formation via the streaming instability are also valid in the case of a zero pressure gradient. Motivated by this result, we chose to use the same mechanism for planetesimal formation at a pressure gradient of zero as we otherwise did.

2.5.2 Code implementation

We implemented planetesimal formation in the code in the following way: (1) we calculated the solid abundance in each grid cell, excluding the already formed planetesimals; (2) we calculated the local pressure gradient in each grid cell; (3) we scaled the metallicity threshold Z_c up and down linearly by the found local pressure gradient divided by the background pressure gradient; (4) we calculated the mean Stokes number in each grid cell; (5) if the criterion for planetesimal formation is reached, then we set the radius of the first super-particle in the grid cell to 100 km; (6) we repeated the process until the criterion was no longer met. We note that the planetesimal size is arbitrary since we do not follow the dynamical evolution of the planetesimals after their formation.

The algorithm described above implies that every time the conditions for the streaming instability are met, planetesimals form. This can be thought of as a limiting maximum case for the planetesimal formation efficiency, and a discussion regarding how accurate this algorithm is can be found in Sect. 7.3.

2.5.3 An alternative model for planetesimal formation

The criterion used in this paper for planetesimal formation via the streaming instability is not the only one that exists. Another commonly used criterion is that the dust-to-gas ratio in the midplane has to be larger than unity. Such a planetesimal formation model is used in Stammler et al. (2019), for example, who conducted a similar study of planetesimal formation in pressure bumps. We implemented this planetesimal formation criterion, which we refer to as the midplane model, in the code and make a comparison of the two planetesimal formation models in Sect. 4.4. For code implementation, we used the same algorithm as above, except that it was applied to the midplane density ratio rather than to the column density ratio.

3 Numerical set up

The code we used for simulations is called PLANETESYS, and it is a modified version of the Pencil Code (Brandenburg & Dobler 2002), which was designed for highly parallel calculations of the evolution of gas and dust particles in protoplanetary discs. The code was developed under the ERC Consolidator Grant “PLANETESYS” (PI: Anders Johansen) and this paper together with the recent paper by Ros et al. (2019) represent the first publications using this tool.

The evolution of the surface density was solved using a first order finite difference scheme with an adaptive time-step. The disc stretches from 1 to 500 au with R_out = 100 au and it was modeled using a linear grid with 4000 grid cells. For the inner boundary condition, we copied the values of the adjacent cells, and for the outer boundary condition we simply set the density to zero at the outer disc edge. This provides the right solution to the viscous disc problem, and fits the analytically derived surface density profile well out to at least 200 au. We used a stellar mass of 1 M_⊙ throughout the simulations and an initial disc accretion rate of Ṁ₀ = 10⁻⁷ M_⊙ yr⁻¹. The accretion rate drops to 2 × 10⁻⁸ M_⊙ yr⁻¹ after 1 Myr, as material drains onto the star. Most simulations were run on 40 cores to speed up the calculations. The typical wall time was 90 h.

Three planets of fixed masses and semimajor axes (except in simulation migration) are included in the simulations, and they were inserted at semimajor axes corresponding to the locations of the major gaps in the disc around the young star HL Tau: at 11.8, 32.3, and 82 au (Kanagawa et al. 2016). The solid population of the disc is represented by 100 000 superparticles (approximately 25 per grid cell). The superparticles were initially placed equidistantly throughout the disc with a radius of 1 μm. The mass of each superparticle was set so as to yield a constant solid-to-gas ratio (also referred to as the metallicity) across the disc. Planetesimal formation was initiated after some time t_plan, which was varied between the simulations depending on how much time it takes for the planets to clear most of their gaps from dust. Furthermore, we only allowed for planetesimal formation interior to 200 au. This is because our numerical surface density profile starts to diverge from the analytically derived one beyond a few hundred astronomical units; however, since we are only interested in the inner ~100 au where the planets are located, this does not affect the results. Finally the system is evolved for 1 Myr. This long running time is a major motivation for simulating in 1D only.

In the nominal model (simulation #1 in Table 3), we used a turbulent viscosity of 10⁻², a turbulent diffusion of 10⁻³, and an initial solid-to-gas ratio of 0.01. We set the planetary masses (M_p) to two times their respective pebble isolation mass, and we kept the planets at a fixed position (V_p= 0). In this simulation, planetesimal formation was initiated after 5000 yr. To explore how the above mentioned parameters affect the planetesimal formation efficiency and the distribution of dust and pebbles, we conducted a parameter study. The parameter values used in the different simulations can be found in Table 3. In simulations #2–#5, we varied the planetary masses; in simulations #6–#8, we lowered the value of the viscosity parameter and turbulent diffusion; in simulation #9, we increased the initial solid-to-gas ratio in the disc; and in simulation #10, we let the planets migrate radially inward with a constant velocity.

4 Results

Results on the disc structure, particle distribution, and efficiency of planetesimal formation in the nominal model are presented in Sect. 4.1. In Sect. 4.2 we show how these results change when we no longer include the pressurescaling for the streaming instability, and when planetesimal formation is removed completely. In Sect. 4.3, we varied different disc and planet parameters and investigate how it affects the results. Finally, in Sect. 4.4, we make a comparison between our model for planetesimal formation and the midplane model.

4.1 Nominal model

The parameters used in the nominal model are found in simulation #1 of Table 3. The evolution of the normalized gas surface density across the disc is shown in the top panel of Fig. 1, and the evolution of the solid-to-gas surface density ratio is plotted in the bottom panel. The solid component was divided into planetesimals and dust+pebbles. The evolution of the particle size distribution and the Stokes number is shown in Fig. 2.

From Fig. 1, it is clear that the part of the disc that is closer to the star than the outermost planet becomes depleted of dust and pebbles already after a few hundred thousand years. Meanwhile, the region’s exterior to the outermost planet maintains a high solid-to-gas ratio for much longer, and at the end of the simulation, it is still at 25% of the initial value. To understand why this is the case, we look at the evolution of the particle sizes in Fig. 2. Considering the inner part of the disc, particle collisions are frequent and quickly result in the formation of millimeter-sized pebbles, which drift toward the star quickly.The same process results in smaller particle sizes and takes more time in the outer part of the disc. Another feature of the particle size distribution is that it is bimodal. The reason for this is that except for very low relative velocities, coagulation is only possible if the projectile is much less massive than the target (compare the two upper panels of Fig. 11 in Güttler et al. 2010). Hence any coherent size distribution evolves to a bimodal distribution as the largest particles grow in mass, while the small are stuck. The bimodal size distribution nevertheless collapses with time to a narrower size distribution as the remaining small particles are finally swept up by the larger pebbles.

The planetary masses used in the nominal simulation lead to relatively deep gaps, which act as hard barriers that stall all particles with a Stokes number that is larger than a certain critical value (see e.g., Pinilla et al. 2012; Bitsch et al. 2018; Weber et al. 2018). The pebbles formed by coagulation beyond the middle planet have Stokes numbers that are well above this critical value, and they are therefore efficiently trapped at the planetary gap edges. Since no particles make it past the gaps, there is no replenishing of the interplanetary disc regions. Combined with fast radial drift, this is the reason why the solid-to-gas ratio in these regions drops to less than 0.001 after only 500 000 yr.

From Figs. 1 and 2, it is evident that planetesimal formation takes place almost exclusively in narrow rings at the location of the gap edges and inside the planetary gaps. These are locations which correspond to places where the pressure gradient is close to zero, that is, places where there is a pressure bump. When the pressure gradient is close to zero the critical density required for the streaming instability to form filaments becomes very low, and because of this essentially everything that enters the pressure bump is turned into planetesimals. So since everything that goes into the pressure bump is turned into planetesimals, and there is no replenishing of the interplanetary regions, the part of the disc that is closer to the star than the outermost planet becomes empty of dust and pebbles. As is seen in Sect. 4.2 Fig. 3, efficient planetesimal formation at the gap edge of the innermost planet is the reason why no pebbles can make it past that gap to the innermost disc region.

The narrowness of the rings in which planetesimals form is also an effect caused by the dependency of the streaming instability on the pressure gradient. The magnitude and steepness of the pressure bump increases with increasing planetary mass, and for a planetary mass of two pebble isolation masses, the pressure gradient is only close to zero in two very narrow regions. In between these regions, the pressure gradient becomes very high, resulting in that very large critical densities are required to form filaments, which is why no planetesimals form in between the location of the gap edges and inside the planetary gaps.

The pressure maxima inside the planetary gaps are unstable equilibrium points due to the fact that tiny displacements from this point are amplified with time due to divergent drift, but we still formed some planetesimals there because we turned on planetesimal formation before the gaps had been entirely cleared of dust and pebbles. Later on during the simulation, all planetesimals form at the gap edges. A histogram showing the total amount of planetesimals that have formed per semimajor axis bin can be viewed in the uppermost panel of Fig. 4. A total of 280 Earth masses of planetesimals have formed at the end of the nominal simulation. The total amount of mass in planetesimals and dust+pebbles at the end of all the simulations is provided in Table 4 for each ring.

Fig. 1 Top panel: evolution of the gas surface density, which is normalized against the initial gas surface density, across the protoplanetary disc for the nominal model with planets (solid lines) and without planets (dotted lines). The vertical dashed linesmark the semimajor axes of the planets, and they coincide with the locations of the three major gaps in the disc around HL Tau (Kanagawa et al. 2016). Bottom panel: evolution of the solid-to-gas surface density ratio across the protoplanetary disc for the nominal model. The solid component is divided into planetesimals (marked with gray) and dust+pebbles. Planetesimals form in narrow rings at the location of the gap edges and inside the planetary gaps (the amount of planetesimals formed inside the gaps is negligible). The interplanetary regions are depleted of dust and pebbles within a few hundred thousand years.

Fig. 2 Top panel: size distribution of particles in the protoplanetary disc at different times during disc evolution for our nominal model. The semimajor axes of the formed planetesimals are indicated at the top of the plot. Particles were initialized with an equal spacing all over the disc and with a radius of 1 μm. Efficient coagulation and high drift velocities in the inner disc result in a depletion of the interplanetary regions after only a few 100 000 yr. These processes occur on a longer timescale in the outer disc. Bottom panel: particle Stokes number versus semimajor axis for the same data as in the top panel.

4.2 The cases with no pressure scaling and no planetesimal formation

How efficiently pebbles are converted into planetesimals at the gap edges has a major effect on the distribution of dust and pebbles in the disc. In the nominal model, planetesimals form whenever the critical density required for the streaming instability to form filaments has been reached. Together with the linear pressure scaling this can be thought of as a maximum limiting case for planetesimal formation. If the dependency of the streaming instability on the pressure gradient is weakened or removed, the critical density increases, resulting in more dust and pebbles being left in the rings. It would also result in planetesimals forming in wider regions around the pressure bump. The case with no planetesimal formation at all is of course the minimum limiting case for planetesimal formation, and the real efficiency should be somewhere in between the minimum and maximum case.

In Fig. 3, the evolution of the solid-to-gas surface density ratio is shown as 2D surface density plots (here referred to as disc images) for the following: the nominal model; the nominal model with no pressure scaling for planetesimal formation; and the nominal model without planetesimal formation. In the case without planetesimal formation (simulation noPlanetesimal), the disc would be seen to have very bright rings at the positions of the outer gap edges, which is similar to AS 209 (Fedele et al. 2018). The amount of dust and pebbles in each ring is written in Table 4, and it is roughly a hundred Earth masses for the two outermost rings, which corresponds to 3–11 Earth masses per astronomical unit. In contrast, the inner ring only has a few Earth masses of dust and pebbles. Comparing that to the tens of Earth masses of planetesimals that form in the ring in the nominal model, it suggests that millimeter-sized pebbles with low Stokes numbers (i.e., close to the star) can drift past our planetary gaps quite efficiently, but that the efficient planetesimal formation in the nominal model prevents them from doing so.

When planetesimal formation without pressure scaling is included (simulation noPscaling) the solid-to-gas surface density ratio drops drastically in the rings. The two rings, which are still clearly visible, are also much thinner than in the case with no planetesimal formation. The amount of dust and pebbles left in these rings are now on the order of a few Earth masses, which corresponds to a quarter of an Earth mass per astronomical unit. When pressure scaling is applied (i.e., the nominal model) the amount of early transport through the gaps decreases, resulting in a faster pebble depletion. After 1 Myr, all rings that were previously visible interior to the outermost planet have disappeared, and we are left with a cavity of roughly 100 au in size in the dust and pebble distribution. A comparison between the amount of planetesimals that form, as well as where they form, in simulations with and without the pressure scaling can be seen in the top panel of Fig. 4. The differences between the simulations in Fig. 3 tells us that two discs that are similar in all aspects except for in the efficiency of planetesimal formation could come across as completely different in observations (see discussion in Sect. 7.3 on the efficiency of the streaming instability).

Fig. 3 2D symmetric disc images of the evolution of the solid-to-gas surface density (excluding the formed planetesimals) for three different versions of the nominal model: the nominal model (top row), the nominal model with planetesimal formation but with no dependency on the pressure gradient (middle row), and the nominal model without planetesimal formation (bottom row). In the nominal model, essentially everything that enters the pressure bump is converted into planetesimals, resulting in a large cavity in the distribution of dust and pebbles. When the pressure dependence is neglected, the critical density required to trigger the streaming instability increases, and we see some rings in the dust and pebble distribution. When planetesimal formation is removed completely, we are left with three rings in which the dust and pebble density is very high.

Fig. 4 Histogram showing the total amount of mass in planetesimals that has formed at different locations in the disc after 1 Myr for each simulation in the parameter study. The dotted vertical lines indicate the semimajor axes of the planets. When the pressurescaling is removed (simulation “noPscaling”), the critical density required for forming planetesimals at the gap edge increases, resulting in less planetesimal formation. The amount of planetesimal formation does not vary a lot between the simulations where the planetary mass is larger than the pebble isolation mass since millimeter pebbles cannot drift past the gaps and planetesimal formation is extremely efficient where the pressure gradient is close to zero. For planetary masses that are lower than the pebble isolation mass, the amount of planetesimals that form rapidly decreases with decreasing planetary mass. Lowering the viscosity parameter and the turbulent diffusion results in faster drift velocities in the viscously expanding part of the disc. Because of this, more pebbles reach the outermost gap edge and are turned into planetesimals. Increasing the metallicity (simulation highMetal) results in sporadic planetesimal formation in the interplanetary regions. When a constant migration speed is added to the planets, the location of the gap edges are shifted inward with time, resulting in the region where planetesimals form to shift inward accordingly.

4.3 Parameter study

To explore how different parameters affect the planetesimal formation efficiency and the distribution of dust and pebbles in the disc, we conducted a parameter study. The values of the parameters which we investigate can be found in Table 3. A histogram showing the total amount of mass locked up in planetesimals at different semimajor axes after 1 Myr is presented inFig. 4 for all simulations in the parameter study. The evolution of the dust and pebble surface density relative to the gas surface density is presented in Fig. 5 as 2D symmetric disc images for the same simulations. The total amount of mass in planetesimals and dust+pebbles in each ring can be found in Table 4. The evolution of the particle size distribution for a few interesting simulations in the parameter study is shown in Fig. B.1.

4.3.1 Planetary mass

The planetary mass is the main controller of the width and depth of the planetary gap as well as the radial pressure gradient. The width of thegap determines where pebbles are trapped, and thus the location of planetesimal formation. The planetesimal formation efficiency is strongly related to the strength of the pressure maxima, both via the scaling of the streaming instability with the radial pressure gradient and via the efficiency of particle trapping. The distribution of dust and pebbles for simulations with varying planetary masses is shown in Fig. 5 (simulations nominal, 0.50 M_iso, 0.75 M_iso, 1 M_iso and 3 M_iso).

For the simulations 0.50 M_iso and 0.75 M_iso, we used a planetary mass that is lower than the pebble isolation mass. In these simulations dust and pebbles are partly transported through the planetary gaps, resulting in a continuous replenishing of the interplanetary regions and the region interior to the innermost planet. There is a small pile-up of pebbles at the gap edges, resulting in a gap-and-ring-like structure. As can be seen in Fig. 4, the amount of planetesimals that form decreases quickly with decreasing planetary mass. By lowering theplanetary mass to 25% below the pebble isolation mass, the amount of pebbles that are converted into planetesimals is halved. Since the radial pressure gradient is shallow around the gap edges, planetesimals form in relatively wide regions as opposed to in two narrow rings, such as in the nominal model. A plot of the size distribution of particles for simulation 0.75 M_iso can be viewed in Fig. B.1.

In the simulations nominal, 1 M_iso and 3 M_iso, the planetary masses are two, one, and three times the pebble isolation mass. For these cases, the amount of mass converted into planetesimals does not change with increasing planetary mass. This is due to two reasons: (1) the planetary gaps act as hard barriers, which prevent any pebbles from passing through; (2) there are pressure maxima outside all gaps, which efficiently turn most or all pebbles into planetesimals. The locations where planetesimals form nevertheless do change a bit, since the widened gap and the steepened pressure gradient result in planetesimals forming further away from the planet and in narrower regions. In summary, all simulations with planetary masses equal to or above the pebble isolation mass appear as discs with large central cavities, with the only difference being a slight dependence of the width of the cavity and the width of the rings on the planetary mass.

4.3.2 Viscous and turbulent α

The viscosity parameter α_visc governs the gas accretion rate onto the central star, and it also enters the particle drift equation via the radial gas velocity v_R. The turbulent diffusion α_turb governs the turbulent speed of particles, which in turn affects the frequency of particle collisions as well as the collision velocities. The values inferred from observations of these parameters vary a lot in the literature. In Pinte et al. (2016), they find a value of a few times 10⁻⁴ for the turbulent diffusion in the disc around HL Tau. They obtained this value by assuming a standard dust settling model, varying the amount of turbulent diffusion and comparing the resulting millimeter dust scale heights to observations. Pinte et al. (2016) further report an upper limit to the viscosity parameter of 10⁻² for the HL Tau disc, a value which was calculated by using an estimate of the disc accretion rate from Beck et al. (2010). Other examples are Flaherty et al. (2017) who report an upper limit of 0.003 for the turbulent diffusion in HD 163296 and Flaherty et al. (2018) who report an upper limit of 0.007 for the viscosity parameter in TW Hya.

In the nominal model, we used a value of 10⁻² for the viscosity parameter and 10⁻³ for the turbulent diffusion. In simulation lowVisc, the viscosity parameter was lowered from 10⁻² to 10⁻³. In the simulation lowTurb, the turbulent diffusion was decreased by a factor of 10 from 10⁻³ to 10⁻⁴. Finally in simulation lowViscTurb, the viscosity parameter was decreased by another order of magnitude to 10⁻⁴, while the turbulent diffusion was kept at the same level as in simulation lowTurb. When lowering the viscosity parameter, the initial disc accretion rate is reduced accordingly, ensuring the same initial disc mass.

Lowering the viscosity parameter in simulation lowVisc results in lower gas accretion rates onto the star, and thus the gas disc evolves on much longer time-scales. Because of this the disc does not expand as much, resulting in a smaller disc size. Another effect on thestructure of the gas disc is that there is a more pronounced pile-up of gas at the inner and outer edges of the planetary gaps. Looking at Fig. 5 and comparing simulation lowVisc to the nominal simulation, we see that lowering theviscosity parameter results in a slower clearing of the gaps and interplanetary regions. We also find that pebbles in the viscously expanding part of the disc drift inward more quickly. A smaller viscosity parameter results in the velocity component directed outward becoming lower, making it easier for particles to drift inward. The faster drift velocities result in more pebbles that reach the outermost planetary gap edge and that are then turned into planetesimals. This is the reason why slightly more planetesimals are formed when using this simulation compared to the nominal one.

In simulation lowTurb, the lowering of the turbulent diffusion results in fewer collisions. Because the particle Stokes numbers in our simulations are low, this also results in lower collisional velocities (if the particle Stokes numbers would have been larger so that the particles start to sediment toward the midplane, this would not have been the case, as the decreasein turbulent velocity would have been compensated for by a decrease in dust scale height). Fewer collisions lead to slower particle growth; however, slower collisional velocities result in a larger maximum pebble size. These larger particles obtain higher drift velocities, resulting in many more pebbles reaching the gap edge of the outermost planet. This is the reason why we see a wide and bright ring in the solid-to-gas surface density at 1 Myr for simulation lowTurb in Fig. 5. The larger particle sizes also require a smaller critical density to trigger the streaming instability, allowing for planetesimals to form further away from the pressure maxima (see Fig. 4).

In simulation lowViscTurb, a fast particle drift toward the star in the outer disc now results in all of the pebbles that still remain at the end of the simulation to be concentrated in a narrow bright ring just beyond the outermost planet. Since more pebbles reach the outermost planetary gap edge, the amount of planetesimal formation at that location has increased compared to the amount in the simulations lowVisc and lowTurb. The gas pile-up at the gap edges is also more pronounced, causing some particles to become trapped at the inner gap edges and leading to some planetesimal formation there (see Fig. B.1 for a plot of the particle size distribution for simulations lowTurb and lowViscTurb).

Fig. 5 2D symmetric disc images of the evolution of the solid-to-gas surface density excluding planetesimals for all simulations in the parameter study. In the nominal simulation where we used a planetary mass of two times the pebble isolation mass fast pebble drift, little or no transport through the planetary gaps and efficient planetesimal formation at the gap edges result in that the part of the disc that is closer to the star than the outermost planet gets depleted of dust and pebbles. In the simulations where the planetary mass is lower than the pebble isolation mass (simulations 0.50 Miso and 0.75 Miso), increased transport past the gaps and less planetesimal formation results in a gap-and-ring-like structure. The width of the region where dust and pebbles are trapped depends on the planetary mass. Therefore, when theplanetary mass is changed to one or three times the pebble isolation mass (simulations 1 Miso and 3 Miso), and we compare to the nominal model, the width of the rings becomes narrower or larger, respectively. Except for this, the results are the same as in the nominal model. In simulation lowVisc, the viscosity parameter is lowered, resulting in a slower clearing of the gaps and interplanetary regions, as well as faster radial drift in the viscously expanding part of the disc. When the turbulent diffusion is lowered (simulation lowTurb), the collisional velocities decrease, resulting in larger particles which drift faster toward the star. This causes the bright ring that can be seen beyond the outermost planet at the end of the simulation. In simulation lowViscTurb, the turbulent diffusion is kept at the same level as in simulation lowTurb, but the viscosity parameter is lowered by an extra order of magnitude compared to simulation lowVisc. The combinationresults in that essentially all solids in the outer disc reach the outermost planetary gap before the end of the simulation, causing the narrow bright ring seen in the dust-to-gas ratio. In simulation highMetal, the initial solid-to-gas ratio in the disc is increased to 2%, which is a change that does not have a big effect on the appearance of the disc at the end of the simulation. Finally in the last simulation, the planets were given a constant radial velocity directed toward the star (simulation migration), resulting in a smaller radius of the cavity.

4.3.3 Metallicity

In simulation highMetal, the initial solid-to-gas ratio was doubled, resulting in an increase in the total amount of formed planetesimals. Since the initial solid abundance is closer to the critical value required for the streaming instability to operate, random local concentrations of pebbles now result in some planetesimal formation in the interplanetary regions. However, since the amplitude of the fluctuations caused by diffusion would be much smaller if a physical number of particles were used, this effect is purely numerical.

4.3.4 Planet migration

In the final simulation of the parameter study (simulation migration) all three planets were given a constant velocity directed toward the star. The migration speed was set to 6.3 AU Myr⁻¹, and thus the semimajor axes of the planets after 1 Myr are 5.48 au, 25.98 au, and 75.676 au. This migration speed is not high enough to significantly perturb the shape of the gap. For faster migrating planets, hydrodynamical studies have shown that the impact on the structure of the disc can be large (see e.g., Li et al. 2009; Meru et al. 2019; Nazari et al. 2019).

Since the migration speed in our simulation is both low enough to preserve the shape of the gap and significantly lower than the particle drift velocity, the amount of planetesimal formation does not change. The location where they form does, however, change. As the pressure maximum moves inward, so does the region of planetesimal formation. For the pebble distribution, the only thing that changes is that the radius of the cavity shrinks with time.

Fig. 6 Histogram showing the total amount of mass in planetesimals that have formed at different locations in the disc after 300 kyr, for the nominal model and the midplane model. The amount of planetesimals formed around the gap edges of the two outermost planets are similar in both models; however, this is not the case at the innermost planetary gap edge. The settling toward the midplane is not efficient enough to counteract the stirring by turbulence in this region, and thus a dust-to-gas density ratio of unity in the midplane is never reached.

4.4 Two planetesimal formation models

The planetesimal formation model used in our simulations was derived from hydrodynamical simulations of particle-gas interactions by Carrera et al. (2015) and Yang et al. (2017). This model tells us whether or not filaments emerge based on the combinationof the particle Stokes number and surface density. As mentioned at the end of Sect. 2.5, some authors use a criterion based on the midplane dust-to-gas ratio instead (the midplane model). We compare this criterion with the one used in this paper by performing two simulations that are identical in all aspects except for the planetesimal formation model. We used linear pressure scaling in both simulations. The results of the comparison are presented in Figs. 6 and 7.

The two models produce very similar amounts of planetesimals at the outermost planetary gap-edge, while there are nosignificant differences in the surface density profiles. Around the second planetary gap-edge, the midplane model produces slightly less planetesimals than the nominal model, resulting in a ring of pebbles that does not exist in the nominal model. At the innermost planetary gap-edge the midplane model fails at producing any planetesimals. This is because the settling toward the midplane is not efficient enough to counteract the stirring by turbulence. Because of this, a large population of dust and pebbles remain at this location. Transport of solids through the planetary gap further results in that the innermost disc region does not get depleted of solids within 300 kyr, which is the case in the nominal model.

In the simulations above, we have taken the dependency of the streaming instability with the pressure gradient into account. If the linear pressure scaling were to be removed, the midplane model does not produce any planetesimals at all. This is because the millimeter-sized pebbles created through coagulation are stirred too much by turbulence, and thus a midplane dust-to-gas ratio of unity is never reached.

5 Lowering the particle size to 100 μm

The dust growth model of Güttler et al. (2010) employed in this work, which is based on a combination of laboratory collisionexperiments and theoretical models, results in the formation of millimeter-sized pebbles in the inner part of the disc with decreasing grain sizes as the semimajor axis increases. These sizes are slightly smaller than the grain size estimates that were obtained from the spectral index of the dust opacity coefficient at millimeter and submillimeter wavelengths, which report maximum grain sizes between 1 millimeter in the outer disc and a few centimeters in the inner disc (e.g., Ricci et al. 2010; Pérez et al. 2012; ALMA Partnership et al. 2015; Tazzari et al. 2016). The estimates based on the spectral index of the dust opacity coefficient nevertheless are not in agreement with the maximum grain sizes that were obtained from observations of polarized emission due to self-scattering, which are consistently around 100 μm, that is, smaller thanthe pebbles in our simulations (e.g., Kataoka et al. 2017; Hull et al. 2018; Ohashi et al. 2018; Mori et al. 2019).

Even when applied to the same source, the maximum grain sizes obtained from the different methods are inconsistent. As an example, weconsider the disc around HL Tau. Carrasco-González et al. (2019) calculated the maximum grain size in HL Tau by fitting the millimeter spectral energy distribution without any assumptions about the optical depth of the emission. By also including the effects of scattering and absorption in the dust opacity, they obtained a maximum grain size ofa few millimeters. Kataoka et al. (2017) instead estimated the maximum grain size in HL Tau to be 100 μm from observations of millimeter-wave polarization.

If the maximum grain size was in fact around 100 μm, as suggested by observations of millimeter-wave polarization, it could have a large impact on our results. In the model of Güttler et al. (2010), particle collisions result in the formation of millimeter-sized particles; however, recently Okuzumi & Tazaki (2019) have shown that dust growth models can result in 100 μm-sized particles if the particles are covered by nonsticky CO₂ ice. By incorporating the composition-dependent sticking into a model of dust evolution, they were able to successfully reproduce the polarization pattern seen in the disc around HL Tau.

Decreasing the particle size results in smaller particle Stokes numbers; however, it should be mentioned that this is not the only way to obtain low particle Stokes numbers. If gas discs are in fact much more massive than the minimum mass solar nebula, then millimeter-sized particles would have smaller Stokes numbers than they have in our simulations (Powell et al. 2019). In such a disc, particle drift would be slower and the critical density required for the streaming instability to form filaments would be higher.

Fig. 7 Evolution of the solid-to-gas surface density ratio across the protoplanetary disc for the nominal model (top panel) and the midplane model (bottom panel). The solid component includes dust and pebbles, whereas planetesimals are excluded. The two models produce very similar amounts of planetesimals around the outermost planetary gap-edge, resulting in very similar solid-to-gas surface density ratios in this part of the disc. Around the second planetary gap-edge, the midplane model is not quite as efficient at forming planetesimals as the nominal model, and it leaves a ring of dust and pebbles behind that does not exist in the nominal model. At the innermost planetary gap-edge, there is no planetesimal formation at all in the midplane model, and thus the solid-to-gas surface density ratios in the two models are very different.

Fig. 8 2D symmetric disc images for simulations where the maximum particle size was limited to 100 μm. For the sake of easy comparison purposes, images of the same simulations but without the maximum limit on the particle size are added on top of the images at 1 Myr. With smaller particle sizes, it becomes harder to trigger planetesimal formation; the result is that more dust and pebbles remain in the rings at the end of the simulation.

5.1 Imposing a maximum particle size of 100 μm in simulations

We study the effect of having a maximum grain size of 100 μm in our simulations by artificially imposing the maximum pebble size to be 100 μm in the coagulation part of our code. In the bottom panel of Fig. B.1, we show the resulting particle size distribution across the disc for the nominal simulation. The solid-to-gas ratios for the following simulations are shown in Fig. 8: nominal, 0.75 M_iso, 1 M_iso, lowVisc, and lowTurb, with 100 μm particles.

In the nominal simulation, 218 M_⊕ of planetesimals were formed. The amount of dust and pebbles left in the rings is larger compared to the nominal simulation with millimeter-sized particles (see Table 4). Regarding the distribution of pebbles, the Stokes number at the location of the outermost planet is still large enough to prevent most pebbles from drifting across the gap. At the locations of the two innermost planets, this is no longer true. However, because of the linear scaling of the streaming instability with the pressure gradient, pebbles at the gap edges are turned into planetesimals before they have time to drift across the gaps. Therefore, the region interior to the middle planet becomes depleted of dust and pebbles.

When the planetary mass decreases to 0.75 M_iso, the global solid-to-gas ratio remains at around 1% throughout the simulation (see second row of Fig. 8). There are still some rings and gaps that are visible in the particle distribution, and the amount of planetesimals that formed is now 80 M_⊕. For a planetary mass of 1 M_iso, it takes more time to clear the interplanetary regions of pebbles than in the nominal simulation, but at the end of the simulation, the solid-to-gas density ratio across the disc looks very similar. The total mass in planetesimals for this simulation is 232 M_⊕.

A lower viscosity (simulation lowVisc) results in faster drift in the viscously expanding part of the disc and more planetesimal formation at the outer gap edge – in total 237 M_⊕ of planetesimals were formed in this simulation. There is no planetesimal formation at the innermost gap edge in this simulation, and instead there are ten Earth masses of pebbles trapped in that ring, corresponding to around five Earth masses per astronomical unit. The amount of pebbles left in the ring outside the middle planet has also increased compared to the nominal model.

When the turbulent diffusion was lowered by an order of magnitude to 10⁻⁴ (simulation lowTurb), the amount of planetesimals that formed decreased to 192 M_⊕. Lower collisional velocities result in particles growing to 100 μm further out in the disc. These particles obtain higher drift velocities and reach the gap edge of the outermost planet within a million years, causing the wide and bright ring seen in the bottom row of Fig. 8.

5.2 The cases with no pressure scaling and no planetesimal formation

Next we study how the distribution of dust and pebbles change when (1) the dependency of the streaming instability on the pressure gradient is removed, and (2) when planetesimal formation is removed completely (analogous to Sect. 4.2 but for the case with a maximum grain size of 100 μm). The results are presented in Fig. 9.

There is little or no transport through the outermost planetary gap in all simulations; however, the 100 μm sized pebbles do drift past the two innermost gaps in the simulation without planetesimal formation (simulation noPlanetesimal). In the simulation without pressure scaling (simulation noPscaling), some of the pebbles that would otherwise have made it past the gaps are now converted into planetesimals instead, resulting in a quicker depletion of the region interior to the middle planet. Apart from this, simulation noPlanetesimal and simulation noPscaling result in relatively similar images, with the major difference being the brightness and width of the rings. When pressure scaling is added (simulation nominal), the picture changes quite a bit. Efficient planetesimal formation now prevents most pebbles from crossing the middle planetary gap, resulting in that the region interior of this becomes void of pebbles.

6 Comparison to observations

6.1 Dust mass estimates in rings

The amount of dust and pebbles remaining in the rings after 1 Myr varies a lot in our simulations (see Table 4). For example, the amount of dust and pebbles remaining in the outermost ring ranges from 1.7 to 13 Earth masses for simulations in the parameter study. We compare these amounts to dust mass estimates by Dullemond et al. (2018) for rings in the DSHARP survey.

Dullemond et al. (2018) find that the amount of dust stored in each ring is of the order tens of Earth masses. For example, AS 209 was estimated to have around 30 Earth masses of dust trapped in the ring at 69−79 au and 70 Earth masses trapped in the ring at 115−125 au. In general the amount of dust stored in the rings ranges from 1 to 10 Earth masses per astronomical unit (see Table 2 in Dullemond et al. 2018). It should be mentioned that there is much uncertainty as to these estimates, mainly due to the uncertainty in the calculation of the dust opacities.

Comparing with our simulations, the only cases where we have more than ten Earth masses of pebbles remaining in one or several narrow rings after 1 Myr are when we either (1) ignore planetesimal formation completely; (2) ignore the pressure scaling; or (3) use a maximum pebble size of 100 μm (we note that the midplane model is discussed separately below). From Table 4, we find that we have between 2 and 11 Earth masses of dust and pebbles per astronomical unit left in the rings when planetesimal formation is neglected (simulation noPlanetesimal). When the pressure scaling is removed (simulation noPscaling), this value decreases to 0–0.65 Earth masses per astronomical units. Another example is simulation lowVisc for 100 μm-sized particles, where we find that five Earth masses of dust and pebbles per astronomical unit is left in the inner ring, and 0.5 Earth masses per astronomical unit is left in the middle ring.

In order to match the dust mass estimations in the DSHARP rings, we would thus either need a very low planetesimal formation efficiency or some mechanism for destroying the planetesimals and thus replenishing the dust population in the rings (this is discussed in Sect. 7.2). One mechanism, which would result in a higher dust population and likely lead to less planetesimal formation, is efficient fragmentation in the pressure bumps; see Sect. 7.1 for a discussion on this.

The dust masses quoted in Table 4 are for simulations in which we used the planetesimal formation criteria from Yang et al. (2017). In Sect. 4.4 where we compare this criteria to the midplane model, it is shown that the midplane model produces less planetesimals and results in a larger amount of dust and pebbles remaining in the rings. More precisely, after 300 000 yr, the amount of dust and pebbles remaining in the ring at the innermost planet is roughly 9 M_⊕ au⁻¹, and the corresponding amount at the middle planet is 3 M_⊕ au⁻¹. For the nominal model after 300 000 yr, the amount of dust and pebbles remaining in the same rings is roughly 0.6−0.7 M_⊕ au⁻¹. Although the values from the midplane model appear to be a better match to the estimates by Dullemond et al. (2018), the midplane criterion for planetesimal formation has not been confirmed by any hydrodynamical simulation that we know of yet. Since a more detailed comparison of the two planetesimal formation models is beyond the scope of this paper, the rest of the paper is only concerned with the simulations done with our nominal model for planetesimal formation.

Fig. 9 2D symmetric disc images of the evolution of the solid-to-gas surface density for three different versions of the nominal model where the maximum grain size was limited to 100 μm : the nominal model (top row), the nominal model with planetesimal formation but no dependency on the pressure gradient (middle row), and the nominal model without planetesimal formation (bottom row). For the sake of easy comparison purposes, images of the same simulations but without the maximum limit on the particle size were added on top of the images at 1 Myr. When the dependency on the pressure gradient is included, efficient planetesimal formation at the gap edges prevents particles from passing through the gaps, resulting in a depletion of the part of the disc that is closer to the star than the middle planet. This does not happen in the other cases, instead several bright rings are seen in the dust-to-gas ratio, and dust is also left in the innermost region of the disc.

6.2 Global dust distribution

Depending on what planet and disc parameters are used, we end up with very different pebble distributions across the disc. For simulations with a maximum grain size of around one millimeter, high drift velocities and little or no dust transportation through the planetary gaps result in the interplanetary regions becoming depleted of pebbles within a few hundred thousand years. This occurs in all simulations except for the ones with planetary masses that are lower than the pebble isolation mass. Combined with efficient planetesimal formation in the pressure bumps, the region interior to the outermost planet becomes devoid of pebbles (see Fig. 5). Such discs with large central cavities resemble transition discs (e.g., Andrews et al. 2011). However, we want to emphasize that we only show the dust-to-gas surface density ratios in this work. We have not looked into how these discs would actually appear in observations of millimeter continuum emission.

If the planetesimal formation efficiency in nature is lower than assumed in our simulations, so that a significant fraction of millimeter-sized pebbles remain in the pressure bumps, then the discs instead evolve a few very narrow and bright rings, which is similar to the structure observed in the protoplanetary disc around AS 209 (Fedele et al. 2018). This can be seen in Figs. 3 and 9 where we present simulations with no planetesimal formation and no dependency on the pressuregradient. Such discs with a high pebble density in the rings could also be created and maintained through cycles of planetesimal formation and planetesimal destruction and/or efficient fragmentation in the rings (see Drazkowska et al. 2019).

Most observed protoplanetary discs with dust rings nevertheless do not appear as AS 209; instead, they have emission coming more evenly from the whole protoplanetary disc (e.g., ALMA Partnership et al. 2015). Such dust distributions could only be obtained in our simulations by using planetary masses that are lower than the pebble isolation mass. Generally, if particles are transported through the planetary gaps efficiently, then the regions between theplanets are continuously replenished as long as there is a large enough repository of solids far out in the disc. In our nominal simulation, the outer disc still holds around 100 Earth masses of solids after 1 Myr. If dust transport through the gaps is efficient, it further reduces the solid-to-gas ratios in the pressure bumps, leading to less planetesimal formation. Possible reasons for why dust transport through planetary gaps could be more efficient than in our simulations are discussed in Sects. 7.2–7.4.

Introducing a maximum grain size of 100 μm results in the interplanetary regions becoming depleted of their pebble mass on a longer time-scale of up to 1 million years. In these simulations, one ring of pebbles remains visible inside the cavity after 1 Myr even for planetary masses that are larger than the pebble isolation mass. As a comparison, the stars in the DSHARP survey have ages between a few hundred thousand to 10 million years (Andrews et al. 2018), so we know that at least some discs must be able to maintain a high pebble density in the rings for a long time.

6.3 Solar System constraints on planetesimal formation

Since the planetesimals in our model form at the edges of planetary gaps, there must have existed an earlier population of planetesimals, which participated in the formation of these gap-opening planets. Those planetesimals should have formed by some other mechanism than the one proposed in our work, for example, through particle pile-ups outside snow lines (Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017). Our planetesimals would thus represent a second generation of planetesimals, which only form once the first gap-opening planet has appeared in the disc, a scenario that fits in well with Solar System observations.

In Kruijer et al. (2017), they used isotope measurements of iron meteorites together with thermal modeling of bodies internally heated by ²⁶Al decay to study the formation of planetesimals in the asteroid belt. From this study, they conclude the following about the parent bodies of noncarbonaceous (NC) and carbonaceous (CC) iron meteorite: (1) they accreted at different times, within 0.4 respective 0.9 Myr after Solar System formation; (2) they accreted at different locations in the disc, with the CC meteorites accreting further out; and (3) they must have remained separated from before 0.9 until 3− 4 Myr after Solar System formation. This picture could be neatly explained by the formation of the CC iron meteorite parent bodies at the edge of Jupiter’s gap.

The first population of planetesimals form early before 0.4 Myr (the NC iron meteorites). Then at some time before 0.9 Myr, Jupiter reaches the pebble isolation mass and shuts off the flow of pebbles to the inner Solar System. The pressure maximum generated at the edge of Jupiter’s gap promotes planetesimal formation and results in a second generation of planetesimals (the CC iron meteorites). Once Jupiter has reached the pebble isolation mass, it continues to grow slowly for a few million years, keeping the NC and CC populations separate. Then at 3− 4 Myr after Solar System formation, something occurs that scatters the population of CC iron meteorites toward the inner Solar System and causes them to mix with the NC population. This could be the onset of runaway gas accretion (Kruijer et al. 2017) or interactions with an outer giant planet (Ronnet et al. 2018).

7 Shortcomings of the model

7.1 A proper handle on fragmentation

We used a collision algorithm where target particles are reduced to the mass of the projectiles in the event of a destructive collision. This is not fully realistic because destructive collisions should result in the formation of multiple fragments, and it prevents us from recovering very small particle sizes. Blum & Münch (1993) showed that when two similar-sized dust particles collide at a velocity higher than the fragmentation threshold, both particles are disrupted into a power-law size distribution. A significant fraction of the mass in such a collision becomes concentrated in the largest particle sizes (Birnstiel et al. 2011; Bukhari Syed et al. 2017). Since most mass is supposed to be tied up in the largest particle sizes, our simplified algorithm is still justified, but it prevents the formation of small dust particles that could make it past the planetary gap. These particles could recoagulate interior to the gap and result in a population of millimeter-sized particles in the interplanetary region.

Drazkowska et al. (2019) used an advanced 2D coagulation model to study the dust evolution in a disc that is being perturbed by a Jupiter-mass planet. They find that fragmentation at the gap edge is indeed important. In their simulations, large grains that are trapped inside the pressure bump fragment and replenish the population of dust, which can then pass through the planetary gaps. This process thus leads to a continuous dust flux through the gaps. Since this process removes solids from the pressure bump, it should also result in less planetesimal formation; however, a comparison between the timescales for drift, fragmentation, and planetesimal formation would be required in order to further assess this.

7.2 Destruction of planetesimals

In our simulations, we only looked at where planetesimals form and how much mass is turned into planetesimals. This means that we did not investigate what happens to the planetesimals once they have formed. The processes which are likely to be important in determining the fate of planetesimals at the edge of planetary gaps are: dynamical interactions with the planet, dynamical interactions with other planetesimals, and sublimation and erosion due to the flow of gas. These processes will be the subject of a follow-up study.

If the planetesimals are not removed from the gap edge directly after formation, then the density of planetesimals in this region should become very high. In such regions planetesimal collisions are likely to be frequent, and as in debris discs such events result in the production of dust (Wyatt 2008). Another process that could result in a replenishing of the dust population in the pressure bumps is planetesimal sublimation due to bow shocks (Tanaka et al. 2013). The heating and sublimation of the planetesimals result in a shrinking of the planetesimal size and the vapor can form dust particles through recondensation. How efficient and relevant these processes are for the production of dust in a pressure bump remains to be studied.

7.3 The streaming instability may not be 100% efficient

We assume in our model that whenever the critical density to trigger the streaming instability is reached, planetesimals form. The planetesimal formation algorithm used in this study results in a maximum efficiency for planetesimal formation. A calculation of the actual formation efficiency would require taking into account the timescale for collapse into planetesimals, the timescale for particles to drift across the gap, turbulence, and many other effects. However, in simulations where there is no transport of pebbles across the gaps, the efficiency should not play a big role. It does not matter if it takes a hundred years or a hundred thousand years to form planetesimals since the particles remain trapped anyways. In simulations where pebbles are able to make it past the gaps, such as in the simulations with 100 micron-sized particles, the efficiency for planetesimal formation becomes much more important.

One mechanism, which would likely result in less planetesimal formation, is efficient fragmentation in the pressure bump, which is discussed in Sect. 7.1. We also stress that the linear scaling with the pressure gradient is an approximation, and if this relation was less steep or leveled out toward a pressure gradient of zero, more pebbles and dust would be left in the pressure bumps. Furthermore, in 1D simulations we do not need to worry about instabilities at the gap edges. However, if the gaps are deep enough in 2D or 3D simulations, the gap edges may become unstable to form vortices (Hallam & Paardekooper 2017). The triggering of vortices could potentially change the efficiency of planetesimal formation; however, planetesimal formation in such an environment is still poorly understood. Finally, the coagulation model from Güttler et al. (2010) results in a bimodal particle size distribution. In Krapp et al. (2019), it is shown that the streaming instability becomes less efficient when multiple particle sizes are involved. However, the difference in the growth rate between single particles and multiple species appears to vanish when the dust-to-gas ratio is above unity. Therefore, we continued to use the mass averaged Stokes number in our planetesimal formation model, and we did not lower the efficiency when the particle size distribution evolved into bimodal.

7.4 Dust filtration through planetary gaps in 1D versus 2D simulations

There could be a systematic difference in the dust filtration by a planet in 2D simulations relative to our 1D simulations. This is shown by Weber et al. (2018) and Haugbølle et al. (2019) who performed detailed studies of dust filtration through planetary gaps. In these works, they used a dust fluid approach in order to track extremely low values of the dust density. They found that dust is morelikely to be transported through the gaps in 2D simulations, although the amount on the interior of the planet orbit is diminished greatly by the filtering. In Drazkowska et al. (2019), they instead found the opposite results when comparing dust filtration in 1D and 2D coagulation simulations. One possible reason for this discrepancy could be that the gap profiles in Drazkowska et al. (2019) were the same in both the 1D and 2D simulations, while in Weber et al. (2018) the density profiles varied between the 1D and 2D simulations. Regardless, it seems clear that dust filtration in 1D simulations does differ from the 2D or 3D case; however, exactly how is still not certain.

7.5 The α-disc model

In the classical α-disc model, a macroscopic viscosity is assumed to drive angular momentum transport throughout the disc. However, the actual origin of this viscosity is not known. Alternatively, the angular momentum may be drained from the protoplanetary disc by strong winds. In such models, mass is primarily removed from the disc surface and not from the midplane. The surface density profile of such wind-driven discs vary a lot in the literature, and while some resemble the classical α-disc model, others have positive density gradients in the inner regions of the disc and multiple density maxima spread across the disc (Gressel et al. 2015; Bai et al. 2016; Suzuki et al. 2016; Béthune et al. 2017; Hu et al. 2019). In such discs, particle drift would be very different from what we use in our model, and our results would therefore change. However we still do not know enough about what drives angular momentum transport in discs, or what the level of turbulent viscosity and wind transport are, to say anything conclusive about which model is correct. Therefore, here, we decided to stick to the well-understood α-model.

8 Conclusions and future studies

In this work we test the hypothesis that dust trapping at the edges of planetary gaps can lead to planetesimal formation via the streaming instability. To study this, we performed 1D global simulations of dust evolution and planetesimal formation in a protoplanetary disc that is being perturbed by multiple planets. We performed a parameter study to investigate how different particle sizes, disc parameters, and planetary masses affect the efficiency of planetesimal formation. We further compare the simulated pebbles’ distribution with protoplanetary disc observations.

The answers we have obtained for the questions posed in the introduction can be summarized as follows:

• 1.

  Do planetesimals form at the edges of planetary gaps? Planetesimal formation occurs in all of our simulations and is almost exclusively limited to the edges of planetary gaps.

• 2.

  How efficient is this process and how does the efficiency vary with different disc and planet parameters? Planets with masses above the pebble isolation mass trap pebbles efficiently, and in the case of millimeter-sized particles essentially all of these trapped pebbles are converted into planetesimals. As long as the pebbles cannot pass through the gaps, the amount of planetesimals that form does not vary between the simulations, although the region in which they form do change a bit. Decreasing the pebble size to 100 micron results in less efficient conversion of pebbles to planetesimals and more transport through the gaps.

• 3.

  What does the distribution of dust and pebbles look like for the different simulations? In the case of millimeter-sized pebbles and planetary masses that are larger than the pebble isolation mass, the region’s interior to the outermost planet gets depleted of pebbles in a few hundred thousand years. For planetary masses lower than this, transport through the gaps leads to a constant replenishment of the interplanetary region, resulting in a gap-and-ring like pebble distribution. In the case where the particle size was lowered to 100 μm, there is always at least one ring of pebbles remaining inside the cavity. When we lower the efficiency of planetesimal formation, by ignoring the drop in the metallicity threshold for planetesimal formation with decreasing pressure support, the discs instead appear to have narrow and bright rings.

• 4.

  How do these distributions compare with observations of protoplanetary discs? Transition discs with large central cavities are known from observations. Similar discs with large cavities are the general outcome of simulations with massive planets, millimeter-sized pebbles, and efficient planetesimal formation. Discs with narrow and bright rings, which are similar to the outer regions of the disc around the young star AS 209, are the outcome of simulations with massive planets but low planetesimal formation efficiency in the rings. A replenishment of the dust population in the rings, through processes such as fragmentation, planetesimal collisions, or planetesimal evaporation and erosion, could result in similar structures and potentially also aid in transporting particles across the gaps (Drazkowska et al. 2019). Setting the maximum grain size to 100 μm results in multiple rings, longer drift time-scales, and a larger variety in disc structures. Generally, the only simulations that could produce images similar to HL Tau, with multiple gap and ring structures but no strong pebble depletion anywhere in the disc, are simulations with planetary masses that are lower than the pebble isolation mass.

In this work, we have focused on studying the efficiency of planetesimal formation and the locations of the formed planetesimal belts, but we have neglected their further evolution after formation. The processes that are likely to be important in determining the fate of planetesimals that formed at the edges of planetary gaps are as follows: dynamical interactions with the gap-forming planets, planetesimal-planetesimal interactions, and erosion and evaporation by the flow of gas. These processes will be the subject of a follow-up study.

Acknowledgements

The authors wish to thank Alexander J. Mustill and the anonymous referee for helpful comments that led to an improved manuscript. The authors further wish to thank Alessandro Morbidelli for inspiring discussions, and Chao-Chin Yang for comments and advice regarding the planetesimal formation model. L.E., A.J. and B.L. are supported by the European Research Council (ERC Consolidator Grant 724 687-PLANETESYS). A.J. further thanks the Knut and Alice Wallenberg Foundation (Wallenberg Academy Fellow Grant 2017.0287) and the Swedish Research Council (Project Grant 2018-04867) for research support. The computations were performed on resources provided by the Swedish Infrastructure for Computing (SNIC) at the LUNARC-Centre in Lund, and are partially funded by the Royal Physiographic Society of Lund through grants.

Appendix A Particle collisions

Particle collisions were performed through a Monte Carlo method. Each particle swarm was assigned a total mass M_i, an individual particle mass m_i, and a number density n_i = M_i∕m_i. The total mass M_i is different for different particles, which is useful in order to resolve a wide range of column densities in the disc.

When two particles collide, we define the larger particle as the target and the smaller particle as the projectile. The rate of interaction (defined below) between the target and the projectile is determined as $$r_{\text{tp}}={\sigma }_{\text{tp}}{v}_{\text{tp}}{n}_{\text{p}}\frac{m_{\text{p}}}{m_{\text{t}}}\max (M_{\text{t}}/M_{\text{p}},1),$$(A.1)

where σ_tp and v_tp are the collisional cross section and relative speed between the target (t) and the projectile (p). If the total mass in the projectile is larger than or equal to the total mass in the target, then the interaction time-scale is defined as the time for each target particle in a swarm to collide with its own mass in projectiles. If the total mass in target particles is larger than the total mass in projectiles, we multiplied it by M_t∕M_p so that the interaction time-scale is instead the time-scale for all projectile particles to have collided with a target particle.

Equation (A.1) can be rewritten as $$\begin{array}{ll}{r}_{\text{tp}}\hfill & ={\sigma }_{\text{tp}}{v}_{\text{tp}}{n}_{\text{t}}\frac{M_{\text{p}}}{M_{\text{t}}}\max (M_{\text{t}}/M_{\text{p}},1)\hfill \\ \hfill & ={\sigma }_{\text{tp}}{v}_{\text{tp}}{n}_{\text{t}}\max (M_{\text{p}}/M_{\text{t}},1).\hfill \end{array}$$(A.2)

Using that the collisional cross section is ${\sigma }_{\text{tp}}=\pi {(s_{\text{t}}+s_{\text{p}})}^2$, we obtain the final equation for the mass doubling rate $$r_{\text{tp}}=\pi {(s_{\text{t}}+s_{\text{p}})}^2{v}_{\text{jk}}{n}_{\text{t}}\max (M_{\text{p}}/M_{\text{t}},1),$$(A.3)

where s is particle radius. The relative speed contains contributions from Brownian motion, differential radial and azimuthal drift, differential reaction to the gas accretion speed, and turbulent speed. The turbulent speed is based on the closed-form expressions of Ormel & Cuzzi (2007), all of the other terms are standard in the literature (see e.g., Brauer et al. 2008).

In order to average over the vertical direction, we assume that the two particle species maintain a Gaussian density profile in the vertical direction and that changes to the target particle by coagulation are immediately diffused over the entire column density of the target particle. Here, we follow a similar approach as Brauer et al. (2008), their Appendix B. The coagulation equation written for a given height z over the midplane is $$r_{\text{tp}}(z)=\pi {(s_{\text{t}}+s_{\text{p}})}^2{v}_{\text{tp}}{n}_{\text{p}}(z).$$(A.4)

The collision rate averaged over n_t(z) is then $$\overline{r_{\text{tp}}}=\frac{{\displaystyle \int r_{\text{tp}}}(z)n_{\text{t}}(z)\text{d}z}{{\displaystyle \int n_{\text{t}}}(z)\text{d}z}.$$(A.5)

Assuming a Gaussian density distribution with a midplane number density n_0,t and scale-heights H_t and H_p, respectively, the integration yields $$\overline{r_{\text{tp}}}=\pi {(s_{\text{t}}+s_{\text{p}})}^2{v}_{\text{tp}}{n}_0,t\times \frac{1}{\sqrt{1+{(H_{\text{t}}/H_{\text{p}})}^2}}.$$(A.6)

In comparing this to Eq. (A.4), we see that the vertical integration of the coagulation equation can be treated as a simple multiplication factor on the rate of collisions in the midplane.

In order to calculate the number density of particles in the midplane, we divided the mass in each superparticle by the area of the annulus where the particle is present (2πrΔr) and then by $1/[\sqrt{2\pi }{H}_{\text{p}}]$ to obtain the midplane density. The time-step contribution from particle coagulation is based on the interaction rate r_ij. The time-step for a particle i is $${\tau }_i=\frac{1}{{\sum }_j(r_{ij})}.$$(A.7)

The Monte Carlo time-step is then ${\min }_i({\tau }_i)$ times a numerical factor, which was chosen to be 0.2. Once the time-step was calculated, we looped over all of the particles in a grid cell and all their unique partners. For each particle pair, we drew a random number, and if that number was smaller than d t × r_ij, we let the swarms interact. We based the outcome of collisions on experimental results by Güttler et al. (2010), and we assume that the particles are porous. The possible outcomes of a collision are sticking, bouncing, bouncing with mass transfer, and fragmentation. Sticking means that the target either doubles its mass or multiplies its mass by (1 + M_p∕M_t), if M_p < M_t. For bouncing with mass transfer, we doubled the mass of the projectile particles and subtracted the projectile particle mass from each target particle. For fragmentation, we set all the target and projectile particles to the mass of the projectile. If there are excess target particles, then they retain their original mass.

Appendix B Particle size distributions

Figure B.1 shows the size distributions of particles at different times during disc evolution for some selected simulations. When a planetary mass is used that is lower than one pebble isolation mass (top panel), the gaps are never completely depleted of dust and pebbles. The pile-up of material at the gap edges is also much less prominent, and since pebbles are now transported through the planetary gaps, the result is a more even distribution of particles throughout the disc. Since there are more particles in the interplanetary regions, we also get more spontaneous concentrations, leading to more planetesimal formation. A comparison with Fig. 4 shows that the amount of planetesimals forming in the interplanetary regions is still negligible compared to the amount that form at the gap edges.

Fig. B.1 Particle size distributions at different times during disc evolution for four different simulations. The semimajor axes of the formed planetesimals are indicated at the top of the plots. Top panel: in the simulation with a planetary mass of 0.75 Miso, planetesimals form in a wider region around the pressure bump than in the nominal simulation. Less efficient pebble trapping at the pressure bump also results in a more even distribution of dust and pebbles in the disc, with no strong depletion at the location of the planets. Second panel: lowering the turbulence diffusion to 10−4 results in slower collisional velocities, which results in slower coagulation, but eventually leads to larger particle sizes. Third panel: when the viscous parameter is lowered to 10−4 as well, we get small bumps in the gas surface density profile at the inner edges of the planetary gaps. Particles become trapped in these bumps, which also result in some planetesimal formation at these locations. Bottom panel: this plot shows the implementation of a maximum grain size of 100 μm.

When the turbulent diffusion is lowered by an order of magnitude to 10⁻⁴ (second panel), the coagulation time-scale increases. Since it takes more time for particles to grow, it also takes more time for the size distribution to become bimodal. The slow particle growth also results in the drift time-scale being longer initially. However, decreasing the amount of turbulence results in lower collisional speeds, which in turn results in larger particle sizes. For such large particles, the time-scale for drift becomes shorter than in the nominal model. The result is that more particles make it from the exponentially tapered outer disc to the inner 100 au where the planets reside. The larger particle sizes also lead to sporadic concentrations becoming more common, which again result in more planetesimal formation in the interplanetary regions.

In the simulation where both the viscous parameter and the turbulence diffusion were decreased to 10⁻⁴ (second panel), trapping at the inner gap edges results in a significant amount of planetesimals being formed at these locations. A small bump in the gas surface density profile is created at the beginning of all simulations, but when the viscous parameter is high, this bump disappears before planetesimal formation is initiated. For a viscous parameter of 10⁻⁴, this pile-up of gas at the inner gap edges is both more prominent and longer lasting than in all other simulations. Since the viscosity is small, the time for gap-clearing is also longer.

In the bottom panel of Fig. B.1, we show the size distribution for particles in the nominal model when a maximum grain size of 100 μm was applied.This constraint does not matter for the particle evolution far out in the disc since particles do not grow that large anyways.

References

All Tables

All Figures

Fig. 1 Top panel: evolution of the gas surface density, which is normalized against the initial gas surface density, across the protoplanetary disc for the nominal model with planets (solid lines) and without planets (dotted lines). The vertical dashed linesmark the semimajor axes of the planets, and they coincide with the locations of the three major gaps in the disc around HL Tau (Kanagawa et al. 2016). Bottom panel: evolution of the solid-to-gas surface density ratio across the protoplanetary disc for the nominal model. The solid component is divided into planetesimals (marked with gray) and dust+pebbles. Planetesimals form in narrow rings at the location of the gap edges and inside the planetary gaps (the amount of planetesimals formed inside the gaps is negligible). The interplanetary regions are depleted of dust and pebbles within a few hundred thousand years.

In the text

Fig. 2 Top panel: size distribution of particles in the protoplanetary disc at different times during disc evolution for our nominal model. The semimajor axes of the formed planetesimals are indicated at the top of the plot. Particles were initialized with an equal spacing all over the disc and with a radius of 1 μm. Efficient coagulation and high drift velocities in the inner disc result in a depletion of the interplanetary regions after only a few 100 000 yr. These processes occur on a longer timescale in the outer disc. Bottom panel: particle Stokes number versus semimajor axis for the same data as in the top panel.

In the text

Fig. 3 2D symmetric disc images of the evolution of the solid-to-gas surface density (excluding the formed planetesimals) for three different versions of the nominal model: the nominal model (top row), the nominal model with planetesimal formation but with no dependency on the pressure gradient (middle row), and the nominal model without planetesimal formation (bottom row). In the nominal model, essentially everything that enters the pressure bump is converted into planetesimals, resulting in a large cavity in the distribution of dust and pebbles. When the pressure dependence is neglected, the critical density required to trigger the streaming instability increases, and we see some rings in the dust and pebble distribution. When planetesimal formation is removed completely, we are left with three rings in which the dust and pebble density is very high.

In the text

Fig. 4 Histogram showing the total amount of mass in planetesimals that has formed at different locations in the disc after 1 Myr for each simulation in the parameter study. The dotted vertical lines indicate the semimajor axes of the planets. When the pressurescaling is removed (simulation “noPscaling”), the critical density required for forming planetesimals at the gap edge increases, resulting in less planetesimal formation. The amount of planetesimal formation does not vary a lot between the simulations where the planetary mass is larger than the pebble isolation mass since millimeter pebbles cannot drift past the gaps and planetesimal formation is extremely efficient where the pressure gradient is close to zero. For planetary masses that are lower than the pebble isolation mass, the amount of planetesimals that form rapidly decreases with decreasing planetary mass. Lowering the viscosity parameter and the turbulent diffusion results in faster drift velocities in the viscously expanding part of the disc. Because of this, more pebbles reach the outermost gap edge and are turned into planetesimals. Increasing the metallicity (simulation highMetal) results in sporadic planetesimal formation in the interplanetary regions. When a constant migration speed is added to the planets, the location of the gap edges are shifted inward with time, resulting in the region where planetesimals form to shift inward accordingly.

In the text

Fig. 5 2D symmetric disc images of the evolution of the solid-to-gas surface density excluding planetesimals for all simulations in the parameter study. In the nominal simulation where we used a planetary mass of two times the pebble isolation mass fast pebble drift, little or no transport through the planetary gaps and efficient planetesimal formation at the gap edges result in that the part of the disc that is closer to the star than the outermost planet gets depleted of dust and pebbles. In the simulations where the planetary mass is lower than the pebble isolation mass (simulations 0.50 Miso and 0.75 Miso), increased transport past the gaps and less planetesimal formation results in a gap-and-ring-like structure. The width of the region where dust and pebbles are trapped depends on the planetary mass. Therefore, when theplanetary mass is changed to one or three times the pebble isolation mass (simulations 1 Miso and 3 Miso), and we compare to the nominal model, the width of the rings becomes narrower or larger, respectively. Except for this, the results are the same as in the nominal model. In simulation lowVisc, the viscosity parameter is lowered, resulting in a slower clearing of the gaps and interplanetary regions, as well as faster radial drift in the viscously expanding part of the disc. When the turbulent diffusion is lowered (simulation lowTurb), the collisional velocities decrease, resulting in larger particles which drift faster toward the star. This causes the bright ring that can be seen beyond the outermost planet at the end of the simulation. In simulation lowViscTurb, the turbulent diffusion is kept at the same level as in simulation lowTurb, but the viscosity parameter is lowered by an extra order of magnitude compared to simulation lowVisc. The combinationresults in that essentially all solids in the outer disc reach the outermost planetary gap before the end of the simulation, causing the narrow bright ring seen in the dust-to-gas ratio. In simulation highMetal, the initial solid-to-gas ratio in the disc is increased to 2%, which is a change that does not have a big effect on the appearance of the disc at the end of the simulation. Finally in the last simulation, the planets were given a constant radial velocity directed toward the star (simulation migration), resulting in a smaller radius of the cavity.

In the text

Fig. 6 Histogram showing the total amount of mass in planetesimals that have formed at different locations in the disc after 300 kyr, for the nominal model and the midplane model. The amount of planetesimals formed around the gap edges of the two outermost planets are similar in both models; however, this is not the case at the innermost planetary gap edge. The settling toward the midplane is not efficient enough to counteract the stirring by turbulence in this region, and thus a dust-to-gas density ratio of unity in the midplane is never reached.

In the text

Fig. 7 Evolution of the solid-to-gas surface density ratio across the protoplanetary disc for the nominal model (top panel) and the midplane model (bottom panel). The solid component includes dust and pebbles, whereas planetesimals are excluded. The two models produce very similar amounts of planetesimals around the outermost planetary gap-edge, resulting in very similar solid-to-gas surface density ratios in this part of the disc. Around the second planetary gap-edge, the midplane model is not quite as efficient at forming planetesimals as the nominal model, and it leaves a ring of dust and pebbles behind that does not exist in the nominal model. At the innermost planetary gap-edge, there is no planetesimal formation at all in the midplane model, and thus the solid-to-gas surface density ratios in the two models are very different.

In the text

Fig. 8 2D symmetric disc images for simulations where the maximum particle size was limited to 100 μm. For the sake of easy comparison purposes, images of the same simulations but without the maximum limit on the particle size are added on top of the images at 1 Myr. With smaller particle sizes, it becomes harder to trigger planetesimal formation; the result is that more dust and pebbles remain in the rings at the end of the simulation.

In the text

Fig. 9 2D symmetric disc images of the evolution of the solid-to-gas surface density for three different versions of the nominal model where the maximum grain size was limited to 100 μm : the nominal model (top row), the nominal model with planetesimal formation but no dependency on the pressure gradient (middle row), and the nominal model without planetesimal formation (bottom row). For the sake of easy comparison purposes, images of the same simulations but without the maximum limit on the particle size were added on top of the images at 1 Myr. When the dependency on the pressure gradient is included, efficient planetesimal formation at the gap edges prevents particles from passing through the gaps, resulting in a depletion of the part of the disc that is closer to the star than the middle planet. This does not happen in the other cases, instead several bright rings are seen in the dust-to-gas ratio, and dust is also left in the innermost region of the disc.

In the text

Fig. B.1 Particle size distributions at different times during disc evolution for four different simulations. The semimajor axes of the formed planetesimals are indicated at the top of the plots. Top panel: in the simulation with a planetary mass of 0.75 Miso, planetesimals form in a wider region around the pressure bump than in the nominal simulation. Less efficient pebble trapping at the pressure bump also results in a more even distribution of dust and pebbles in the disc, with no strong depletion at the location of the planets. Second panel: lowering the turbulence diffusion to 10−4 results in slower collisional velocities, which results in slower coagulation, but eventually leads to larger particle sizes. Third panel: when the viscous parameter is lowered to 10−4 as well, we get small bumps in the gas surface density profile at the inner edges of the planetary gaps. Particles become trapped in these bumps, which also result in some planetesimal formation at these locations. Bottom panel: this plot shows the implementation of a maximum grain size of 100 μm.

In the text