© ESO 2018

1 Introduction

A theory of planet formation should be able to explain the variety of planetary systems discovered within a coherent framework. In particular, the presence of gas giants poses a fundamental constraint. Within the core accretion scenario, the interstellar dust grains need to grow from μm size to a 10–20 M_⊕ planetary core before gas accretion sets in. This growth must happen before the star can photoevaporate the gas disc, which occurs on a timescale of ~3 Myr (Hillenbrand 2008). Moreover, during the growth period, planets are embedded in the ambient disc and their orbital evolutions are determined by planet–disc and planet–planet interactions. Some planets may end up accreted onto the star or ejected from the system if no other physical mechanisms intervene to stop them (see e.g. Alexander & Pascucci 2012; Ercolano & Rosotti 2015). Planet formation and evolution are determined by the structure of the protoplanetary discs in which they form. The observations of these discs can give some information about their masses, rotation, and density profile (Williams & Cieza 2011). The observed diversity in the sample of extrasolar planets indicates that the evolution of a planet may depend on variations in the initial conditions or random (external or internal) events occurring during this crucial phase. An important initial condition is the stellar environment of the growing planetary system, which can strongly affect the disc lifetime by tidally truncating the outer regions of its birth disc and the dynamical evolution of the planetary system (see e.g. Picogna & Marzari 2015).

We can place some constraints on the initial conditions and the giant planet formation models by studying their current physical and chemical properties. The Galileo mission measured the abundances of various elements in the outer layers of Jupiter. Young (2003) found that they were in the range of 2–4 times solar, with a predicted core mass in the range 0–18 M_⊕, strongly dependent on the assumed equation of state (Fortney & Nettelmann 2010). The internal composition has also been derived for hot Jupiters such as HAT-P-13 b, where Buhler et al. (2016) used the analysis of secondary eclipses of the planet to infer a core mass of M_c < 25M_⊕ with a most likely value of 11 M_⊕. These observations can be explained by a bottom-up model of planet formation, such as the core accretion model (Pollack et al. 1996), which predicts an enriched solid composition respect to the solar one. Within this framework, we are interested in studying the process that can explain how the minimum solid core mass, necessary to rapidly build up a massive gaseous envelope, can be accreted within the disc lifetime.

The solid materials accreting onto the forming planetary core can have different origins based on the local size distribution (and Stokes number) of the solid disc. One solid reservoir consists of the planetesimals that are gravitationally perturbed by the planetary embryo. If they can cross the mean motion resonances with the planet and enter into its gravitational influence zone, they can end up being accreted onto it. This model of solid core accretion via planetesimals can explain a certain class of gas giants within few au from the central star, but the timescales needed to form the observed planets at tens or hundreds of au are prohibitive.

One possibility to overcome this limitation is the rapid accretion of pebble-sized particles (Ormel & Klahr 2010; Lambrechts & Johansen 2012). In this context pebbles are centimetre-to-meter-sized objects that strongly interact with the gas via drag force. Planetary embryos with an increased gas density in their proximity have an enhanced sphere of influence to accrete solid particles, and they can interact with an higher flux of pebble-sized particles owing to their fast drift speed. Adopting this accretion channel, the timescale of giant planet formation can be lowered, making the build up of gas giants at tens of au possible, if a significant reservoir of pebbles is present. In 2D hydrodynamical simulations with embedded particles, this fast accretion was confirmed (Morbidelli & Nesvorny 2012). The limiting factor for this accretion process is given by the need of a continuous resupply of material from the outer disc because the drift timescale of this pebble-sized particles is very short. The formation of a strong pressure gradient created, for example, by the growing planet can filtrate particles, thereby preventing these particles from reaching the planet or inner parts of the disc; this process sets in at the so-called pebble isolation mass of the core and can explain a class of the observed transition discs. In this paper, we address these limiting factors by studying the evolution of a variety of particle sizes in a global turbulent disc, deriving their accretion rates onto the planet, and obtaining new estimates on the pebble isolation mass.

We consider the evolution of particles in discs made turbulent by a purely hydrodynamical process, the vertical shear instability (VSI) as described in Nelson et al. (2013). Recently, Stoll & Kley (2016) showed that the dust dynamics in VSI-turbulent discs, following the gas behaviour, has a drift speed directed inwards at the disc midplane and outwards in its upper layers in a similar way as for global MHD simulations (Flock et al. 2011). This is exactly opposite to the meridional flow observed for laminar viscous discs. This phenomenon can have important effects on the planet formation process, resupplying materials to the outer disc regions and explaining the observed chondrule population in the outer regions of the solar system (Bockelée-Morvan et al. 2002). Moreover, Stoll & Kley (2016) found that the strong vertical motions induced by the VSI were able to collect pebble-sized particles in rings of high surfacedensity and low relative velocity, potentially aiding the planet formation process through streaming instability (Youdin & Goodman 2005; Auffinger & Laibe 2018). In this work, we extend our recent analysis of the planet–disc interaction in a laminar and turbulent disc (Stoll et al. 2017b) by adding dust particles to the simulations and study their accretion dynamics on the planet. We perform two sets of models. In the first set, the dust is embedded in VSI turbulent discs (with no explicit viscosity added) and in the second series corresponding viscous discs models are performed. This allows us to disentangle the effect of turbulence on the planet–dust interaction and the resulting accretion rate of solid particles. Recently, Xu et al. (2017) studied the accretion of pebbles on small cores in turbulent discs driven by the magneto-rotational instability (MRI). We compare our results to their study.

In Sect. 2 we describe the various forces acting on dust particles embedded in the protoplanetary disc, and then focussing in Sect. 3 on the models for planet–disc interaction and solid particle accretion. In Sect. 4 we describe the adopted set-up for the numerical analysis and the main results obtained in Sect. 5. Finally, we discuss the obtained solid accretion rate onto a growing planetary core in Sect. 6 and draw the main conclusions in Sect. 8.

2 Dust dynamics

We consider a thin vertically isothermal gaseous disc with an embedded protoplanet of mass M_p orbiting around a Sun-like star. Additionally, we follow simultaneously the motion of dust particles of various sizes whose motions are determined by the star, planet, and turbulent gas. In the VSI turbulent disc models, the particles experience the normal drag forces due gas–particle interaction; see Sect. 2.2. On the other hand, in the viscous disc models, the effect of the underlying turbulence is modelled via additional stochastic kicks on the particles as described in Sect. 2.3, in addition to the regular drag forces.

2.1 Equations of motion

A dust particle immersed in the disc is subject to (i) the gravitational force of the central star and the protoplanet, (ii) the drag force due to the varying velocity between the dust orbiting with a Keplerian speed, and the gas, which rotates with a slightly sub-Keplerian speed (Whipple 1964), due to the radial pressure gradient that partially supports it against the stellar gravity, (iii) gas turbulent motion, which influences the dust dynamics by radially and vertically spreading small particles; (iv) photophoretic gas pressure and radiation pressure, which we do not take into account since we focus mainly on a region close to the disc midplane where the efficiency of these processes is expected to be low, (v) growth/fragmentation (depending on their composition/relative speed) due to collisions between grains (see e.g. Testi et al. 2014). We do not consider this because we study the dynamics of particles of varying sizes remaining agnostic about the dust size distribution.

We define a reference frame in spherical coordinates, centred at the location of the star with mass M_⋆ = 1 M_⊙ and co-rotating with constant angular velocity Ω_f, following a protoplanet of mass M_p and fixed position r_p. Then, we can describe the equation of motion of a dust particle of mass m_d as $$m_{\text{d}}{\stackrel{\mathrm{..}}{r}}_{\text{d}}=F_{\text{grav}}+F_{\text{drag}}+F_{\text{turb}}+F_{\text{nonin}}.$$(1)

In this equation, the first term is the gravitational interaction with the star and planet, $$F_{\text{grav}}=-\frac{G{M}_{\star }{m}_{\text{d}}}{|r_{\text{d}}{|}^3}{r}_{\text{d}}+\frac{G{M}_{\text{p}}{m}_{\text{d}}}{|r_{\text{p}}-r_{\text{d}}{|}^3}(r_{\text{p}}-r_{\text{d}}),$$(2)

the second term is the drag force (see Sect. 2.2), the third term is the turbulence force (see Sect. 2.3), and the last is the non-inertial term imparted to the star by the planet, $$F_{\text{nonin}}=-\frac{G{M}_{\text{p}}{m}_{\text{d}}}{|r_{\text{p}}{|}^3}{r}_{\text{p}}.$$(3)

We ignore the self-gravity of the disc since we model a low mass disc in which a large planetary core has already formed.

2.2 Drag force

The drag force acting on a particle depends strongly on the physical condition of the gas and the shape, size, and velocity of the particle. We limit ourselves to spherical particles, for which the drag force always acts in the direction opposite to the relative velocity. The drag regime experienced by a dust particle is described by three non-dimensional parameters as follows:

• 1.

  The Knudsen number, K = λ∕(2s), is the ratio of two characteristic length scales of the system: the mean free path of the gas molecules λ and the particle size, where s denotes the particle radius.

• 2.

  The Mach number, M = v_r∕c_s, is the ratio of the relative velocity between dust and gas, v_r, to the gas sound speed c_s.

• 3.

  The Reynolds number is given by $$Re=\frac{2v_{\text{r}}s}{{\nu }_{\text{m}}},$$(4)

  where ν_m is the gas molecular viscosity defined as $${\nu }_{\text{m}}=\frac{1}{3}\left(\frac{m_0{\overline{v}}_{\text{th}}}{\sigma }\right),$$(5)

  and m₀ and ${\overline{v}}_{\text{th}}=\sqrt{\pi /8}{c}_{\text{s}}$ are the mass and mean thermal velocity of the gas molecules, and σ is their collisional cross section.

2.2.1 Drag law

We adopt a law that can model the drag force for a broad range of Knudsen numbers, using the approach implemented by Woitke &Helling (2003), who used a quadratic interpolation between the Epstein and Stokes regimes $$F_{\text{drag}}={\left(\frac{3K}{3K+1}\right)}^2{F}_{\text{drag},\text{E}}+{\left(\frac{1}{3K+1}\right)}^2{F}_{\text{drag},\text{S}}.$$(6)

For large Knudsen numbers, the first term dominates reducing the drag to the Epstein regime (Baines et al. 1965; Kwok 1975), $$F_{\text{drag},\text{E}}=-\frac{4}{3}\pi {\left(1+\frac{9\pi }{128}{M}^2\right)}^{1/2}{\rho }_{\text{g}}(r_{\text{p}})s^2{\overline{v}}_{\text{t}}{v}_{\text{r}},$$(7)

where ρ_g(r_p) is the gas density at the particle location. For small Knudsen numbers, the second term dominates leading to the Stokes regime $$F_{\text{drag},\text{S}}=-\frac{1}{2}{C}_{\text{D}}\pi s^2{\rho }_{\text{g}}(r_{\text{p}})v_{\text{r}}{v}_{\text{r}},$$(8)

where the drag coefficient C_D for low Mach numbers is (Whipple 1972; Weidenschilling 1977) $$C_{\text{D}}\simeq \{\begin{array}{ll}24R{e}^{-1}\hfill & Re<1\hfill \\ 24R{e}^{-0.6}\hfill & 1<Re<800\hfill \\ 0.44\hfill & Re>800.\hfill \end{array}$$(9)

For more information, see Picogna & Kley (2015, and references therein).

2.2.2 Stopping time

A fundamental parameter to determine the strength of the drag force is the stopping time, t_s, defined as $$F_{\text{drag}}=-\frac{m_{\text{d}}}{t_{\text{s}}}{v}_{\text{r}}.$$(10)

It approximates the timescale on which the embedded gas particles approach the velocity of the gas. In the Epstein regime, the stopping time takes the form $$t_{\text{s}}=\frac{s{\rho }_{\text{s}}}{{\rho }_{\text{g}}{\overline{v}}_{\text{th}}},$$(11)

where ρ_s is the internal particle density. It is also useful to derive a dimensionless stopping time (or, hereafter, Stokes number) as $${\tau }_{\text{s}}=t_{\text{s}}{\Omega }_{\text{K}}(r),$$(12)

where Ω_K is the Keplerian orbital frequency. The Stokes number τ_s (sometimes abbreviated by St) describes the effect of a drag force acting on a particle independent of its location within the disc. With our definition of the stopping time in Eq. (11), the Stokes number is defined in the midplane of the disc.

2.3 Turbulence

Turbulence in the gas acts to stir up well-coupled solid particles, preventing the settling process into a thin layer at the disc midplane. In general, the source of this turbulence is unknown (either driven by MHD or purely hydrodynamic processes as in our case), but it is responsible for both angular momentum and particle transport within the disc (Armitage 2010). By equating the gravitational force in the vertical direction |F_grav|_z with the drag force |F_drag|_z, we can derive a characteristic settling speed $$v_{\text{set}}=t_{\text{s}}{\Omega }_{\text{K}}^{2}z.$$(13)

The condition for which the turbulence strength can counteract the vertical settling of small dust particles is then obtained by comparing the settling time $$t_{\text{set}}=\frac{z}{v_{\text{set}}}=\frac{1}{t_{\text{s}}{\Omega }_{\text{K}}^2},$$(14)

to the time t_diff the turbulence needs to erase the spatial gradients in the particle concentration $$t_{\text{diff}}=\frac{z^2}{D_{\text{d}}},$$(15)

where D_d is the turbulent diffusion coefficient of the particles (dust). If one assumes that D_d equals, to a first approximation, the diffusion coefficient for the gas D_g and that we can write D_g ≃ αc_sh, assuming that the turbulence acts like an effective viscosity (Shakura & Sunyaev 1973), then one can derive the minimum α value required to prevent dust settling at one scale height, z = h $$\alpha \gtrsim {\tau }_{\text{s}}.$$(16)

In the following we use a more complex turbulent diffusion model that distinguishes between D_d and D_g.

Turbulent diffusion model. The source of turbulence in planet-forming discs is unknown. It can depend strongly on the environment, and different sources might be dominant in the various regions and during the evolution of the disc. In the laminar disc simulation, we do not consider the origin of the turbulence and use a simplified turbulence diffusion model to evolve the dust population. The basic idea is to mimic turbulent transport as a diffusive process (through a Brownian motion; Dubrulle et al. 1995; Youdin & Lithwick 2007; Charnoz et al. 2011) with a stochastic term in the equation of dust motion to account for the kicks induced by the turbulent gas velocity field. We model the kick on the particle position as a random Gaussian variable δr_d,T with mean ⟨δr_d,T⟩ and variance ${\sigma }_{\text{d},\text{r}}^2$ depending on the dust diffusion coefficient D_d as follows: $$\delta r_{\text{d},\text{T}}=\{\begin{array}{l}\langle \delta r_{\text{d},\text{T}}\rangle =\frac{D_{\text{d}}}{{\rho }_{\text{g}}}\frac{\partial {\rho }_{\text{g}}}{\partial x}\text{d}t,\hfill \\ {\sigma }_{\text{d},\text{r}}^2=2D_{\text{d}}\text{d}t,\hfill \end{array}$$(17)

where dt is the time step and ∂∕∂x is the spatial derivative along the considered direction. The relation between particle and gas diffusion can be written as $$D_{\text{d}}=\frac{D_{\text{g}}}{Sc},$$(18)

where Sc is the Schmidt number (Youdin & Lithwick 2007): $$Sc=\frac{1+{\Omega }_{\text{K}}^2{\tau }_{\text{s}}^2}{1+4{\Omega }_{\text{K}}{\tau }_{\text{s}}}.$$(19)

This prescription for the turbulent diffusion assumes that the diffusion coefficients in the vertical and radial direction are identical. This is a crude assumption as we have shown in Stoll et al. (2017a) because the α parameter can differ in the radial and vertical directions by more than two orders of magnitude; see also the discussion in Youdin & Lithwick (2007). Nevertheless, we use this simplified description to model the vertical and radial spread of dust particles in the viscous disc simulations and find results comparable to those created self-consistently by VSI turbulence as we show later.

3 Planet–solid disc interaction

3.1 Particle accretion

The idea that gas plays a pivotal role in the accretion of solids by planetary cores was first introduced in the Kyoto model (Hayashi et al. 1977; Nakazawa & Nakagawa 1981; Nakagawa et al. 1983). The basic concept was that the orbital decay experienced by planetesimals is size dependent, occurring at a lower rate for larger bodies. In this way, a large embryo can grow as drag feeds it with dust and small planetesimals. Later, Weidenschilling & Davis (1985) found that orbital resonances with the growing core can effectively filter a significant fraction of planetesimals for which the drag force is not strong enough to allow them to cross those stable regions. However, since eccentricities are pumped up at resonances, collisions between large planetesimals become more frequent, increasing the fraction of smaller bodies that can cross the resonances and accrete onto the planetary core.

On the other hand, Kary et al. (1993) showed that even if a body is small enough to cross all the resonances, it can avoid being accreted. The impact probability typically ranges between 10% and 40% but can be higher if the core possesses an extended atmosphere (D’Angelo et al. 2014). More generally, the accretion rate is inversely proportional to the strength of the drag force and the inclination of the planetesimal. Moreover, Kary et al. (1993) found that for cores with mass ratio q = M_p∕M_⋆ > 10⁻⁵, the material approaching the planet can be captured into a stable orbit around the planet, thereby forming an accretion disc around it.

This strong perturbation in the local environment of the protoplanet creates pressure gradients that impact the evolution of dust and planetesimals (Paardekooper & Mellema 2004; Paardekooper 2007). In particular, small protoplanets, depending on their surface and temperature profiles, can carve a gap in the dust disc even if there is no gas gap (Picogna & Kley 2015; Rosotti et al. 2016; Dipierro et al. 2016).

3.2 Resonances

A particle that migrates within the disc feels a stronger (regular) gravitational interaction with a planetary companion when it reaches specific locations in the disc where its mean motion n_d = 2π∕T, where T is its orbital period, is a multiple of the planet mean motion n_P $$\frac{n_{\text{d}}}{n_{\text{P}}}=\frac{(l+m)}{l},$$(20)

where l and m are integer numbers. These are called mean motion orbital resonances (MMR), where m gives the order of the resonance. These MMRs can effectively excite the eccentricity and inclination of particles, potentially halting their drift process. Their strength grows for decreasing values of m and increasing l. Thus, focussing on first order MMR (m = 1), the larger l, the smaller the particles that can be stopped from accreting onto the planet. The resonances are yet not able to halt all the particles because they become more and more closely spaced as l grows until the point at which they overlap leading to a chaotic behaviour of the dust particles that can cross the higher order resonances (Wisdom 1980). The minimum size s_min of a particle for which the resonant perturbations due to a planet with mass ratio q are stronger than the drag force is (Weidenschilling & Davis 1985) $$s_{\text{min}}=\frac{{\rho }_{\text{g}}h{a}_{\text{d}}}{3{\rho }_{\text{d}}qC(l)l^{3/2}},$$(21)

where C(l) is an increasing function of l, a_d is the particle semi-major axis, and the region of chaotic behaviour close to the planet location starts at (Duncan et al. 1989) $$|r-a_{\text{d}}|\simeq 1.5q^{2/7}.$$(22)

This relation depends strongly on the local gas properties. When the planet opens up a gap in the gaseous disc, reducing the gas surface density, the particle Stokes number increases; thus the inner resonances can halt a larger fraction of incoming particles, which are less coupled to the gas.

For a planet with q = 10⁻⁴, Paardekooper & Mellema (2004) found three visible regimes. Particles with Stokes number less than τ_s ≃ 0.1 are well coupled to the gas, and they always reach the planet surface. On the other hand, particles with τ_s > 10 are trapped in external resonances, and their accretion rate is very low. Finally, the intermediate regime is reaching the co-orbital region of the planet, but not all of them are accreting as predicted by Kary et al. (1993).

4 Set-up

We used the PLUTO code (Mignone et al. 2007) and modified it to take into account the evolution of partially coupled particles. The main parameters of the reference simulations are summarised in Table 1. The simulations analysed in this work are the same as in Stoll et al. (2017b), where we analysed the dynamics of a planet embedded in a VSI turbulent disc without particles, so we briefly describe the set-up, focussing only on the dust part. For a more detailed description of the initial conditions of the gaseous disc, see Stoll et al. (2017b).

4.1 Gas component

The initial disc profile is axisymmetric and extends from 2.08 to 13 au (0.4–2.5 in code units, where the unit of length is 5.2 au). The gas moves with azimuthal velocity given by the Keplerian speed around a 1 M_⊙ star, corrected for the pressure term and rotational velocity of the coordinate system that rotates here with the orbital speed of the planet. The total disc mass is 0.01 M_⊙ and the density distribution, created by force equilibrium, is given in cylindrical coordinates (R, Z, ϕ) by $${\rho }_{\text{g}}(R,Z)={\rho }_{\text{g},0}{\left(\frac{R}{R_{\text{p}}}\right)}^p\exp \left[\frac{G{M}_{\text{s}}}{c_{\text{s}}^2}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$$(23)

where ρ_g,0 is the gas midplane density at R = 1 and p = −1.5 is the density exponent. In our case ρ_g,0 = 2.07 × 10⁻¹¹g cm⁻³ such that the vertically integrated surface density at R = 1 is Σ_g = 200 g cm⁻². The disc is modelled with a locally isothermal equation of state, and we assume a constant aspect ratio H∕R = 0.05, which corresponds to a radial temperature profile with an exponent q = −1 and T(R_p) = 121 K. For the inner and outer radial boundary, we apply reflective conditions, while outflow conditions are implemented for the vertical boundaries and periodic conditions in the azimuthal direction. We perform two sets of simulations. One set has an inviscid disc in which thesource of turbulence is given entirely by the VSI and the other uses a viscous disc in which the viscosity is given by ν = 2∕3αc_sH, where we use a constant α viscosity as derived from the VSI simulation, which is α = 5 × 10⁻⁴ (Stoll et al. 2017b).

4.2 Dust component

The solid fraction of the disc is modelled with 10⁶ Lagrangian particles divided into ten size bins as reported in Table 1. This approach has the great advantage of modelling a broad range of Stokes numbers (see Eq. (12)) self-consistently using the same model particles. The trade-off is that in the regions of low density, the resolution of the dust population is lower. However, for our study this is not a problem since we are mainly interested in the dynamical evolution of dust particles; thus we do not take into account collisions between particles or the backreaction of the dust onto the gas. We study particles with sizes from 0.1 mm up to 1 km and internal density ρ_d = 1 g cm⁻³. The particlesizes and corresponding Stokes numbers are quoted in Table 1, where the Stokes numbers are evaluated atthe planet location. The particle sizes are chosen to cover a wide range of different dynamical behaviour. The initial surface density profile of the dust particles is $${\Sigma }_{\text{d}}(r)\propto R^{-1}.$$(24)

This choice was made to have a larger reservoir of particles in the outer disc. This particle distribution leads to equal number particles in each radial ring as the grid is spaced logarithmically in the radial direction. The dust particles are placed initially at the disc midplane in the disc model with active VSI driven turbulence because the particle stirring is obtained via the turbulent mechanism. For the laminar disc, we start with a vertical distribution given by the local disc scale height and the dust diffusion coefficient. By comparing Eqs. (14) and (15), we find that, for our initial profiles, particles larger than 1 mm are going to settle to the disc midplane. The particles are introduced at the beginning of the simulation, and they are evolved with two different integrators depending on their Stokes numbers. Following the approach by Zhu et al. (2014), we adopt a semi-implicit leapfrog-like (drift-kick-drift) integrator in spherical coordinates for larger particles and a fully implicit integrator for particles well coupled to the gas. We include in Appendix B the detailed implementation of the two integrators. We do not consider the effect of the disc self-gravity on the particle evolution. Particles that leave the computational domain at the inner boundary are re-entered at the outer boundary. Accreted particles are flagged but are otherwise kept in the simulations.

4.3 Planets

We embed a planet, with a mass in the range [5, 10, 30, 100] M_⊕, orbiting a solar mass star on a circular orbit with semi-major axes a_p = 1 in code units (5.2 au). The planet does not migrate and its mass is kept fixed. To prevent a singularity close to the planet location, its gravitational potential is smoothed with a cubic expansion inside a sphere centred on the planet location with a radius given by the smoothing length d_rsm = 0.5R_H (Klahr & Kley 2006; Stoll et al. 2017b), where R_H denotes the radius of the Hill sphere: $$R_{\text{H}}=R_{\text{p}}{\left(\frac{1}{3}q\right)}^{1/3}.$$(25)

After the dust component has been evolved for 20 orbits in the computational domain, the planetary mass is slowly increased over additional 20 orbits to allow for a smooth initial phase. Each simulation was run over 200 orbital periods of the planet when the disc structure had reached a quasi-stationary state.

5 Global dust motion

In this section, we analyse the overall behaviour of the dust particles in the presence of the planet in combination with the disc turbulence. Of particular interest are the changes in the spatial distribution of the particles as induced by the planet. However, before we focus on the action of the planet we comment briefly on the dust dynamics in the disc.

5.1 Vertical dust distribution

The action of the turbulence in the disc works against the tendency of the dust particles to settle towards the midplane of the disc and leads to a vertical spreading of the particles. Concerning this particle stirring in turbulent discs, we present in Fig. 1 the vertical scale height of the particles (in units of the gas scale height H) as a function of their Stokes number, τ_s for various models. The two cases studied in this work are shown by the blue and green crosses for the VSI turbulent case and the laminar disc with stochastic particle kicks, respectively. We analysed the particle distribution for the 5 M_⊕ case at a radius of 1.8r_p averaged between 150 and 200 orbital periods. From our previous work in Stoll & Kley (2016), we know that the timescale for spreading the particles vertically in the presence of fully developed VSI turbulence is about 100 orbital periods, so near the end of our simulations (150–200 orbits) the particle distribution has reached a quasi-stationary state. Furthermore, we checked the particle vertical distribution from 100 to 150 orbital periods at the same distance and found the same profile, confirming that an equilibrium was reached. Additionally, we show the results of Stoll & Kley (2016) for locally isothermal discs using the data taken from their Table 1 (labelled SK16) and Fromang & Nelson (2009) who studied particle settling in global ideal MHD disc displaying MRI turbulence (labelled FN09).

Overplotted to the data is an approximation by Youdin & Lithwick (2007, their Eq. (28)), which can be written as (neglecting a correction factor of order unity) $$\frac{H_{\text{p}}}{H}=\sqrt{\frac{{\alpha }_z}{{\alpha }_z+{\tau }_{\text{s}}}},$$(26)

where α_z measures the vertical diffusion of the gas; see Youdin & Lithwick (2007). In Fig. 1 we use α_z = 1.737 × 10⁻³ for the fit. Equation (26) accounts for the fact that for small τ_s the particles are well coupled to the gas and the two scale heights agree, H_p = H, while larger particle settle more to the midplane of the disc and have a smaller thickness. For large τ_s the slope becomes $\propto {\tau }_{\text{s}}^{-0.5}$ as can be inferred from Eq. (26).

For the small particle sizes, our distribution is similar to that of Fromang & Nelson (2009). For their investigated particle sizes with τ_s = (10⁻⁴, 10⁻³, 10⁻²), they find a scaling $H_{\text{p}}/H\propto {\tau }_{\text{s}}^{-0.2}$ in rough agreement with our findings. For the larger particles, the slope becomes steeper than the expected $\propto {\tau }_{\text{s}}^{-0.5}$ scaling because we have reached the resolution limit in our simulations such that the particle scale height cannot be resolved anymore.

In our previous simulations of particles embedded in VSI turbulent discs (Stoll & Kley 2016), we find for the mean vertical velocity at 5 au $<v_z^2>=5\times {10}^{-6}{v}_{\text{K},\text{1au}}^2$ (normalising to the Keplerian velocity at 1 au). Using this value and H∕r = 0.05, we find for the mean vertical Mach number M_z ≈ 0.1. In Stoll & Kley (2016), we quote for the (dimensionless) eddy turnover timescale τ_e ≈ 0.2. From these we can calculate a vertical diffusion coefficient of (Youdin & Lithwick 2007) $${\alpha }_z={\tau }_{\text{e}}{M}_z^2.$$(27)

Hence, from Stoll & Kley (2016), we find τ_e ≈ 0.2, which is consistent with the value obtained by the fit for Fig. 1.

Fig. 1 Measured particle scale height, Hp, in units of the gas scale height, H, as a function of the particle Stokes number for the simulations of the turbulent VSI, and the laminar disc plus stochasistic kicks at R = 1.8rp, averaged between 150 and 200 orbital periods. Shown are the results of the runs used in this paper (labelled with blue and green crosses), and of Stoll & Kley (2016) and Fromang & Nelson (2009) indicated with light blue and black circles, respectively. Additionally, we overplot the fit of the VSI particle scale heights by Youdin & Lithwick (2007; see Eq. (26)).

5.2 Dust filtration

In Fig. 2 we plot the radial distribution of three representative size particles as a function of time. Shown are the results for both the viscous α- and turbulent VSI-disc models for all planet masses from 5 M_⊕ (top) to 100 M_⊕ (bottom). The particle sizes increase from left to right from 10.0 cm to 10 m, which corresponds to the Stokes numbers 8 × 10⁻², 1.23, and 67, respectively. Clearly visible are the different radial drift velocities of the particles that are a function of the Stokes number, τ_s. Indeed, the speeds found in our simulations are in good agreement with the theoretical expectation of Nakagawa et al. (1986), which is given by $$v_{\text{drift}}=\frac{\partial \ln p}{\partial \ln R}{\left(\frac{H}{R}\right)}^2\frac{u_{\text{K}}}{{\tau }_{\text{s}}+{\tau }_{\text{s}}^{-1}}\equiv -2\eta \frac{v_{\text{K}}}{{\tau }_{\text{s}}+{\tau }_{\text{s}}^{-1}}.$$(28)

Equation (28) indicates that the maximum speed, reached for τ_s = 1, is given by v_drift = −ηv_K, where η is typically of order ${(H/R)}^2$ and v_K is the Keplerian azimuthal velocity. The results on the drift speed for the VSI turbulent and viscous disc are very similar for all cases studied.

For the low mass 5 M_⊕ planet and small particles, one notices small disturbances near the planet (first two rows on the left), but the planets are not able to stop the particles from crossing their location. The same behaviour is also found in the 10 M_⊕ case displayed in the second two rows of Fig. 2. Focussing on the middle column, we can see the evolution of 1 m particles, which have a Stokes number of order unity (τ_s = 1.23). Their drift speed is so high that they can cross the whole computational domain in ~150 orbits, in agreement with Eq. (28); see also Stoll et al. (2017b). A change in the drift speed is also visible as the particles cross the planet co-orbital region because the Stokes number suddenly increases because of the drop in the gas density. The only exception is given by the planetesimal-sized objects (10 m, τ_s = 67, right column), which do not feel a strong gas drag, such that the planetary core can effectively perturb their orbits depleting its co-orbital region. This regime is described in Dipierro & Laibe (2017), who found that for a Stokes number number greater than a critical value $${\tau }_{\text{s},\text{crit}}\simeq 2.76\left(\frac{-\zeta }{1+ϵ}\right)\simeq 6.83,$$(29)

where ϵ = 0.01 is the dust-to-gas ratio, and ζ = ∂ lnp∂lnR = −2.5, the minimum mass to open a gap in the solid disc is $$M_{\text{crit}}\simeq 1.38{\left(\frac{-\zeta }{1+ϵ}\right)}^{3/2}{\tau }_{\text{s}}^{-3/2}{\left(\frac{H}{R}\right)}^3{M}_{\odot }.$$(30)

For our parameter space, we find that this transition happens between the particles with s = 300 cm (τ_s = 6.91) and a critical mass of 12 M_⊕ and the particles with s = 1000 cm (τ_s = 67.2) and a critical mass of 0.4 M_⊕. We see that only the second sample of particles is depleted from the co-orbital region of the small mass planets, confirming their analytical prescription.

The third two rows of Fig. 2 show the particle evolution for a 30 M_⊕ planet. This planetary mass can change the final particle distribution dramatically. A gap is already visible for the 10 cm particles, while for the particles with Stokes number of order unity (central column) the planet acts as a barrier and can filtrate the dust in the outer disc. This effect is due to the formation of a pressure maximum beyond the planet where small particles are trapped (Paardekooper & Mellema 2006). After ~ 100 orbits the particles are located either in the pressure bump close to the planet position or in the outer disc. This dust filtration leads to a strong reduction of particles inside the planetary orbit, and hence the number of particles leaving through the inner radius is strongly diminished. As those particles are re-entered at the outer boundary, eventually this results in a shut-off of the flow of particles in the outer disc. We decided not to have a constant inflow of particles because they would only end at the pressure bump as all the others. The particle concentration could become a sweet spot to have a second generation of planets due to streaming instability, but it is beyond the purpose of this paper to study this high dust density regions in more detail.

For the planetesimal-sized objects, the 30 M_⊕ planet can open a deeper and wider gap compared to the small mass cases. Only planets greater than 10 M_⊕ are able to filter the pebble-sized particles efficiently. This result confirms the value obtained by Lambrechts et al. (2014) in which they defined the pebble isolation mass around 20 M_⊕ for similar initial disc condition. As we have seen, a planet that modifies the pressure profile in the disc can effectively stop the inward drift of certain size of dust particles. In a related scenario, such a dust filtration is believed to explain the observation of a class of transition discs (Type 2), where the dust is highly depleted in the inner region of the protoplanetary disc while the gas accretion rate onto the star remains high (for a recent review on the subject, see Ercolano & Pascucci 2017). The last two rows show the particle evolution for the 100 M_⊕ planet where the dynamical behaviour of the dust particles is very similar to the 30 M_⊕ case but thegap opens earlier, and so the simulations reach a stationary state on a shorter timescale.

Fig. 2 Spatial distribution of the dust particles as a function of time for the different planetary masses and for three representative particle sizes. The Stokes numbers from left to right are 0.08, 1.23, and 67.2.

5.3 Gap opening

A planet can open a gap in the dust disc even if no clear gap appears in the gas distribution (see e.g. Picogna & Kley 2015). We analysed the radial distribution of the dust population by splitting the computational domain into 400 logarithmically spaced bins and following its evolution with time. In Fig. 3 we plot the distribution at the end of the simulation (200 planetary orbits) for the same three representative size particles as in the previous plot. In the left column, the 10 cm dust particles (corresponding to a Stokes number of τ_s = 0.08) do not show a strong perturbation by the presence of the small mass planets (in the first two rows, corresponding to 5 and 10 M_⊕). A clear gap appears starting with the 30 M_⊕ planet (third row), where the influence of the VSI (represented by a blue line) favours the formation of deeper gaps.

The intermediate case of meter-sized particles, which corresponds to a Stokes number of τ_s = 1.2, represents the fastest evolving particles in the simulation; see Eq. (28). As shown in the central column, the distribution is strongly affected by the vertical motion of the VSI where a bunching behaviour can be noticed (Stoll & Kley 2016). This feature cannot be reproduced by the viscous α-disc model. In the 30 M_⊕ planet case (third column) we see that the VSI also leads to a faster dust filtration process, which is already completed for the 100 M_⊕ case in which the inner disc is practically devoid of particles. For the massive planets, one notices that a large number of small and large particles remain at the co-rotation location. These are particles that collect near the two Lagrange points L4 and L5.

For a planetesimal-sized object of 10 m in the third row (corresponding to a Stokes number of τ_s = 67) the gas influence is negligible. A gap is visible in the distribution already for the small mass planets due to their gravitational interaction with the particles. In this case, the VSI does not affect the evolution of planetesimals and yields results identical to the viscous case.

The effect of filtering and gap formation can be understood in terms of the angular velocity distribution in the disc around the planet. The onset of gap formation leads to a super-Keplerian flow just outside of the planet, which coincides with a maximum in the radial pressure distribution. The angular velocity is shown in Fig. 4, where the ratio of v_ϕ∕v_kep is shown forthe four different planetary masses. Outside of the planet the ratio is slightly smaller than one owing to the pressure support of the gas. For the planet masses displayed the super-Keplerian motion begins to show for the 30 M_⊕ case. Hence, as expected the filtering process is directly related directly to the maximum in the angular velocity. The property that particles (with unit Stokes number) cannot be accreted above a critical planet mass is referred to as the pebble isolation mass. In Appendix A we present additional simulations to confirm that our simulations have been run long enough to draw this conclusion about the super-Keplerian motion.

5.4 Planet–solid disc interaction

The planet is not only able to open a gap in the dust and gas disc, but it can also generate non-axisymmetric features in their distribution that might be observable with modern observational facilities. In Fig. 5 we show the surface density distribution of the dust population after 80 planetary orbits. The spiral arms that are typically generated by embedded planets are only (barely) visible for the most massive planet (bottom row) and the smallest particles (left column) that are well coupled to the gas dynamics. For the 100 M_⊕ planet, strong vortices are created in the gas disc for the VSI case and less so for the laminar case (Stoll et al. 2017b). Because they are pressure maxima, particles tend to accumulate in these vortices, which is reflected in the corresponding particle distributions as seen in the bottom left part of the plot.

Also visible in the middle column is the strong effect that the VSI has on pebble-sized particles creating regions where particles are collected, as shown in Stoll & Kley (2016). For the large planets we see, as shown in more detail before, in the bottom part of the middle column that a sort of transition disc is formed since the planet has reached the pebble isolation mass, stopping the influx of meter-sized particles (as seen also in Ayliffe et al. 2012). On the other hand, for the planetesimal-sized objects, shown in the right panel, we observe ripples close to the planet location due to the excitation of the eccentricity in the dust particles by the planet that the gas is not able to effectively damp on a short timescale. Furthermore, the planetesimals are collecting in the Lagrangian points (L4 and L5) in front and behind the planet location, and their density is enhanced at the outer 2:1 mean motion resonance with the planet, visible in the upper part of the last column for the 100 M_⊕ planet.

The eccentricity distribution indicated in Fig. 6 shows only small differences between the α- and VSI-disc models, while they are much more pronounced for the inclination distribution shown in Fig. 7. The excitations at the resonance locations are visible for the biggest size objects in the lower panel. From Eq. (22), we find that the region where the chaotic behaviour prevents the resonances from stopping planetesimal objects starts at a radial distance of 0.1478 from the planet location for the 100 M_⊕ planet and 0.1048 for 30M_⊕. These values roughly correspond to the 5:4 and 7:6 MMRs with the planet. From Fig. 6, where the location of the major MMRs are plotted, we can confirm this finding. Planetesimal objects cannot reach the region inside the 5:4 resonance with the 100 M_⊕ planet, while for the30 M_⊕ planet the bodies are able to reach the 7:6 resonance where the chance of being accreted by the planetary core is much higher.

In Table 2 we also report the mean values of all the solids in the disc at the end of the simulation. The eccentricity and inclinations are on average higher for the VSI discs with respect to the laminar discs in the small mass planets. This effect can be explained by the highly anisotropic turbulence nature of the VSI, which is able to stir up dust particles exciting their orbitalelements more efficiently. This trend, however, is inverted for the high mass planets, possibly because the VSI strength is partly reduced by the presence of massive planets, and the production of vortices in which the particles tend to be collected.

The viscous α-disc model cannot correctly reproduce the distributions observed in the VSI disc because we assumed a constant α-value throughout the whole computational domain, and did not distinguish between radial and vertical angular momentum transport. However, as shown recently, the VSI turbulence behaves strongly anisotropic with a large difference between radial and vertical transport (Stoll et al. 2017a). Since the turbulent kicks in the particle motion are generated based on the constant alpha value, they over- or under-predict the turbulence efficiency in the laminar disc resulting in a different particle scale height, and thus inclination. Nevertheless, this does not seem to play a crucial role in the solid accretion rate to the planet from small turbulent velocities.

Fig. 3 Histogram of the dust surface number density distribution as a function of radius after 200 planetary orbits for three representative dust sizes in the turbulent (blue) and laminar (green) case. The Stokes numbers from left to right are 0.08, 1.23, and 67.2.

Fig. 4 Azimuthal gas speed in units of the Keplerian speed as a function of radius for the different planetary masses and models (VSI = red, alpha-disc = green). When the gas speed becomes super-Keplerian outside the planet location, the dust-filtration process occurs and the pebble isolation mass is reached.

Fig. 5 Surface density distribution of the dust particles after 80 planetary orbits for the different planetary masses and for three representative particle sizes. The Stokes numbers from left to right are 7.79 × 10−5, 1.23, and 67.2.

Fig. 6 Eccentricity distribution for 3 representative size particles in the VSI (blue) and alpha-disc (green) model at the end of the simulation. The first six inner and outer first-order MMRs are overlaid.

6 Solid accretion rate

In order to detect the particles that accrete onto the planet, we adopted two different approaches, depending on the ratio of their Stokes number and the time they spend inside the Hill sphere which, for pebble-sized particles with τ_s is t_enc = R_H∕Δv, where Δv is the relative velocity between the particle and the planet. For particles with stopping time shorter than t_enc, their trajectory close to the planet location is determined primarily by the drag force. Whether the particle is then accreted depends on the relative strength of the gravitational attraction and the drag force. If the drag force dominates, it can sweepout a particle even if it is gravitationally bound. Thus, we checked if the timescale for gravitational attraction t_g = Δv∕g is shorter than four times the timescale for the stopping time. The factor of four stems from the results of Ormel & Klahr (2010), who had found in their numerical simulations that only a small deflection of a fourth of the velocity is needed to accrete the particle.

Larger particles lose only a small amount of momentum through drag when they cross the Hill sphere. Thus, we checked for whether particles inside the Hill sphere are bound by the gravity of the planet, which is the case if the particle has not enough kinetic energy to leave the Hill sphere, that is $$e_{\text{kin}}+e_{\text{grav}}<e_{\text{grav}}(R_{\text{H}}).$$(31)

Both approaches also agree in the transition region where the stopping time is similar to t_enc. We checked for these conditions every tenth of a planetary orbit and we flagged as accreted the particles that fulfil the previous criteria without removing these particles from the computational domain.

In Fig. 8 we show the number of accreted particles per orbit for the different planet masses as a function of time. In the initial phase of the simulations, the number of accreted particles increases while the mass of the inserted planets grows to their final value (within 20 orbits). After that, the number of accreted particles drops continuously and settles roughly to a constant value for the lower mass planets, as there is at least for the faster drifting particles a continuous supply from the outer disc (see Fig. 2 middle column). For the larger mass planets, the accretion rates are further reduced as they have reached their isolation mass for the particles with Stokes number around unity. Additionally, for the large planets, the small and large particles drift very slowly and, after the particles within the horseshoe region have been accreted, the new inflow from the outer disc is very slow.

In Fig. 9 we show the number of accreted particles of various sizes for different planet masses summed over 50 planetary orbits, from t = 100 to 150. Several interesting features can be noted. The number of accreted particles peaks for particles in the range between 30 and 300 cm (corresponding to pebble-sized objects with Stokes number of order unity) for the small planetary masses, while these particles areeffectively accreted and filtered for the two higher mass planets. The total number of accreted particles adds up to about 5000 for each of the lower mass planets in agreement with the results shown in Fig. 8. Moreover, although the effect of VSI seems marginal it shows in nearly all cases a slightly higher solid accretion rate than for the viscous α-disc model. As the difference is rather small we may conclude that our modelling of the stochastic motion of particles in discs also gives reasonable results for the accretion rates of the particles onto embedded planets.

To obtain an actual mass accretion rate onto the planet, we need then to convolve this result with a dust size distribution.Birnstiel et al. (2012) showed that the size distribution is very steep for Stokes numbers less than 0.1. At that point, there is a gap due to the so-called meter-sized barrier. This effect removes the peak point of our accreted particles for the small mass planets and renders the accretion growth for smaller dust particles even steeper. For the bigger size objects, there are far fewer constraints from models. The leading theory of streaming instability predicts that the peak of formed planetesimals is 10–100 km-sized objects with a tail in the distribution also to kilometre-sized objects, which represent our bigger size objects in the simulation (Simon et al. 2016, 2017).

Efficiency of pebble accretion. Very important for the mass growth of a planetary core is the efficiency, P_eff, of the accretion process, i.e. the number of accreted particles onto the planet divided by the number of particles that would otherwise drift across the location of the planet in an unperturbed disc. Following Ormel & Klahr (2010), we define this efficiency as $$P_{\text{eff}}=\frac{{\dot{M}}_{\text{acc}}}{{\dot{M}}_{\text{drift}}},$$(32)

where Ṁ_acc is the actual accretion rate of solids onto the protoplanet and Ṁ_drift is the particle drift rate through the disc, given by $${\dot{M}}_{\text{drift}}=2\pi r{\Sigma }_{\text{p}}{v}_{\text{drift}},$$(33)

where Σ_p is the particle surface density and v_drift is given by Eq. (28). The quantity P_eff in Eq. (32) is, in fact, the probability that a particle that drifts through the disc is accreted by the protoplanet.

When analysing the data in this way, we encountered the problem that for the particles with very small and very large Stokes numbers, the drift velocities v_drift are very small (in agreement with Eq. (28)) such that the calculated efficiencies became very high because of the relatively large amount of particles accreted. The reason for this lies in the fact that the accreted small and large particles originate primarily from the horseshoe region and did not migrate to the planet, which is only a transient effect visible in the initial phaseof the simulations. Hence, we decided to use an alternative way of measuring the efficiency of particle accretion from our simulation, which is illustrated in Fig. 10. To measure the accretion efficiency of particles on to a growing planet we monitor the evolution of particles that are initially in a ring just outside the planet. We use the radial range from x_s to 2x_s, where x_s is the horseshoe half-width as defined, by Paardekooper & Papaloizou (2009), as $$x_{\text{s}}=1.68R_{\text{p}}{\left(\frac{q}{h}\right)}^{1/2}.$$(34)

Radial drift brings the particles into the co-orbital region of the planet. Some of the particles are accreted (and marked so), while othersare able to cross the horseshoe region and are not accreted. These latter particles enter the inner region of the domain and are called the survivors. The results of using this procedure for our simulations are shown in Fig. 11 for the different planet masses and particle sizes.

For the two larger planet masses (30 M_⊕ and 100 M_⊕), the results are not very meaningful because the total number of accreted particles is very small as they have already reached their isolation masses. Hence, for the growth of planets, we focus on the two smaller mass planets (5 M_⊕ and 10 M_⊕). In both cases the lowest accretion efficiency is reached for particles with Stokes number τ_s ≈ 1. For smaller and larger particles, the efficiency rises but due to the very slow drift speeds becomes unreliable for very small (τ_s < 10⁻²) and large particles (τ_s > 10−100). The result shows that particles with fast radial drift (τ_s around unity) have small accretion efficiencies because a high percentage of particles can cross the planetary orbit. For Stokes number τ_s = 1, we find P_eff ≈ 1.6% for the 5 M_⊕ and around 3% for the 10 M_⊕ planet.

We can compare our measurements to the results of previous estimates using particle trajectories in the vicinity of the planet. We use our set-up with the protoplanet located at R_p = 5.2 au. The radial drift speed is then v_drift = 30 m s⁻¹ and for the particle density, we have with 1% in solids Σ_p = 2 g cm⁻². To calculate the efficiency we use Eq. (37) from Ormel & Klahr (2010) with log₁₀P_col = 0.5 for the (dimensionless) collision rate, and drop the 3D correction. The specific collision rate P_col is given in this case by (see their Eq. (3)) $${\dot{M}}_{\text{col}}=P_{\text{col}}{\Sigma }_{\text{p}}{R}_{\text{H}}{v}_{\text{H}},$$(35)

where $$v_{\text{H}}={\Omega }_{\text{K}}{R}_{\text{H}}$$(36)

is the Hill velocity. For an alternative comparison, we use the accretion rate in the Hill regime as given in Lambrechts & Johansen (2012), also in the 2D version: $${\dot{M}}_{\text{Hill}}=2{\Sigma }_{\text{p}}{R}_{\text{H}}{v}_{\text{H}},$$(37)

i.e. the ratio of these two is given by 10^0.5∕2 ≈ 1.6. To compare directly to other estimates, we can transform our obtained accretion rates in terms of the Hill accretion rate and define $$f_{\text{Hill}}\equiv \frac{{\dot{M}}_{\text{acc}}}{{\dot{M}}_{\text{Hill}}}=P_{\text{eff}}\frac{{\dot{M}}_{\text{drift}}}{{\dot{M}}_{\text{Hill}}}=P_{\text{eff}}\pi \frac{R_{\text{p}}{v}_{\text{drift}}}{R_{\text{H}}{v}_{\text{H}}}.$$(38)

Applying the definitions of v_drift, R_H and v_H from Eqs. (28), (25) and (36), one obtains with τ_s = 1 and $\eta ={(H/R)}^2$ $$f_{\text{Hill}}\approx 6.5P_{\text{eff}}{\left(\frac{h}{q^{1/3}}\right)}^2.$$(39)

Using this equation and our findings for the particle accretion, we obtain the results quoted in Table 3. We notice that our estimates are about 50% lower than the 2D approximation for the Hill case, but here we also conside the third dimension, which can lower the amount of particles within the reach of an embedded planet. In Fig. 11 we compare also the accretion efficiency with the prescription from Lambrechts & Johansen (2014), that is $$P_{\text{eff},\text{LJ}}\simeq 0.034{\left(\frac{{\tau }_{\text{s}}}{0.1}\right)}^{-1/3}{\left(\frac{M_{\text{c}}}{M_{\oplus }}\right)}^{2/3}{\left(\frac{r}{10\text{au}}\right)}^{-1/2},$$(40)

where we see that although our result are slightly lower, the scaling with the Stokes number and planetary masses for the intermediate τ_s, for which our approach was reliable, is consistent.

From our simulation, we may estimate the mass doubling time for our low mass planets, for which we use $$t_{\text{double}}=\frac{M_{\text{core}}}{{\dot{M}}_{\text{acc}}}.$$(41)

With our results on P_eff we find t_double = 20 000 yr for the small mass planets. The region in the disc that can supply this amount of solid material within the time t_double extends to roughly 36 au assuming a constant surface density of solids, Σ_p = 2 g cm⁻².

Fig. 7 Inclination distribution for 3 representative size particles in the VSI (blue) and alpha-disc (green) model at the end of the simulation. The first six inner and outer first-order MMRs are overlaid.

Fig. 8 Number of accreted particles per orbit over time for the four different planet masses. Shown is the total number summed over all size bins.

Fig. 9 Accreted particles as a function of the particle stopping time, where their approximate size is also reported in the top x-axis, integrated over 50 planetary orbits (from t = 100 to 150). The VSI (solid line) and alpha disc model (dashed line) are compared.

Fig. 10 To measure the accretion (and survival) efficiency, the evolution of particles initially within a ring outside the planet co-orbital region (defined as in Eq. (34)) is monitored.

Fig. 11 Efficiency of accreted (red line) and survived particles (green line) as a function of Stokes number for different planetary masses. The laminar disc run is represented with solid lines, while the VSI run with dashed lines. The fit from Lambrechts & Johansen (2014) is overplotted for intermediate τs values with a dashed grey line (see Eq. (40)).

7 Discussion

After having presented our main findings we now compare our results to other studies of particles embedded and accreted onto a planetary core in MHD turbulent discs. Then we shall discuss possible limitations of our simulations because in performing the simulations we had to use several approximations to make them feasible.

7.1 Comparison to MHD simulations

Recently, Xu et al. (2017) studied the accretion of particles onto small planets embedded in discs exhibiting magnetically driven turbulence. These authors considered a local shearing box centred at the growing core and studied three different types of discs: a fully turbulent disc using ideal MHD, a less turbulent disc with ambipolar diffusion, and a non-turbulent hydrodynamic disc. For all three cases, particles with various Stokes numbers were injected to the flow after reaching equilibrium, and the accretion rate of particles onto the core was measured. The measured particle accretion rates were then compared to that in the 2D Hill regime as estimated by Lambrechts & Johansen (2012). In all cases, these authors found for Stokes numbers around unity that their measured rates agree very well (within ≈ 10%) with the 2D accretion rate of Eq. (37) in Lambrechts & Johansen (2012) and argue for a high accretion efficiency. In our simulations we measure the absolute accretion efficiency as defined in Eq. (32) and this is much lower than one for 10⁻² ≤ τ_s ≤ 1, with an additional drop towards τ_s = 1; see Fig. 11. To compare directly to Xu et al. (2017), we can use Eq. (38) and the values quoted in Table 3, which shows that our calculated accretion rates are about 50% smaller than theirs.

The higher rates obtained by Xu et al. (2017) may be a result of the relatively small vertical and radial extent of the disc in their simulations (only one H in each direction) and the fact that they do not consider a stream of particles through their domain. Hence, they measure the transient accretion rate of particles present initially in the computational box. In our case we measure the accretion rate of particles in an equilibrium situation. The absolute accretion efficiency, as defined in Eq. (32), is more important; we measure this rate directly. This is actually very low, but the combination with the large drift speed allows us to tap a larger reservoir of particles making the accretion of τ_s = 1 particles very useful in the overvall growth process.

7.2 Planet migration

In our simulations, the location of the planet is kept fixed at 5.2 au. However, planets interact gravitationally with their protoplanetary disc, and this usually results in an inwards migration through the disc, which may modify our conclusions about the capture efficiency of particles. To check the impact of planet migration, we compare the particle drift rates to the expected planet migration for which we use the 3D results of D’Angelo & Lubow (2010) for low mass planets in the linear regime, as stated in Kley & Nelson (2012): $$t_{\text{mig}}=C\frac{M_{\star }^2}{M_{\text{p}}{\Sigma }_g(r_{\text{p}})r_{\text{p}}^2}{\left(\frac{H}{r}\right)}^2{\Omega }_{\text{K}}^{-1},$$(42)

where C = 1∕(1.36 + 0.62 β_Σ + 0.43 β_T). The value β_Σ = p − 1 is the coefficient of the surface density profile, while β_T = q is the coefficient of the temperature radial profile. In our models we used β_Σ = 0.5 and β_T = 1.0 (see Sect. 4.1) and then we find C = 0.48. In Fig. 12 we compare the dust drift velocity measured at the beginning of the simulations, when the planet has not yet perturbed the disc structure, to the theoretical drift speed (as obtained from Eq. 28; black solid line), and the migration speed of the two small planets (two dashed lines). We can see that only the tails of the particle size distribution drift slower than the planets. However, pebble-sized particles have a migration speed orders of magnitude faster than the planet. Hence, we conclude that any planet migration does not influence our results for pebble accretion. For larger mass planets that migrate with Type II migration, the drift speed slows down considerably, and their masses are well above the pebble isolation mass. The impact of non-circular planetary orbits on the pebble accretion efficiency was studied recently by Liu & Ormel (2018) who found that it can be increased slightly for moderately eccentric orbits.

Fig. 12 Comparison between the average drift speed of the different dust size particles. These values are calculated at the end of the simulation at R = 1.8 with the analytical drift speed (black solid line), the migration speed of a 5 M⊕ (dotted line), and a 10 M⊕ (dashed line) mass planet. The error bars reflect the impact of radial diffusion for the various stopping times. The smaller particle sizes in the VSI simulation are migrating outwards rather then inwards due to the inverse meridional circulation (Stoll et al. 2017a).

7.3 Equation of state

In our simulations, we used an isothermal equation of state and now briefly discuss a possible impact of including radiative effects. The inclusion of radiative transfer leads to finite cooling times of the gas that lowers the efficiency of the VSI-driven turbulence (Nelson et al. 2013). In full simulations that include radiative transfer, it has been shown that in irradiated discs an efficiency of α ≈ 10⁻⁴ can be reached (Stoll & Kley 2014, 2016), while Flock et al. (2017) find a somewhat smaller value. All those simulations apply to larger distances from the star and it remains to be seen what the VSI-efficiency is at shorter distances from the star. In any case, a reduced turbulence level leads to a concentration of the dust particles in the midplane, which might enhance the accretion process. On the other hand, the inclusion of radiative transfer allows for additional disc heating by the planet (by the spiral waves), which enhances the disc temperature and might lead to partial evaporation of the particles. However, for the lower mass planets, for which the dust accretion efficiency is higher, the effect on the disc is not be that strong and we do not expect a large impact. Additionally, the dust clearing around the planet alters the opacity of the medium and hence the radiative transport. These impacts of radiative transport and the link to observations have to be investigated in more detail by future simulations.

7.4 Dust feedback

In our simulations we have neglected the backreaction of the dust onto the gas. Within a disc without an embedded planet the particle concentrations are such that the dust density remains typically smaller than the gas density, given an initial dust to gas ration of 1/100. In the presence of a planet, this is not true in the case of filtering because then the dust density can equal the gas density near the pressure maximum. This situation has recently been explored by Weber et al. (2018), who showed that dust feedback can potentially displace the gas density maximum, and thus the pressure maximum, outwards. Additional dust diffusion (for example from disc turbulence) can smooth the density peak of the dust distribution, altering the dust filtration process for particles with Stokes numbers around unity. However, concerning the filtration ability, which also affects pebble accretion onto the planet, they did not observe any difference by adding dust feedback. Hence, we conclude that dust feedback does not impact our results significantly. In a realistic scenario, where a dust size distribution is present, this effect is even less pronounced. The impact of dust feedback onto the dust dynamics in the dust trap itself will have to be investigated in more detail in the future.

7.5 Numerical convergence

Our hydrodynamical simulations are performed with one numerical resolution as given in Table 1. This is based on our results in Stoll et al. (2017b) in which we studied the effect of doubling the grid resolution on the VSI and found no noticeable differences in the disc dynamics. The calculated α_SS close to the inner boundary (see their Fig. 1) was increased marginally, however it had no effect at the location of the planet. We also checked in Stoll et al. (2017b) that the torque acting onto the planet reached a constant value (see their Fig. 6), guaranteeing that the model was run long enough for the disc-planet system to have reached a quasi-stable state. The obtained torque distributions on the planet were also identical for the standard resolution (used here) and the simulations with doubled resolution. The dust dynamics is not affected by numerical resolution since we did not take into account their backreaction on the gas. Thus, we do not expect that our results on accretion efficiencies and dust dynamics in the vicinity of the planet are impacted by resolution that is too low.

8 Conclusions

In this study we have modelled the dynamics of a broad range of solid particles, ranging from 100 μm dust particles to kilometer-sized planetesimals, interacting with a growing planetary core in a 3D globally isothermal disc. By modelling a global disc, we were able to take into account the effect of MMRs for planetesimal-sized objects and the enhanced particle density in spiral arms and vortices. The turbulence driving the disc evolution has been modelled both self-consistently through the VSI instability and with an alpha parameter derived from the VSI simulation where the turbulence has been recreated in the particle dynamics by adding random kicks to their motion. We determined the solid accretion rate onto the planet after it reaches a stable state averaging the values over 50 planetary orbits in Fig. 9.

The actual growth rates in particles that a planet can achieve depends on the particle size distribution. In our study, we sample the particle dynamics in ten different size bins. One can convolve this result with a model of dust size distribution to obtain a mass accretion rate onto the planet. We observed a peak in the absolute number of accreted particles in the range of pebble-sizedobjects (100 cm) with Stokes number of order unity, but the strength of this effect depends strongly on the chosen dust size distribution. Concerning the accretion efficiency, we find that the minimum efficiency is reached for particles with τ_s = 1, where P_eff= 0.016 and 0.03 for the planets withmasses 5 M_⊕ and 10 M_⊕, respectively. For smaller and larger particles, the efficiency rises but due to the rapid inwards drift of particles with τ_s = 1, we find that the optimal particle size for pebble accretion for our massive cores is about 1 m at the orbit distance of about 5 au. If all the solid material in the disc was this size range, the mass doubling time would be around 20 000 yr. We find that the obtained accretion efficiencies are very similar for the VSI turbulent disc and the laminar disc models, one has to keep in mind however to add the stochastic kicks to the particles for the viscous model. This similarity can be attributed to the fact that the overall turbulence generated by the VSI is relatively weak such that the disc structures are very similar, despite the occurrence of vortices in the VSI case.

The accretion efficiency found in our simulations agrees reasonably well with previous results, for example, the 2D approximations of Ormel & Klahr (2010) and Lambrechts & Johansen (2012) or the 3D turbulent simulations of Xu et al. (2017) who found similar results for particles with τ_s = 1. To obtain the efficiencies of very small or large particles exactly, one needs longer integration time due to the very slow radial drift. Concerning the pebble isolation mass of a growing planet we confirm that the occurrence of a pressure maximum in the gas created by the planet is sufficient to filter particles with Stokes numbers of unity efficiently, at least for the relatively weak VSI turbulence. Hence, using purely hydrodynamical studies the dependence of the isolation mass on the viscosity and pressure scale height of the disc has been examined recently to obtain scaling relations (Bitsch et al. 2018).

The treatment of the turbulence adopted for the particles in the laminar disc produced accretion rates in good agreement with those of the self-consistent VSI treatment. The impact of radiative transfer within the disc and the migration of the planet through the disc were not treated. These topics need to be addressed in future work.

Acknowledgements

We thank the anonymous referees for useful comments and suggestions. The very helpful discussions with Chris Ormel are gratefully acknowledged. G. Picogna acknowledges the support through the German Research Foundation (DFG) grant KL 650/21 within the collaborative research programme “The first 10 Million Years of the Solar System”. M.H.R.S. acknowledges the support through the (DFG) grant KL 650/16. This work was performed on the computational resource ForHLR I funded by the Ministry of Science, Research and the Arts of Baden-Württemberg, and the DFG. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”.

Appendix A Long-term two-dimensional integrations

In Fig. 4 above we displayed the ratio of the angular velocity of the gas to the Keplerian velocity for the differentplanet masses and compared the turbulent case to the viscous laminar case. From this, we argued that the occurrence of super-Keplerian flow can be taken as an indication of having reached the isolation mass for that particular planet. However, this argument is only valid if the equilibrium of the flow has already been reached anddoes not change in time significantly anymore. For a viscous disc with a kinematic viscosity ν, the viscous timescale is given by $${\tau }_{\nu }=\frac{\Delta r^2}{\nu },$$(A.1)

Fig. A.1 Ratio of the azimuthal gas velocity to the Keplerian value for a 2D disc with an embedded 10 M⊕ planet at different times after insertion of the planet. The simulation in the left panel uses a viscosity of α = 10−4 and the right α = 5 × 10−4.

where Δr² is the spatial region under consideration. Assuming an α-type viscosity with ν ~ αΩ_KH² and Δr = f_HH one finds forthe viscous timescale $${\tau }_{\nu }=\frac{f_{\text{H}}^2}{2\pi \alpha }{P}_{\text{K}},$$(A.2)

where P_K is the Keplerian period. The maximum of Ω occurs roughly at a distance of Δr = 2H in our case and for f_H = 2 Eq. (A.2) gives an equilibration time of over 1200 orbits. To run our computations in full 3D for such a long timescale would have been too costly and we investigated this issue by performing comparison 2D simulations of planets embedded in flat discs. For these, we used a 10 M_⊕ planet and two different effective viscosities, α = 5 × 10⁻⁴ and a lower viscosity case using α = 10⁻⁴. The results are shown in Fig. A.1. While for the low viscosity case the flow becomes super-Keplerian, the model with α = 5 × 10⁻⁴ remains sub-Keplerian throughout. From this, we infer that in VSI turbulent discs with an effective α = 5 × 10⁻⁴ the isolation mass is indeed above 10 M_⊕ as found in the full 3D simulations presented above.

Appendix B Integrator

We used two different integrators to evolve the Lagrangian particles, based on their coupling with the gas dynamics.

B.1 Semi-implicit integrator in polar coordinates

The dynamics of particles well coupled to the gas, which have a stopping time much smaller than the time step adopted to evolve the gas dynamics, is described by adopting the semi-implicit Leapfrog (Drift-Kick-Drift) integrator described in Zhu et al. (2014) in polar coordinates. This method guarantees the conservation of the physical quantities for the long-term simulations performed in this paper and, at the same time, it is faster than an explicit method. The variables are updated beginning with a first half drift: $$\begin{array}{llllll}{L}_{\text{r},n+1}\hfill & =L_{\text{r},n},\hfill & L_{\theta ,n+1}\hfill & =L_{\theta ,n},\hfill & L_{\varphi ,n+1}\hfill & =L_{\varphi ,n}.\hfill \\ r_{n+1}\hfill & =r_n+L_{\text{r},n}\frac{\text{d}t}{2},\hfill & {\theta }_{n+1}\hfill & ={\theta }_n+\frac{1}{2}\left(\frac{L_{\theta ,n}}{r_n^2}+\frac{L_{\theta ,n+1}}{r_{n+1}^2}\right)\frac{\text{d}t}{2},\hfill & {\varphi }_{n+1}\hfill & ={\varphi }_n+\frac{1}{2}\left(\frac{L_{\varphi ,n}}{R_n^2}+\frac{L_{\varphi ,n+1}}{R_{n+1}^2}\right)\frac{\text{d}t}{2},\hfill \end{array}$$

followed by a kick step: $$\begin{array}{ll}{r}_{n+2}\hfill & =r_{n+1},{\theta }_{n+2}={\theta }_{n+1},{\varphi }_{n+2}={\varphi }_{n+1}.\hfill \\ L_{\varphi ,n+2}\hfill & =L_{\varphi ,n+1}+\frac{\text{d}t}{1+\frac{\text{d}t}{2t_{\text{s},n+1}}}\left[-{\left(\frac{\partial \Phi }{\partial \varphi }\right)}_{n+1}+\frac{L_{\varphi ,\text{g},n+1}-L_{\varphi ,n+1}}{t_{\text{s},n+1}}\right],\hfill \\ L_{\theta ,n+2}\hfill & =L_{\theta ,n+1}+\frac{\text{d}t}{1+\frac{\text{d}t}{2t_{\text{s},n+1}}}\left[\frac{1}{2}\frac{\cos ({\theta }_{n+2})}{\sin ({\theta }_{n+2})}\left({\left(\frac{L_{\varphi ,n+1}}{R_{n+1}}\right)}^2+{\left(\frac{L_{\varphi ,n+2}}{R_{n+2}}\right)}^2\right)-{\left(\frac{\partial \Phi }{\partial \theta }\right)}_{n+1}+\frac{L_{\theta ,\text{g},n+1}-L_{\theta ,n+1}}{t_{\text{s},n+1}}+\right],\hfill \\ L_{\text{r},n+2}\hfill & =L_{\text{r},n+1}+\frac{\text{d}t}{1+\frac{\text{d}t}{2t_{\text{s},n+1}}}\left[\frac{1}{2r_{n+2}}\left({\left(\frac{L_{\varphi ,n+1}}{R_{n+1}}\right)}^2+{\left(\frac{L_{\varphi ,n+2}}{R_{n+2}}\right)}^2+{\left(\frac{L_{\theta ,n+1}}{r_{n+1}}\right)}^2+{\left(\frac{L_{\theta ,n+2}}{r_{n+2}}\right)}^2\right)-{\left(\frac{\partial \Phi }{\partial r}\right)}_{n+1}+\frac{L_{\text{r},\text{g},n+1}-L_{\text{r},n+1}}{t_{\text{s},n+1}}\right],\hfill \end{array}$$

and, for the laminar disc case, also a random kick, i.e. $$\begin{array}{llllll}{r}_{n+2}\hfill & =r_{n+2}+\delta r_{\text{d},\text{T}},\hfill & {\theta }_{n+2}\hfill & ={\theta }_{n+2}+\delta {\theta }_{\text{d},\text{T}},\hfill & {\varphi }_{n+2}\hfill & ={\varphi }_{n+2}+\delta {\varphi }_{\text{d},\text{T}}.\hfill \end{array}$$

Finally, a second half drift follows as the first half drift.

B.2 Fully implicit integrator in polar coordinates

For particles with stopping time much smaller than the numerical time step, the drag term can dominate the gravitational force term, causing the numerical instability of the integrator. Thus, it is necessary to adopt a fully implicit integrator following Bai & Stone (2010) and Zhu et al. (2014).

We begin with a predictor step for the particle positions: $$\begin{array}{llllll}{r}^{\prime }\hfill & =r_n+L_{\text{r},n}\text{d}t,\hfill & {\theta }^{\prime }\hfill & ={\theta }_n+\frac{L_{\theta ,n}}{r_n^2}\text{d}t,\hfill & {\varphi }^{\prime }\hfill & ={\varphi }_n+\frac{L_{\varphi ,n}}{R_n^2}\text{d}t,\hfill \end{array}$$

followed by a shift for the momenta $$\begin{array}{ll}{L}_{\varphi ,n+1}\hfill & =L_{\varphi ,n}+\frac{\text{d}t/2}{1+\text{d}t\left(\frac{1}{2t_{\text{s},\text{n}}}+\frac{1}{2t_{\text{s},\text{n+1}}}+\frac{\text{d}t}{2t_{\text{s},\text{n}}{t}_{\text{s},\text{n+1}}}\right)}\cdot [-{\left(\frac{\partial \Phi }{\partial \varphi }\right)}_n-\frac{L_{\varphi ,n}-L_{\varphi ,\text{g},n}}{t_{\text{s},n}}\hfill \\ \hfill & +(-{\left(\frac{\partial \Phi }{\partial \varphi }\right)}_{n+1}-\frac{L_{\varphi ,n}-L_{\varphi ,\text{g},n+1}}{t_{\text{s},n+1}})\left(1+\frac{\text{d}t}{t_{\text{s},n}}\right)],\hfill \\ L_{\theta ,n+1}\hfill & =L_{\theta ,n}+\frac{\text{d}t/2}{1+\text{d}t\left(\frac{1}{2t_{\text{s},\text{n}}}+\frac{1}{2t_{\text{s},\text{n+1}}}+\frac{\text{d}t}{2t_{\text{s},\text{n}}{t}_{\text{s},\text{n+1}}}\right)}\cdot [-{\left(\frac{\partial \Phi }{\partial \theta }\right)}_n-\frac{L_{\theta ,n}-L_{\theta ,\text{g},n}}{t_{\text{s},n}}+\frac{\cos ({\theta }^{\prime })}{\sin ({\theta }^{\prime })}{\left(\frac{L_{\varphi ,n}}{R}\right)}^2\hfill \\ \hfill & +(-{\left(\frac{\partial \Phi }{\partial \theta }\right)}_{n+1}-\frac{L_{\theta ,n}-L_{\theta ,\text{g},n+1}}{t_{\text{s},n+1}}+\frac{\cos ({\theta }^{\prime })}{\sin ({\theta }^{\prime })}{\left(\frac{L_{{\varphi }^{\prime }}}{R^{\prime }}\right)}^2)\left(1+\frac{\text{d}t}{t_{\text{s},n}}\right)],\hfill \\ L_{\text{r},n+1}\hfill & =L_{\text{r},n}+\frac{\text{d}t/2}{1+\text{d}t\left(\frac{1}{2t_{\text{s},n}}+\frac{1}{2t_{\text{s},n+1}}+\frac{\text{d}t}{2t_{\text{s},n}{t}_{\text{s},n+1}}\right)}\cdot [-{\left(\frac{\partial \Phi }{\partial r}\right)}_{\text{n}}-\frac{L_{\text{r},n}-L_{\text{r},\text{g},n}}{t_{\text{s},n}}+\frac{1}{r_n}\left(\frac{L_{\varphi ,n}^2}{R_n^2}+\frac{L_{\theta ,n}^2}{r_n^2}\right)\hfill \\ \hfill & +(-{\left(\frac{\partial \Phi }{\partial r}\right)}_{\text{n+1}}-\frac{L_{\text{r},n}-L_{\text{r},\text{g},n+1}}{t_{\text{s},n+1}}+\frac{1}{r^{\prime }}\left(\frac{L_{\varphi ,n+1^{\prime }}^2}{{R^{\prime }}^2}+\frac{L_{\theta ,n+1}^2}{{r^{\prime }}^2}\right))\left(1+\frac{\text{d}t}{t_{\text{s},n}}\right)],\hfill \end{array}$$

a turbulent kick for the laminar disc case $$\begin{array}{llllll}{r}_n\hfill & =r_n+\delta r_{\text{d},\text{T}},\hfill & {\theta }_n\hfill & ={\theta }_n+\delta {\theta }_{\text{d},\text{T}},\hfill & {\varphi }_n\hfill & ={\varphi }_n+\delta {\varphi }_{\text{d},\text{T}},\hfill \end{array}$$

and finally a corrector step for the particle positions: $$\begin{array}{llllll}{r}_{n+1}\hfill & =r_n+\frac{1}{2}(L_{\text{r},n}+L_{\text{r},n+1})\text{d}t,\hfill & {\theta }_{n+1}\hfill & ={\theta }_n+\frac{1}{2}\left(\frac{L_{\theta ,n}}{r_n^2}+\frac{L_{\theta ,n+1}}{r_{n+1}^2}\right)\text{d}t,\hfill & {\varphi }_{n+1}\hfill & ={\varphi }_n+\frac{1}{2}\left(\frac{L_{\varphi ,n}}{R_n^2}+\frac{L_{\varphi ,n+1}}{R_{n+1}^2}\right)\text{d}t.\hfill \end{array}$$

References

All Tables

All Figures

Fig. 1 Measured particle scale height, Hp, in units of the gas scale height, H, as a function of the particle Stokes number for the simulations of the turbulent VSI, and the laminar disc plus stochasistic kicks at R = 1.8rp, averaged between 150 and 200 orbital periods. Shown are the results of the runs used in this paper (labelled with blue and green crosses), and of Stoll & Kley (2016) and Fromang & Nelson (2009) indicated with light blue and black circles, respectively. Additionally, we overplot the fit of the VSI particle scale heights by Youdin & Lithwick (2007; see Eq. (26)).

In the text

	Fig. 2 Spatial distribution of the dust particles as a function of time for the different planetary masses and for three representative particle sizes. The Stokes numbers from left to right are 0.08, 1.23, and 67.2.
In the text

	Fig. 3 Histogram of the dust surface number density distribution as a function of radius after 200 planetary orbits for three representative dust sizes in the turbulent (blue) and laminar (green) case. The Stokes numbers from left to right are 0.08, 1.23, and 67.2.
In the text

	Fig. 4 Azimuthal gas speed in units of the Keplerian speed as a function of radius for the different planetary masses and models (VSI = red, alpha-disc = green). When the gas speed becomes super-Keplerian outside the planet location, the dust-filtration process occurs and the pebble isolation mass is reached.
In the text

	Fig. 5 Surface density distribution of the dust particles after 80 planetary orbits for the different planetary masses and for three representative particle sizes. The Stokes numbers from left to right are 7.79 × 10−5, 1.23, and 67.2.
In the text

	Fig. 6 Eccentricity distribution for 3 representative size particles in the VSI (blue) and alpha-disc (green) model at the end of the simulation. The first six inner and outer first-order MMRs are overlaid.
In the text

	Fig. 7 Inclination distribution for 3 representative size particles in the VSI (blue) and alpha-disc (green) model at the end of the simulation. The first six inner and outer first-order MMRs are overlaid.
In the text

	Fig. 8 Number of accreted particles per orbit over time for the four different planet masses. Shown is the total number summed over all size bins.
In the text

	Fig. 9 Accreted particles as a function of the particle stopping time, where their approximate size is also reported in the top x-axis, integrated over 50 planetary orbits (from t = 100 to 150). The VSI (solid line) and alpha disc model (dashed line) are compared.
In the text

	Fig. 10 To measure the accretion (and survival) efficiency, the evolution of particles initially within a ring outside the planet co-orbital region (defined as in Eq. (34)) is monitored.
In the text

	Fig. 11 Efficiency of accreted (red line) and survived particles (green line) as a function of Stokes number for different planetary masses. The laminar disc run is represented with solid lines, while the VSI run with dashed lines. The fit from Lambrechts & Johansen (2014) is overplotted for intermediate τs values with a dashed grey line (see Eq. (40)).
In the text

Fig. 12 Comparison between the average drift speed of the different dust size particles. These values are calculated at the end of the simulation at R = 1.8 with the analytical drift speed (black solid line), the migration speed of a 5 M⊕ (dotted line), and a 10 M⊕ (dashed line) mass planet. The error bars reflect the impact of radial diffusion for the various stopping times. The smaller particle sizes in the VSI simulation are migrating outwards rather then inwards due to the inverse meridional circulation (Stoll et al. 2017a).

In the text

	Fig. A.1 Ratio of the azimuthal gas velocity to the Keplerian value for a 2D disc with an embedded 10 M⊕ planet at different times after insertion of the planet. The simulation in the left panel uses a viscosity of α = 10−4 and the right α = 5 × 10−4.
In the text