© ESO 2018

1 Introduction

Protoplanetary discs consist of gas and approximately 1% dust grains. These grains collide and grow to millimeter and even centimeter sizes (Brauer et al. 2008; Güttler et al. 2010). These particles are often referred to as pebbles. Pebbles interact with the gas disc through gas drag and drift inwards (Weidenschilling 1977; Brauer et al. 2008). Pebbles can become concentrated in pressure bumps and through the streaming instability (see Johansen et al. 2014 for a review), leading to planetesimal formation by gravitational collapse of the filaments. Planetesimals formed by the streaming instability have characteristic sizes of 100 km (Johansen et al. 2015; Simon et al. 2015).

In classical planet formation models, the cores of the giant planets form through mutual collisions between these planetesimals (Pollack et al. 1996). However, to achieve a core mass high enough (≈10 M_E) to attract a gaseous envelope, a surface density of planetesimals of a few times the Minimum Mass Solar Nebula (MMSN) is needed. Additionally, the growth timescale increases steeply with orbital distance, making the formation of the ice giants in the solar system basically impossible to achieve with planetesimal accretion alone. Gravitational stirring of the planetesimals by a set of growing protoplanets decreases the growth rates even more (Levison et al. 2010).

In recent years, a new paradigm of solid accretion has emerged: pebble accretion (Johansen & Lacerda 2010; Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012). When a pebble enters the planetary Hill sphere, it is subject to gas drag, which robs the pebble of angular momentum, resulting in an inward drift of the pebble onto the planet. When the largest planetesimals have grown to a few hundred kilometers in size by accreting other planetesimals, rapid pebble accretion allows further growth to cores of ten Earth masses well within the lifetime of the protoplanetary disc (Ida et al. 2016; Visser & Ormel 2016; Johansen & Lambrechts 2017).

A growing planet opens a partial gap in the protoplanetary gas disc, which influences the motion of solids in the disc (Paardekooper & Mellema 2006; Rice et al. 2006). Pebble accretion stops when the gap carved by the planet generates a pressure maximum outside of its orbit, which stops the inward flux of pebbles (Morbidelli & Nesvorny 2012; Lambrechts et al. 2014). This is referred to as the pebble-isolation mass. As the influx of pebbles is stopped, the planet’s gas envelope loses its hydrostatic support, and the envelope can then contract to form a planet with an extensive gaseous atmosphere (Lambrechts et al. 2014). In the solar system, pebble isolation is a potential mechanism to explain the dichotomy between the ice and gas giants, with ice giants never reaching the pebble-isolation mass and hence not able to undergo gas accretion within the lifetime of the protoplanetary disc (Lambrechts et al. 2014; Venturini & Helled 2017; Frehlik & Murray-Clay 2017).

Lambrechts et al. (2014) used hydrodynamical simulations to infer a pebble-isolation mass given by (1)

where H∕r is the discs aspect ratio. Lambrechts et al. (2014) found a weak dependence on the viscosity parameter. This is not surprising, since a dependence of gap opening on viscosity was also reported by Crida et al. (2006). However, because the dependence on the viscosity was not explored in detail in Lambrechts et al. (2014), an explicit mapping of the pebble-isolation mass as a function of viscosity is of crucial importance. Additionally, the headwind felt by the particles depends on both the disc aspect ratio, as evident in Eq. (1) above, and on the radial, initially unperturbed, pressure gradient ∂ ln P∕∂lnr.

Here, we investigate the dependence of the pebble-isolation mass on all local disc parameters, namely the disc aspect ratio H∕r, the viscosity ν, and the pressure gradient of the disc ∂lnP∕∂lnr. In order to probe this large parameter space, we adopt 3D isothermal simulations executed with the FARGOCA code (Lega et al. 2014; Bitsch et al. 2014). As particles with different sizes are coupled in different ways to the gas disc, we additionally integrate the trajectories of single pebbles with various sizes (and therefore Stokes numbers) inthe gas disc to probe which particle sizes can be trapped in the pressure bump as a function of planet mass.

The paper is organised as follows. In Sect. 2 we present our hydrodynamical set-up and discuss the different parameters that influence the pebble-isolation mass, to which we provide a fit in absence of turbulent diffusion. In Sect. 3 we integrate the trajectories of single pebbles in discs with embedded planets and infer the pebble sizes that are trapped in the pressure bumps. We then discuss turbulent diffusion of dust particles through the pressure bump. We also present the fitting formula for the pebble-isolation mass including turbulent diffusion, which is useful for planet formation simulations involving pebble accretion, in Sect. 3.3. In Sect. 4 we show the influence of the new-found pebble-isolation mass on simulations of planet formation, where we compare our results to Bitsch et al. (2015b). We additionally discuss implications of our results in Sect. 5 and finally summarise in Sect. 6.

2 Hydrodynamic simulations

2.1 Simulation set-up

In order to simulate the 3D disc-planet interaction, we used the 3D hydrodynamical code FARGOCA (Lega et al. 2014; Bitsch et al. 2014) in a locally isothermal configuration, where the radial temperature profile remains fixed throughout the simulation. We used the locally isothermal configuration because it allows a fast probing of parameter space in α, H∕r, ∂ ln P∕∂lnr, and planetary masses, which is needed to constrain the pebble-isolation mass. Here α is related to the viscosity through ν = αH²Ω_K (Shakura & Sunyaev 1973), where Ω_K denotes the Keplerian rotation. Additionally, locally isothermal simulations allow an easier probing of the aspect ratio because it can be set as an input parameter in contrast to simulations with heating and cooling, where the aspect ratio is set ultimately by the opacity profile that determines the cooling. Nevertheless, in Sect. 2.7 we test the predictions made with the isothermal simulations against simulations with heating and cooling.

We used for our simulations a 3D grid in spherical coordinates (r, ϕ, θ) with 315, 720, and 32 grid cells. Our grid ranged from 0.4 to 2.5 in radius, where the planet is located at 1, and spanned the full azimuthal range. We used evanescent boundary conditions for the radial boundaries to damp out the spiral waves exerted by the planet in order to avoid disturbances due to reflections of the spiral waves caused by the planet. We simulated different values of the viscosity parameter α, the aspect ratio H∕r, the pressure gradient ∂ ln P∕∂lnr, and a range of planetary masses (5–120 M_E ).

The planetary potential was modelled with a cubic potential (Kley et al. 2009). We used a smoothing length of r_sm = 0.6r_H, where r_H denotes the planetary Hill radius. This is the same smoothing length as in Lambrechts et al. (2014). The smoothing length has no influence on our results because the pressure bump generated by the planet outside of its orbit lies well beyond the smoothed zones. We tested different smoothing lengths, for instance, r_sm = 0.8r_H and r_sm = 0.4r_H, and the pressurebump outside of the planetary orbit did not change compared to our standard smoothing length of r_sm = 0.6r_H.

2.2 Measurement of the pebble-isolation mass

Pebbles in protoplanetary discs are subject to radial drift that is due to the headwind they feel from the gas (Weidenschilling 1977; Brauer et al. 2008). The gas orbits at a slightly sub-Keplerian speed because of the force exerted by the radial pressure gradient in the protoplanetary disc. This velocity difference is expressed as (2)

where (3)

and P is the pressure in the protoplanetary disc. In isothermal discs the pressure is given by , with c_s being the isothermal sound speed c_s = HΩ_K and ρ_g the gas volume density.

If η is lower than 0, the azimuthal gas velocity becomes higher than the Keplerian velocity and thus the particles feel a net outwards acceleration that will stop the inward motion of pebbles. As the planet grows, it carves a (partial) gap in the gas distribution around it by pushing material away from its orbit. This will eventually accelerate the gas outside of the planetary orbit to super-Keplerian velocities (Lambrechts et al. 2014). We calculate in the following an azimuthally averaged value of η. The azimuthally averaged η quantity gives a good handle on the generation of the pressure bump (see Appendix A). The planetary mass at which the created pressure bump stops the radial inward flow of pebbles through radial drift is called the pebble-isolation mass without diffusion.

In Fig. 1 we display the η parameter outside of the planetary orbit (the planet is fixed at r = 1) for several planet masses in a disc with H∕r = 0.05, α = 0.001 and Σ _g ∝ r^−0.5, where Σ _g denotes the gas surface density of the disc. A negative value of η means that the gas velocity is super-Keplerian, and inwards-drifting pebbles are stopped and cannot reach the planet any more. The pressure gradient is calculated for an azimuthally averaged pressure in the protoplanetary disc. In this case, a pebble-isolation mass of 25 M_E is found, in rough agreement with Lambrechts et al. (2014).

Fig. 1 Pressure gradient parameter η as a functionof orbital distance from the planet for α = 0.001, H∕r = 0.05, and differentplanetary masses. The location of the planet is fixed at r = 1. A negative η parameter indicates the formation of a pressure bump and with it super-Keplerian gas velocities that stop the inwards drift of pebbles. Here a mass of about 25 Earth masses is needed to generate the pressure bump outside of the planetary orbit. Inside of the planetary orbit, η is always positive.

2.3 Dependence on viscosity and aspect ratio

In Fig. 2 we present the pebble-isolation mass determined by 3D hydrodynamical simulations with different viscosities and aspect ratios. We define the pebble-isolation mass as the planetary mass at which the pressure bump outside of the planetary orbit becomes large enough to turn η negative. The pebble-isolation mass increases with α and with H∕r, as predictedby Lambrechts et al. (2014). We show a fit for our obtained data that scales with the aspect ratio, similar to Lambrechtset al. (2014), (4)

Additionally, the fit includes a dependency on α, which varies the pebble-isolation mass by a factor of 2–3 between low and high α values. Our fit is therefore also a function of α in the following way: (5)

where α₃ = 0.001.

In contrast, the original formula from Lambrechts et al. (2014), which is at the base of the new, refined formula proposed here, is similar to the expression of the critical mass for the wake to shock in a disc (Goodman & Rafikov 2001; Rafikov 2002). This highlights the fact that opening a deep gap around the orbit of the planet and creating a small gap that just reverses the pressure gradient are not similar processes and obey different physics.

To emphasise this effect, we calculated the depth of the gap following the formula for giant planet gap depths by Crida & Morbidelli (2007) given as (6)

The parameter is given by Crida et al. (2006) as (7)

Here q is the star-to-planet mass ratio, r_H the planetary Hill radius, and the Reynolds number given by . In our simulations we have checked that the left term of Eq. (7) varies from 2 to 10 for planets that have reached the pebble-isolation mass. Thus Eq. (6) predicts a gap depth of between 15% and 60%, while our 3D simulations yield gap depths of only 10–20%. The partial gap opened at the pebble-isolation mass thus corresponds to a different regime compared to the gap opening mass for giant planets of Crida et al. (2006) and the gap depth of Crida & Morbidelli (2007). Hence, the expression of Crida & Morbidelli (2007) for the depth of the gap should not be extrapolated to the regime of large (see Eq. (7)), low planet masses, and shallow gaps, and cannot be used to estimate the pebble-isolation mass.

2.4 Global pressure gradient

By changing the background gradient in surface density, the global pressure gradient changes. We varied the background surface density gradient Σ_g ∝ r^s in the disc from s = 0.5 to s = −1.5 and determined the pebble-isolation mass for discs with H∕r = 0.05 and α = 0.001 with the same method as above. A global inversion of the gas surface density gradient does not already imply a pressure bump in 3D simulations, in contrast to 2D simulations, because a pressure gradient inversion in 3D isothermal discs can only be reached by an inversion of the volume density gradient, where ρ ∝ r^s−1 for a radially constant H∕r. The difference between 2D and 3D simulations regarding the pebble-isolation mass is discussed in more detail in Appendix B. In Fig. 3 we show the pebble-isolation mass as a function of the background value of ∂ ln P∕∂lnr of the unperturbed disc.

The dependence on the background gradient of surface density (and thus on the background value of ∂ ln P∕∂lnr) is not very strong. We approximated the dependency of the pebble-isolation mass on ∂ ln P∕∂lnr with the expression (8)

where the reference value ∂lnP∕∂lnr = −2.5 corresponds to an unperturbed disc with H∕r =const. and Σ_g ∝ r^−0.5.

Fig. 3 Pebble-isolation mass as a function of ∂lnP∕∂lnr of unperturbed discs with different surface density gradients. Steeper surface density slopes s result in more negative ∂lnP∕∂lnr values. All simulations have been performed for planets in discs with H∕r = 0.05, f = 0 and α = 0.001.

2.5 Flared discs

Using Σ_g ∝ r^s and H∕r ∝ r^f , one can derive a dependency of η on the orbital distance r (using ), (9)

Here H₀∕r₀ indicates the aspect ratio at r = 1. This equation indicates that η only varies radially in discs with non-constant H∕r. As the pressure bump generated by the planet is located about 2H outside of the planet position, one could imagine that a change in η with orbital distance might influence how the pressure bump is generated. We therefore tested the influence of the flaring index on the pebble-isolation mass in isothermal discs with α = 0.001, Σ _g ∝ r^−0.5 and H∕r = 0.05r^f, where f spans from − 0.42 to + 0.42. We did not find any dependence on the pebble-isolation mass in discs with different flaring index. The pebble isolation mass in this case is determined only by the local unperturbed ∂lnP∕∂lnr value, α and H∕r at the location of the planet, but not by the flaring of the disc itself.

2.6 Pebble isolation mass without diffusion

To summarise the results of Sects. 2.3–2.5, we find that the pebble-isolation mass without diffusion is given by (10)

where (11)

with α₃ = 0.001.

2.7 Radiative simulations

In reality, discs have complex radial temperature profiles and a non-isothermal vertical structure. Different heating sources (viscous heating, stellar heating) are balanced by radiative cooling, which can alter the disc structure quite severely compared to simple power laws (Bitsch et al. 2015a). In order to test the predictions of the pebble-isolation mass (Eq. (10)), we studied the pebble-isolation mass in discs with heating and cooling.

In the adiabatic (and radiative) case, the sound speed changes by a factor of compared to the isothermal sound speed. This leads to a difference in the scale height of the protoplanetary disc for the isothermal and adiabatic configuration, which are related in the following way: (12)

In the radiative configuration, we therefore used the adiabatic scale height to estimate the pebble-isolation mass in a radiative disc. The disc set-up was similar to before, where we now additionally included radiative cooling and viscous heating (with α =6 × 10⁻³) as described in Kley et al. (2009). At the planet location, H_{adi, pla} = 0.0414 and ∂ ln P∕∂lnr = −3.26 (the flaring index of the disc is f =−0.38 at this location), which leads to a pebble-isolation mass of 28.6 M_E according toEq. (10). The results of our 3D simulations in discs with heating and cooling are shown in Fig. 4. Clearly, Eq. (10) matches the 3D simulations of discs with heating and cooling well.

Fig. 4 Pressure gradient parameter η as a function of orbital distance in a disc set-up with heating and cooling for two different planetary masses. The planet is placed at r = 1. The pebble-isolation mass is reached at ≈28.3 ME in the simulations, in good agreement with Eq. (10).

2.8 Application of the new fitting formula

The pebble-isolation mass depends not only on the disc aspect ratio H∕r, but also on the viscosity and the radial pressure gradient of the protoplanetary disc (Eq. (11)). For a fixed H∕r, a change in α from 10⁻⁴ to 10⁻² increases the pebble-isolation mass by a factor of ≈3 (Fig. 2), while anincrease in ∂lnP∕∂lnr from − 3.5 to − 1.5 decreases the pebble-isolation mass by about ≈30%. Making use of Eq. (10), we calculated the pebble-isolation mass in a disc with H∕r = 0.05 for different values of α and ∂lnP∕∂lnr and show the resulting pebble-isolation mass in Fig. 5.

Clearly, high values of viscosity ( ≈10⁻²) increase thepebble-isolation mass significantly, where the pebble-isolation mass can reach over ≈50 M_E for our nominal ∂lnP∕∂lnr = −2.5. Nevertheless, the strongest dependence of the pebble-isolation mass is on the disc aspect ratio H∕r. The disc aspect ratio is determined by the heating of the disc, either through viscosity or stellar irradiation. As the disc evolves in time, the aspect ratio decreases in the inner part of the disc as a result of reduced viscous heating and in the outer parts as a result of a decreasing stellar luminosity (Bitsch et al. 2015a). These effects reduce the pebble-isolation mass in time and only discs with high viscosities can maintain a high pebble-isolation mass in the outer parts of the disc, as the disc evolves in time.

Fig. 5 Pebble-isolation mass as a functionof the pressure gradient ∂lnP∕∂lnr and α in a disc with a constant H∕r = 0.05. The two black lines mark 20 and 50 ME. Clearly, higher values of viscosity and ∂lnP∕∂lnr result in significantly higher pebble-isolation masses.

3 Drift of small pebbles through the bump

Small particles (τ_f ≪ 1) are stronglycoupled and move with the radial gas accretion flow, while larger particles (τ_f ≫ 1) are only weakly affected by gas drag. Here τ_f denotes the Stokes number of the pebbles. The acceleration of a pebble in a gas disc is given by (13)

where r denotes the vector between the central star and the pebble, and t_f is the friction time, which is related to the Stokes number with τ_f = t_fΩ_K and Δv = ηv_K. The variables v_peb and v_gas are the pebble and gas velocities, respectively.

In Fig. 6 we show the radial and azimuthally averaged gas velocities in locally isothermal discs with embedded planets. The pressure bump generated by the planet outside of its orbit is clearly visible in the azimuthal velocity pattern, where the gas can reach speeds higher then the Keplerian value. Particles entering this pressure bump can be trapped, depending on their size.

A negative radial velocity indicates an inward flow of the gas, while a positive radial velocity indicates an outward movement of the gas. The planet generates a radial outward flow of gas close to its vicinity, but limited to the region in front of the pressure bump. This outward flow is related to the gap-opening process, where the planet pushes the material away from its orbit. This material can then move upwards to maintain hydrostatic equilibrium and again falls in from the top regions of the disc onto the planet. The same meridional flow was also observed for gap-opening giant planets (Morbidelli et al. 2014). In this region, inflowing particles might in principle be trapped as well.

After the hydrodynamical simulations shown in Sect. 2 reached an equilibrium state, we integrated the movement of test particles in the steady-state gas distribution to determine the dependence of the pebble-isolation mass on the Stokes number of the particles. The steady-state surface density profile and η profile of the disc are shown in Appendix A, where the pressure bump is clearly visible for a 25 M_E planet.

We integrated the pebble trajectories in 2D planes of the protoplanetary disc, where we mainly focused on integration in the disc midplane. Pebbles injected at higher altitudes in the disc are additionally subject to vertical settling, which moves them quickly towards the midplane. Integrating the pebble trajectories in a 2D plane above the disc midplane revealed, however, that pebbles are stopped at all altitudes when the planet has reached isolation mass. This implies that at higher altitudes pebbles cannot drift through the generated pressure bump either.

In Fig. 7 we show the trajectories of pebbles in the midplane with constant τ_f = 1.0 in the gas velocity field generated by planets with 10 M_E and 25 M_E. As the 10 M_E planet does not generate a pressure bump outside of its orbit (Fig. 1), the pebbles drift through towards the inner disc. The 25 M_E planet generates a pressure bump in the disc (Fig. 1) where the pebbles can be trapped.

In the top panel of Fig. 8, we show the time evolution of the orbital distance of integrated pebble trajectories with different Stokes number τ_f in a disc with an embedded 25 M_E planet. Pebbles with τ_f > 0.005 are trappedin the pressure bump and do not drift inwards any more. Particles with τ_f < 0.005 are well enough coupled with the gas to move through the pressure bump towards the system interior to the 25 M_E planet. The pressure bump generated by the 25 M_E planet is even quite weak (Fig. 1), which explains why small particles can still drift through.

The pressure bump generated by the planet increases in strength with planetary mass (Fig. 1), which allows the trapping of pebbles with smaller τ_f in it for higher planetary masses (bottom panel in Fig. 8 for 30 M_E planet). Increasing the planetary mass by 20% reduces the Stokes number of particles that can still drift through the generated pressure bump by more than an order of magnitude. Now particles with τ_f > 5 × 10⁻⁴ are trapped inside the pressure bump and can no longer reach the inner system.

Very small particles, however, with τ_f ≈ 1 ×10⁻⁴, are so strongly coupled with the gas that they completely follow the gas flow (Brauer et al. 2008) and are thus no longer blocked by the pressure bump located at r = 1.13 (see Sect. 3.1). Instead, the pebbles drift through the pressure bump, but are then caught just outside of the planetary orbit because of the radial outward flow of the gas (Fig. 6), which prevents further inward drift. For the 25 M_E planet, the pebbles also drift through the radial outward flow of the gas because it was not strong enough to keep the pebbles from the inner disc (see Sect. 3.1).

Fig. 6 Midplane gas velocities as a function of orbital distance in discs with H∕r = 0.05, α = 0.001, and Σ ∝ r−0.5. Top plot: radial velocity in units of the sound speed at r = 1.0. A negative velocity indicates an inward flow of the gas. The red line indicates the radial velocity of the unperturbed disc. Bottom plot: azimuthal velocity in units of the Keplerian velocity. If vθ is larger than 1, the gas orbits super-Keplerian, indicating the pressure bump in the disc.

Fig. 7 Trajectories of pebbles with τf = 1.0 that started at (1.5;0) moving in the gas velocity field generated by planets with MP = 10 ME and MP = 25 ME . We plot the pebble trajectories in the rotating frame, and the planet position is marked by the black dot. The 10 ME planet does not generate a pressure bump, so the pebble drifts through all the way towards the inner boundary of the computational domain. The 25 ME planet, on the other hand, generates a pressure bump outside of its orbit and thus blocks the flow of pebbles.

Fig. 8 Evolution of the orbital distance of pebbles as they drift through a disc with an embedded 25 ME planet (top) and an embedded 30 ME planet (bottom). Even though the planet has generated a pressure bump outside of its orbit (maximum at r = 1.123, see Fig. 1), small pebbles with τf < 0.005 can drift through the pressure bump for the 25 ME planet, while the Stokes number of the particles has to be an order of magnitude smaller to drift through the pressure bump generated by the 30 ME planet.

3.1 Dependence on the radial gas velocity

The radial velocity of the gas in an α accretion disc is determined directly by the viscosity of the disc, (14)

However, there is a strong debate in the literature about the causes of the turbulence and about the size of its magnitude (Turner et al. 2014). As the drift velocity of the particles depends on the gas velocity, we investigate in this section how a change in the radial gas velocity influences particle drift through the disc. We are particularly interested in how a change in the radial gas velocity allows or hinders particles from drifting through the pressure bump generated by the azimuthal gas velocity changes induced by the planet outside of its orbit (Fig. 6). We artificially modified the radial velocity pattern to a fixed value, but kept the azimuthal gas velocity profile of a disc perturbed by a planet (Fig. 6). In this way, we mimicked the effects of different levels of turbulence without simulating discs with magneto-rotational instability (MRI) or vertical-shear instability turbulence.

In Fig. 9 we show the trajectories of pebbles embedded in discs with fixed radial gas velocities, but with azimuthal gas velocity profiles that correspond to Fig. 6 for the 25 M_E planet. In the top panel we show the trajectories of pebbles in discs with radial gas velocities lower than in Fig. 8, while in the bottom panel we show the trajectories of pebbles in discs with radial gas velocities higher than in Fig. 8. Clearly, a lower radial gas velocity allows for more efficient trapping of smaller pebbles compared to a higher radial gas velocity.

The reason is that the radial pressure gradient is negative everywhere in the disc except at the centre of the outer pressure bump generated by the planet, where it is zero for a planet that has just reached pebble-isolation mass. Hence any tiny radial gas velocity can transport dust particles of any size across the pressure bump. Pebbles are only safe when the planet is more massive and the inner edge of the pressure bump has a positive radial gradient in pressure. As the radial gas velocity is determined by the viscosity in an α-accretion disc, the movement of particles is determined by viscosity for discs with high viscosity and by drift in discs with low viscosity (de Juan Ovelar et al. 2016).

Pebbles with a (positive) terminal velocity (15)

high enough to compensate for the (negative) radial gas velocity v_r,g can be trapped in the pressure bump. This can be reformulated as (16)

Using this equation, we can estimate the Stokes number of particles that are blocked by the pressure bump generated by the planet by just looking at the velocity fields. The radial gas velocity is ≈1 cm s⁻¹ in our disc model, and the measured Δv is shown in Fig. 10. This indicates that pebbles with τ_f > 0.005 should be blocked by a planet of 25 M_E and particles with τ_f > 0.0003 should be blocked by a planet of 30 M_E, in agreement with our simulations (Fig. 8). We have marked the minimal Stokes number of particles that can be blocked by the pressure bump generated by the planet with blue circles in Fig. 10 for a radial velocity of 1 cm s⁻¹ . The lowest Stokes number depends linearly on the radial gas velocity, which is slightly different for the higher planetary masses (Fig. 6). For a 50 M_E planet, the radial gas velocity is roughly ≈2 cm s⁻¹ (Fig. 6), higher than for the 25 M_E planet, because the planet influences the velocity pattern of the disc. This implies that pebbles with τ_f > 1.2 × 10⁻⁴ can be stopped by a 50 M_E planet, but we note that the blue dots in Fig. 10 correspond to v_g = 1 cm s⁻¹. This is also in agreement with our simulations. These results do not take diffusion into account, which we discuss in the next section.

Fig. 9 Evolution of the orbital distance of pebbles as they drift through a disc with an embedded 25 ME planet, where the radial velocity is fixed to − 1 × 10−5cs,0 (top) and − 3 × 10−5cs,0 (bottom). As the drift speed of the particles depends on the gas velocities, particles of different sizes cannot be blocked or drift through compared to the nominal case shown in Fig. 8. In particular, a slower gas flow allows for a more efficient particle trapping in the pressure bump generated by the planet.

Fig. 10 Velocity perturbation Δv as a function of planet mass at the location of the pressure bump. A negative Δ v value indicates that the gas speed is super-Keplerian. The blue circles mark the lowest Stokes number of particles that can be stopped at the pressure bump assuming a radial gas velocity of 1 cm s−1 and following Eq. (16).

3.2 Diffusion of dust particles

Our hydrodynamical simulations do not include any turbulent motion of the gas, such as those seen in magnetohydrodynamics (MHD) simulations (Bai 2017) or in simulations with the vertical shear instability (Nelson et al. 2013; Stoll & Kley 2016). These turbulent motions in the gas velocities can act on the movements of the pebbles, giving them random kicks. Several authors have considered the effects of diffusion on particles embedded in discs with planets. Paardekooper & Mellema (2006) included diffusion into the motion of dust particles in gas discs in the presence of 30 M_E planets and estimated this effect to be of the order of 1%, indicating that diffusion of dust particles does not play a role in opening a gap in dust distribution of protoplanetary discs. Pinilla et al. (2016), on the other hand, studied dust filtration by giant planets in the context of transition discs. Giant planets open deep gaps in protoplanetary discs that prevent dust from drifting through. However, in their 2D simulations, the authors found that a planet of 1 M_Jup still does not stop all dust particles and small dust grains (τ_f ≈ 10⁻³ ) can drift through the gap generated by the planet at 20 AU. Their disc set-up would lead to a pebble-isolation mass of 27.5 M_E according toEq. (10), where the pebble-isolation mass is reached in 2D disc simulations at lower masses (see Appendix B) than in 3D simulations.

Pinilla et al. (2016) considered turbulent mixing of dust particles, where the dust diffusivity follows the prescriptions by Youdin & Lithwick (2007), which depend on the Stokes number and the gas diffusivity (assumed to be equal to the disc viscosity). If then the pressure gradient is not steep enough, the particles can be released from the pressure bump, where particles with τ_f > α are trapped. Smaller particles are diffused out of the pressure bump and dragged by the gas. Without diffusion, the trapping of particles is more efficient (Pinilla et al. 2012). This mechanism allowed the small particles to move across the pressure bump generated by the planet in Pinilla et al. (2016), while our simulations show an effective trapping of small particles as a result of the lack of diffusion. The difference of Pinilla et al. (2016) to Paardekooper & Mellema (2006) is probably related to different prescriptions of diffusion. We therefore estimate the effects of diffusion in the following for the pebble-isolation mass.

The equilibrium between radial advection of dust particles and turbulent diffusion is achieved when (as also stated in Pinilla et al. 2012) (17)

where v_r,p is the radial velocity of the pebbles, ρ_p the pebble density, ρ_g the gas density, D the diffusion coefficient parameterised by D = αc_sH_g , and ϵ = ρ_p∕ρ_g is the dust-to-gas ratio. The radial velocity of the pebbles is given by (18)

To efficiently trap dust in a pressure bump despite turbulent diffusion, the equilibrium flux stated above must be obtained for a radial dust scale-length (H_p ) that is roughly equal to the extent of thepressure bump (H_b). As the latter has an approximate width of one gas scale-height, this leads to (19)

We estimatethe scale-height of the dust in the pressure bump from Eq. (17) to obtain (20)

Expanding now D, we arrive at (21)

which leads to (22)

where Π = Δv∕c_s. Particles with Stokes numbers larger than this (Eq. (22)) are still trapped in a pressure bump, particles smaller than this can diffuse through the pressure bump. We show this critical Stokes number in Fig. 11 for a disc with H∕r = 0.05.

We show Δv∕c_s as a function of planetary mass in a disc with α = 0.001 and H∕r = 0.05 in Fig. 12. We can now compare directly to Fig. 11 to determine when and how pebbles inside the pressure bump are affected by diffusion. For a 30 M_E planet in a disc with α = 0.001, the minimum Stokes number of particles that are not affected by diffusion is ≈0.025, while for the 50 M_E planet, this particle size is ≈0.004, which corresponds to the particle size stopped by the 25 M_E planet in the case of no diffusion (Fig. 8). This indicates that the pebble-isolation mass for these Stokes numbers in case of diffusion is a factor of ≈2 higher than in the case withoutdiffusion for H∕r = 0.05 and α = 0.001, depending on τ_f.

Inside the pressure bump generated by a Jupiter-mass planet in 2D simulations of a disc with α =0.001 and H∕r = 0.05, the maximum Δ v∕c_s = 0.66. This leads to a critical Stokes number of τ_f = 7.5 × 10⁻⁴, in agreement with previous studies including diffusion, see for example Pinilla et al. (2016), where very small particles can drift through the pressure bump generated by large planets. Adding the effects of particle diffusion in 2D disc simulations, Ataiee et al. (2018) found that diffusion can increase the pebble-isolation mass, in agreement with Pinilla et al. (2012) and our estimates.

The difference between the pebble-isolation mass derived from pure hydrodynamical simulations compared to simulations taking diffusion into account also depends on the level of turbulence in protoplanetary discs. For example, blocking particles with τ_f = 0.001 requires for α =2 × 10⁻⁴ an increase in pebble-isolation mass by about a factor of 1.5, while for α = 0.001 an increase of much more than a factor of 2 is needed (a Jupiter-mass planets blocks pebbles with τ_f > 7.5 × 10⁻⁴ in the case of diffusion).

However, the level of turbulence in discs is not very well constrained, where α values from 10⁻⁴ to 10⁻¹ can be reached in simulations of the MRI (see Turner et al. 2014 for a review). Simulations with hydrodynamical instabilities in regions of the disc that are not subject to MRI-driven turbulence show α values of a few times 10⁻⁴, which probably sets a lower limit on turbulence (Klahr & Bodenheimer 2003; Nelson et al. 2013; Stoll & Kley 2016). Recent simulations indicate that disc winds could be the main driver of disc accretion. In these simulations, the midplane regions only have a very low level of turbulence. For such a low level of turbulence, the effect of turbulent diffusion on the pebble-isolation mass is quite low compared to discs with high viscosity.

The dominant particle size in the planet formation simulations of Bitsch et al. (2015b) with α =0.0054 presented below is of the order of 0.1 (Fig. 13), and about 0.01–0.2 when taking fragmentation into account (see below). From Fig. 11 it can be inferred that taking turbulent diffusion into account, the pebble-isolation mass for these Stokes numbers is higher than predicted in the case without diffusion by about a factor of 2.

Fig. 11 Critical Stokes number of particles that are affected by turbulent diffusion as a function of absolute Δ v∕cs in the pressure bump and α. Particles with Stokes numbers lower than the critical value easily diffuse and can drift through a pressure bump. The black lines denote particles with Stokes numbers of 10, 1, 0.1, 0.01, and 0.001 from left to right. The white dots correspond to different planetary masses, where we have taken the maximum value of the pressure bump generated by the planet for simulations with different α. The planetary masses are marked with the white numbers next to the dots. The horizontal dashed black line markes our nominal α = 0.001.

Fig. 12 Δv∕cs as function of orbital distance in discs with α = 0.001 and H∕r = 0.05 and embedded planets of different masses. Note that Δv∕cs = η∕(H∕r). The Δ v∕cs values shown here correspond to the white dots at α = 0.001 (dashed line) in Fig. 11.

Fig. 13 Stokes number of the pebbles in the simulations presented in Figs. 15 and 16 at a given time t and orbital distance r in the disc. The black lines mark Stokes numbers of 0.05 to 0.4. Note that each planet growth trajectory stars at a given point in r-t and the planet the moves to higher time values, indicating that the Stokes number of the pebbles accreted by the planet increases. Additionally, the planet migrates in the disc to smaller orbital distances, which also increase the Stokes number of the accreted particles. The over-plotted blue lines correspond to the growth tracks shown in Fig. 14, where the solid line corresponds to solid (pebble) accretion, the dot marks the pebble-isolation mass, and the dashed line corresponds to the gas accretion phase. The planetary growth track moves upwards to increasing time.

3.3 Pebble-isolation mass including diffusion

The pebble-isolation mass does not only depend on H∕r, α and ∂ ln P∕∂lnr, but also, as shown in the previous subsection, on the turbulent diffusion of particles. The critical pressure gradient parameter Π_crit to block particles of Stokes number τ_f is given by (23)

From our hydrodynamical simulations we can measure how Π = Δv∕c_s in the pressure bump generated by a planet that has already reached changes with planetary mass as (24)

where (25)

where f_fit is defined in Eq. (11). This fit only applies to planets that have already reached the pebble isolation mass without diffusion , because λ gives the slope of the change of Π inside the pressure bump generated by the planet, where is the minimum mass needed to invert the radial pressure gradient ∂ ln P∕∂lnr in the disc. We note that λ is only valid until M_p ≈ 2.5 , when λ changes, because the growing planet slowly transitions into the gap depth regime predicted by Crida & Morbidelli (2007). When setting Π _crit = Π, we can define the pebble isolation mass with diffusion M_iso as (26)

Using Eq. (26), we study in the next section the effect of this new-found pebble isolation mass on the formation of planetary systems and the core masses of the formed planets.

4 Influence on planet formation

The pebble isolation mass determines the final mass of the planetary core in the pebble accretion scenario because the reduced accretion luminosity facilitates the accretion of gas (Lambrechts et al. 2014; Lambrechts & Lega 2017) and the planet can eventually grow to become a gas giant. The formation pathway of the growing planet is determined by the growth rate and size of the planetary core because it influences its gas accretion rates (Piso et al. 2015) and migration behaviour (Baruteau et al. 2014). By reaching a different pebble isolation mass, the planet can undergo a different formation pathway.

We therefore investigate in this section the influence of the pebble isolation mass on planet growth by comparing planet growth simulations to the new pebble isolation mass (Eq. (26)) with simulations with the pebble isolation mass measured by Lambrechts et al. (2014), who only inferred the dependence on H∕r. For this we make use of the planet growth simulations including planet migration and disc evolution presented in Bitsch et al. (2015b).

4.1 Planet growth and migration model

The planet growth and migration model is described in great detail in Bitsch et al. (2015b), therefore we only repeat the essential points here. The planet growth and migration rates strongly depend on the disc structure. We used here the disc structure model of Bitsch et al. (2015a). This semi-analytical disc model features bumps and dips in the inner disc structure caused by transitions in the opacity at the water ice line, which can act as planet traps for low-mass planets (Bitsch et al. 2015a; Bitsch & Johansen 2016) and evolves in time. We used a disc lifetime of 3 Myr.

The growth rate of the planet depends on the pebble surface density Σ _peb at the planet location of the protoplanetary disc (27)

where r_H is the planetary Hill radius and v_H the Hill speed at which the particles enter, given by v_H = r_HΩ_K. In the drift-limited growth of dust particles to pebbles (Birnstiel et al. 2012), the pebble surface density depends on the pebble flux Ṁ_peb (Lambrechts & Johansen 2014) in the following way: (28)

Here ϵ_P = 0.5 (Lambrechts& Johansen 2014). We note that the nominal pebble flux used in Bitsch et al. (2015b) was overestimated by a factor of ≈10 and that we used here a modified pebble growth model presented in Bitsch et al. (2018).

After the planet has reached its pebble-isolation mass, it can contract a gaseous envelope (Piso & Youdin 2014), and as soon as the mass of the gaseous envelope is higher than the core mass of the planet, runaway gas accretion can start (Machida et al. 2010).

Growing planets interact with their natal protoplanetary disc and migrate in it. Low-mass planets do not perturb the disc significantly and migrate in type I migration, which depends mainly on the disc viscosity and on the radial gradients of surface density, temperature, and entropy (Paardekooper et al. 2011). The disc structure is therefore of crucial importance in determining the migration rates. Planets growing further (e.g. as a result of rapid gas accretion) push the gas away from their orbit (or accrete it, Crida & Bitsch 2017) and open a gap in the protoplanetary disc and migrate in type II migration. This migration rate depends on the viscosity of the protoplanetary disc. For a review on planet migration, see for example Baruteau et al. (2014).

To calculate the torque Γ exerted by the disc on the planet, we followed Paardekooper et al. (2011) for type I migration and the viscous evolution for type II migration. The orbital migration time t_m is given as (29)

where J is the angular momentum of the protoplanet. Migration here is thus defined in terms of the total torque exerted on the orbit. The migration time so defined is positive when the total torque is negative. For constant eccentricity, the time taken to migrate to the centre is t_m ∕2 (Papaloizou & Larwood 2000). This factor of 2 was absent in the original simulations of Bitsch et al. (2015b) and was now added here.

4.2 Stokes numbers in the planet formation model

The Stokes numbers of the pebbles in our model are shown in Fig. 13 as a function of distance and time. The blue lines indicate the growth tracks of the planets shown in Fig. 14 as the planets grow and migrate. The solid lines correspond to solid accretion, which stops when the planet reaches pebble-isolation mass (marked as a dot), while the dashed lines indicates gas accretion.

Typically, the pebbles accreted by planets in our model have Stokes numbers in the range of 0.05 < τ_f < 0.4, as those are the Stokes numbers of the pebbles dominated by radial drift (Birnstiel et al. 2012; Lambrechts & Johansen 2014). As shown in Sect. 3.2, only small pebbles can drift through the pressure bump (with τ_f < 0.01 when the planet just reached pebble-isolation mass in the absence of diffusion, ). In our model, α = 0.0054, indicating a strong turbulent diffusion efficiency (Fig. 11), but even turbulent diffusion cannot carry the particles of the dominant Stokes number across the pressure bump because of the large pebble sizes (Fig. 13). Therefore the pebble-isolation mass is only slightly increased by the effects of turbulent diffusion in our model. This means that the bulk of the pebbles and therefore also the bulk of the solid mass is blocked by a planet that has reached the pebble-isolation mass if the particle size is dominated by radial drift. Additionally, the blocking of pebbles becomes more efficient as the planet grows and starts to accrete a gaseous envelope.

The final pebble sizes in our model were determined by radial drift alone, where we did not take the effects of fragmentation (Brauer et al. 2008; Birnstiel et al. 2011), bouncing (Zsom et al. 2010), or condensation (Ros & Johansen 2013; Schoonenberg & Ormel 2017) into account. Fragmentation and bouncing can lead to smaller pebble sizes than in the drift-limited case, while condensation of volatiles onto already existing pebbles can increase their size.

The sizes of pebbles determined by the fragmentation limit are given in Birnstiel et al. (2015) as (30)

where ρ_peb is the density of the pebble itself (set to 1.5 g cm⁻³ for water ice) and v_f is the fragmentation velocity limit, where water-ice particles have a higher fragmentation velocity of 10 m s⁻¹ than silicate grains (Gundlach & Blum 2015). When we use only this fragmentation limit for icy particles, the Stokes numbers of the particles in our disc model are 0.01–0.2, which is slightly smaller than in the drift-limited case (shown in Fig. 13). Lower Stokes numbers will result in higher pebble-isolation masses (Eq. (26)) and thus higher core masses. However, condensation at ice lines could lead to even larger particles in these regions. Future models of planet formation have to take these effects into account in order to calculate more realistic grain sizes.

In the following, we use Eq. (26) to include the effects of diffusion in the calculations of the pebble-isolation mass and to simulate the growth of planets through pebble accretion, where we use the Stokes numbers of the drift-limited solution shown in Fig. 13.

Fig. 14 Growth tracks of planets starting at several different initial positions in a disc that is already 1.5 Myr old for pebble-isolation masses given by Lambrechts et al. (2014) in Eq. (1) and for the new-found pebble isolation mass Miso in Eq. (26). The initial growth phase is the same for both isolation masses, so the growth tracks diverge only when the pebble-isolation mass (Eq. (1)) is reached, marked by the dots. The solid lines indicate pebble accretion, while the dashed lines mark gas accretion.

4.3 Growth tracks

Growth and migration depend on the structure of the protoplanetary disc, where we follow the disc model of Bitsch et al. (2015a), and use a dust metallicity of Z_dust = 0.5% to set the disc opacity. The disc viscosity is α = 0.0054. The planetary growth rate depends crucially on the amount of available pebbles that can be accreted by the planet. In the remainder of the paper, we use Z_peb = 1.0% as in Bitsch et al. (2015b). We set the disc lifetime to be 3 Myr.

In Fig. 14 we show the growth tracks of planetary seeds starting at different locations in a disc that is already 1.5 Myr old, meaning that the planets evolve for 1.5 Myr to reach a disc lifetime of 3 Myr. For each orbital distance we ran two simulations, where the only difference was the final pebble-isolation mass, determined either by Lambrechts et al. (2014), Eq. (1), or by the new-found pebble-isolation mass, Eq. (26). This means that the initial growth is the same for the two simulations, until the planet in one simulation reaches the pebble-isolation mass and gas accretion starts. This is generally the case for simulations following the pebble-isolation mass of Eq. (1), which is generally lower than Eq. (26), especially for the given disc model with α = 0.0054.

The final core mass also determines the contraction phase of the envelope, where (Piso & Youdin 2014), which in turn determines how fast the planet transitions into runaway gas accretion and can then open a gap and transition into the slower type II migration phase.

In the inner parts of the protoplanetary disc, the pebble-isolation mass is low because of the low H∕r, while the pebble-isolation mass is high in the outer parts of the disc because H∕r is high. For planets forming in the inner regions of the disc (r < 15 AU), the difference between Eq. (1) and Eq. (26) is not that large, allowing planets to arrive at similar total masses and orbital distances. However, the core masses of planets formed using Eq. (26) are higher. In the outer disc, the pebble flux is quite low, so that the differences in the pebble-isolation mass result in a slower growth of the planets, where M_iso is determined by Eq. (26), because pebble accretion is slower than gas contraction for these pebble densities. This also results in further inward migration before the planet opens a deep gap and transitions into type II migration. In total, the differences regarding the final orbital position and the final planetary mass seem not very great. However, the differences in the core masses can be up to 30% (see Fig. 16), which is crucial for the formation of the ice giants in our solar system, which had low core masses in Bitsch et al. (2015b).

Fig. 16 Final core masses of planets Mc as a function of formation distance r0 and formation time t0 in the disc. The white lines correspond to core masses of 10, 15, and 20 ME (top to bottom) and are marked by white numbers. The black and blue lines have the same meaning as in Fig. 15. The top plot corresponds to the pebble-isolation mass stated in Lambrechts et al. (2014), while the bottom plot corresponds to the pebble-isolation mass in this work (Eq. (10)). Clearly, Eq. (10) allows for higher core masses, which is crucial for the formation of ice giants in our own solar system.

4.4 Global picture

We now extend the approach of the growth tracks to probe the planetary growth for starting positions of the planetary seeds from r₀ = 0.2 AU to 50 AU and from t₀ = 100 kyr to 3 Myr. In Fig. 15 we show the final total planetary mass of planets as a function of r₀ and t₀ for pebble-isolation masses following Eq. (1) (top) and Eq. (26) (bottom). The white crosses mark the growth tracks shown in Fig. 14.

At first glance, the difference in final orbital position r_f and final planetary mass M_p is not that large compared for the different pebble- isolation masses, in agreement with Fig. 14. For 0.1 Myr < t₀ < 2.0 Myr and r < 20 AU, the formation of close-in planets that have reached the inner edge of the disc at 0.1 AU is triggered. These planets form too close to the central star, so that planetary migration drives them towards the inner edge of the disc during the lifetime of the protoplanetary disc for our migration rates in discs with high viscosity. As the pebble-isolation mass is higher, the planetary cores with low-mass gaseous envelopes (that have not reached runaway gas accretion) become too large to be contained in the region of outward migration (which can only hold planets of a few Earth masses after about 1 Myr; Bitsch et al. 2015a). These planets then drift inwards as rock-dominated planets (bottom panel of Fig. 15 in contrast to the top panel of Fig. 15).

Planets forming in the outer part of the protoplanetary disc reach higher pebble-isolation masses owing to the higher aspect ratio. In the top panel of Fig. 15, the pebble-isolation mass is reached earlier, but the core masses are high enough ( ≈10 M_E ) to allow a transition into runaway gas accretion and thus gas giant formation. However, in the bottom panel of Fig. 15, the pebble-isolation mass is higher, which prolongs the core formation timescale and thus results in gas accretion at later stages. The overall differences in r_f and M_p are not verylarge, however.

Clearer differences can be seen with respect to the core masses of these planets (Fig. 16), where Eq. (26) delivers core masses that are about ≈30% higher, compared to Eq. (1) in our model with α = 0.0054. Figure 2 shows the clear dependence on H∕r for the pebble-isolation mass, explaining why the core masses increase with orbital distance of the formed planets, because H∕r increases outwards in a stellar irradiated disc (Bitsch et al. 2013, 2015a). However, the 10 M_E core mass line seems to be constant as a function of t₀ at early times in the inner regions of the disc (r < 15 AU) of Fig. 16: the growth time in the outer parts of the disc is much longer because of the lower pebble flux and the larger pebble scale height. The planetary growth additionally competes with planetary migration (driving the planet inwards to parts of the disc with lower H∕r and thus lower pebble-isolation mass) and disc evolution, where the disc aspect ratio decreases with time (Bitsch et al. 2015a).

These higher core masses are more consistent with the planetary structure of the solar system¹ , where Uranus and Neptune have not reached the pebble-isolation mass (Lambrechts et al. 2014) and thus stayed at 15–20 M_E without accreting a large gaseous envelope. The final core mass could additionally be increased if more pebbles are available, allowing a faster growth of the planets.

Fig. 15 Final masses of planets (total mass MP = Mc + Menv) as a function of formation distance r0 and formation time t0 in the disc. Planets that are below the dark blue line have reached pebble-isolation mass and can accrete gas. All planets that are below the white line have Mc < Menv, indicating that they have undergone runaway gas accretion. The black lines indicate the final orbital distance rf of the planet, namely 0.1, 0.5, 1.0, 5.0, 10.0, and 20.0 AU. The top plot corresponds to the pebble-isolation mass found in Lambrechts et al. (2014), while the bottom plot corresponds to the pebble-isolation mass in this work (Eq. (10)). The higher pebble-isolation mass suppresses the formation of gas giants in the very outer parts of the disc at late times.

Fig. 2 Pebble-isolation mass as a function of α and for different aspect ratios H∕r. The pebble-isolation mass is fitted through different α values with asimple fit that also scales with .

5 Discussion

5.1 Planet migration versus pebble-drift speeds

Planets embedded in protoplanetary discs interact gravitationally with the disc and move through the disc (Paardekooper & Mellema 2006; Baruteau & Masset 2008; Kley et al. 2009; Paardekooper et al. 2011). The migration timescale τ_mig of low-mass planets is estimated in Tanaka et al. (2002) and given as (31)

Here, r_pl and are the orbital distance and the aspect ratio at the planet location. The constant C reflects the migration speed through the disc surface density profile and disc temperature profile, given by C = 1∕(2.5 + 1.7β_T − 0.1α_Σ) (Paardekooper et al. 2011), where α_Σ is given by Σ_g = Σ₀r^−α_Σ with Σ₀ = 350 g/cm² and β_T by T ∝ r^−β_T. For our standard disc with α_Σ = 0.5 and β_T = 1, the pre-factor is C = 0.24.

We can now compare this with the radial drift speed of particles (Brauer et al. 2008), which is given by (32)

The radial speed of the gas v_r,gas in an α disc is estimated by Takeuchi & Lin (2002) as (33)

The quantity v_r,gas that describes the radial drift of individual dust particles is given by Weidenschilling (1977) as (34)

where Δv is the maximum drift velocity, which is calculated as (35)

In Fig. 17 we show the radial drift speed of particles as a function of Stokes number, the radial gas velocity, and the migration speed of planets with 25 M_E (corresponding directly to the pebble-isolation mass without diffusion) and 50 M_E in a disc with α = 0.001. Clearly, the particles drift faster than the planet migrates when the pebble-isolation mass is reached. Even for the 50 M_E planet, particles with τ_f > 0.001 drift faster than the planet migrates, indicating that particles drifting inwards from the outer disc will be trapped in the pressure bump outside of the planetary orbit.

However, planetary migration is more complicated than the simple estimate provided in Eq. (31). The corotation torque can change the migration speed and also the direction of migration (Paardekooper & Mellema 2006; Baruteau & Masset 2008; Kley et al. 2009; Paardekooper et al. 2011). This can result in regions of the disc where the planet does not migrate at all, so-called zero-torque regions (Bitsch et al. 2013, 2014, 2015a; Baillié et al. 2015). However, the corotation torque is prone to saturation, which depends on the local viscosity of the protoplanetary discs (Paardekooper et al. 2011), where a lower viscosity allows an easier torque saturation, preventing outward migration. Even in these cases, however, the pebbles with τ_f > 0.001 drift faster than the planet migrates and will thus accumulate outside of the planetary orbit in the generated pressure bump.

Fig. 17 Total inward velocity of a dust particle (red) as a function of Stokes number τf and the corresponding gas velocity for a disc with α = 0.001 (magenta). Over-plotted are also the radial speeds of a 25 and a 50 ME planet in the same discs, assuming pure type I planetary migration. Even for the 50 ME planets, particles with τf > 0.001 drift faster than the planet and would thus accumulate at the pressure bump outside of the planetary orbit.

5.2 Massloading in the pressure bump

As soon as the planet reaches the pebble-isolation mass, pebbles drifting inwards from the outer disc are stopped in the pressure bump generated by the planet. As the flux of pebbles from the outer disc continues, pebbles accumulate in the pressure bump and the pebble-to-gas ratio increases. However, an increased pebble-to-gas ratio will trigger the streaming instability (Bai & Stone 2010; Carrera et al. 2015), transforming the pebbles into planetesimals. For the streaming instability to occur, a vertically integrated pebble-to-gas ratio of a few percent is needed.

Pebbles can also accumulate in vortices generated outside of the gap carved by the planet, where they would form planetesimals. Raettig et al. (2015) showed that the accumulated pebbles could destroy a vortex in a disc, but in our case, the vortex is fedby the presence of the planet itself, which was not taken into account in their work. Auffinger & Laibe (2018) studied the linear growth regime of the streaming instability in pressure bumps in discs, and they found that streaming instability can occur within the pressure bump. The accumulated pebbles inside the pressure bump therefore turn into planetesimals, which do not affect the gas velocities and thus do not disrupt the pressure bump. However, how the presence of a planet influences the streaming instability in a pressure bump is still subject to investigation.

Nevertheless, the current evidence occurrence of planetesimal formation inside the pressure bump before mass loading with pebbles influences the gas dynamics of the pressure bump itself. This makes the pressurebump outside of the planetary orbit an interesting candidate for subsequent planet formation.

5.3 Particle filtering by proto-Jupiters

The pebble accretion scenario does not only allow for fast accretion of planetary cores at large distances, it also gives potential solutions to (a) the great dichotomy between the terrestrial planets and the gas giants (Morbidelli et al. 2015), (b) the inward motion of the water ice line as the protoplanetary disc evolves in time and crosses the orbit of the Earth in less than 1 Myr (Morbidelli et al. 2016), and (c) explain the difference between the non-carbonaceous and carbonaceous meteorites through different isolated reservoirs (Kruijer et al. 2017). The solution to all these problems could be related to the growth of the Jupiter core and to the amount of pebbles that can drift past it after reaching pebble-isolation mass.

In these scenarios, the Jupiter core forms in the cold part of the protoplanetary disc (r > r_ice), where the pebbles are large. This makes the accretion very efficient because larger pebbles can be accreted more efficiently, allowing Jupiter to grow faster than the bodies in the terrestrial region (Morbidelli et al. 2015). Additionally, as soon as Jupiter reaches its pebble-isolation mass, the inward flux of large pebbles (τ_f > 10⁻¹) is stopped. The small pebbles drifting through are accreted very inefficiently (unless they grow again through coagulation), slowing down the growth of the bodies in the terrestrial region.

Morbidelli et al. (2016) reported that the inward flux of icy pebbles was stopped by a proto-Jupiter that had reached pebble-isolation mass, thus fossilizing the water ice line at ≈3 AU because the bulk of the material was stopped outside of the proto-Jupiter. However, a small fraction of water ice is needed to explain ordinary chondrites (and to a lesser extent, even enstatite chondrites) as they show evidence for some water alteration. The amount of water available had to be much lower than expected from solar proportion, however. This shows that these meteorites formed in a cold environment, but the availability of icy grains was somehow strongly reduced. The passage of small grains ( <10μm) through the barrier at the Jupiter pressure bump coupled with the inefficient accretion of such small grains explains these observations and is in agreement with the blocking efficiency of planets at pebble-isolation mass found in this study.

Kruijer et al. (2017) showed through meteoritic evidence that the reservoir between non-carbonaceous meteorites and carbonaceous meteorites was spatially separated in the protoplanetary disc around the young Sun at about ≈1 Myr. This separation can be achieved by a growing planet that stops the inward flux of particles, which corresponds to the core of Jupiter. Kruijer et al. (2017) gave the mass of the Jupiter core as ≈20 M_E . This is in agreement with thepebble-isolation mass we found here, but the exact mass at which a growing Jupiter generated a pressure bump outside of its orbit depends on the disc properties (viscosity, aspect ratio, and pressure gradient).

5.4 Ice giant formation

As soon as the planet reaches its pebble-isolation mass, the envelope of the planets is no longer heated by infalling pebbles, and gas accretion can start (Lambrechts et al. 2014). This initial contraction of the gaseous envelope depends on the cooling for the envelope and with that on the opacity inside it (Piso et al. 2015). At low temperatures (below 1000 K), the opacity is dominated by the dust grains, where a larger amount of dust grains increases the opacity and thus prolongs the contraction of the planetary envelope. However, in the pebble accretion scenario, the planet blocks the influx of new pebbles, and only very tiny dust grains can reach the planet.

As soon as the main flux of pebbles onto the planet is stopped, the core stops to grow, but very small grainsmight still enter the planetary atmosphere and thus keep the opacity high, prolonging gas envelope contraction (Lambrechts & Lega 2017). However, even in the case of no diffusion, particles with τ_f < 0.005 can reach the planet (Fig. 8) and keep the opacity high. Only for larger planets can the pebble flux sufficiently be reduced. Additionally, this depends on the viscosity of the protoplanetary disc because of the diffusion of particles through the pressure bump, where higher viscosities allow a more efficient diffusion and the planet has to reach higher masses to block pebbles with the same Stokes number compared to planets in low-viscosity environments. This initial growth stage might be very important for the growth of ice giants, preventing them from immediately entering into rapid gas accretion and thus explain why Uranus and Neptune have envelopes of 10–15% of their mass without entering into runaway gas accretion.

Lambrechts et al. (2014) found that Uranus and Neptune may never have reached pebble-isolation mass and that this prevented them for accreting gas. Our findings indicate that these planets might have reached the pebble-isolation mass without diffusion limit (to stop efficient growth of the core), but the influx of small particles prevented an efficient cooling of the atmosphere and thus runaway gas accretion. However, this could only have occurred if the viscosity of the protoplanetary disc was very low (α ≈ 10⁻⁴ ) because otherwise the pebble-isolation mass is higher than the mass of Uranus and Neptune for typically expected disc aspect ratios (H∕r > 0.04) in the outer disc at late disc evolution stages. We note that the outer disc structure is dominated by stellar irradiation, so that it is independent of the disc viscosity and H∕r becomes smaller through the reduced stellar irradiation as the system ages (Bitsch et al. 2015a).

6 Summary

We have conducted 3D hydrodynamical simulations to measure the pebble-isolation mass as a function of the disc structure and turbulence strength. In particular, we investigated the dependence on the disc aspect ratio H∕r, the disc viscosity (parametrised through α), the radial pressure gradient ∂lnP∕∂lnr, and the particle size described by the Stokes number τ_f. Our findings generally agree with the results presented in Lambrechts et al. (2014), who inferred the cubic dependence on the disc aspect ratio H∕r, but we refined the pebble-isolation mass to more disc parameters and also confirmed our results in fully radiative discs with heating and cooling. In Eq. (26) we provide the pebble-isolation mass as a function of our investigated disc parameters, which is useful for planet formation simulations involving pebble accretion (Bitsch et al. 2015b; Levison et al. 2015; Chambers 2016; Matsumura et al. 2017).

Our findings result in a pebble-isolation mass that is up to a factor of 2–3 higher than found in Lambrechts et al. (2014) in the high-viscosity case (α ~ 0.01) and a factor of ≈1.5 higher than in 2D simulations, see Appendix B. A higher viscosity and a steeper radial pressure gradient both result in a higher pebble-isolation mass. For very low viscosities, our simulations match the results of Lambrechts et al. (2014). This results in higher core masses in planet formation simulations compared to previous simulations (Sect. 4) in discs with α =0.0054. Discs with higher viscosity thus better match the heavy element content of the ice giants in our own solar system, because even as the disc evolves and H∕r decreases, the pebble isolation mass stays high enough so that the ice giants did not reach pebble isolation mass. Additionally, discs with higher viscosity can more easily match overall the heavy element content of giant exoplanets (Thorngren et al. 2016).

We also investigated the penetration of particles of various Stokes number τ_f through the pressure bump by radial advection with the gas and through turbulent diffusion. In the absence of turbulentdiffusion, a planet that has reached pebble-isolation mass can readily block pebbles with τ_f > 0.005, while a mass higher by a factor two is necessary to block pebbles with Stokes numbers as low as τ_f ~ 10⁻⁴. Including turbulent diffusion of particles due to viscosity changes this picture. Depending on viscosity and the particle size, the generated pressure bump needs to be stronger. To block particles of τ_f = 0.05, a typical size in drift-limited pebble growth models (Birnstiel et al. 2012), the planetary mass has to be increased by up to a factor of ≈2 compared to the pebble-isolation mass without turbulent diffusion (Eq. (10)) for high-viscosity discs. In low-viscosity discs, blocking of particles with τ_f = 0.05 requires a much smaller increase of the planetary mass than in the pebble-isolation mass without diffusion (Fig. 11).

Acknowledgements

B.B. was supported by the Knut and Alice Wallenberg Foundation (grant 2012.0150). A.J. thanks the Knut and Alice Wallenberg Foundation (grants 2012.0150, 2014.0017, 2014.0048), the Swedish Research Council (grant 2014-5775) and the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) for their financial support. The computations were done on the “Mesocentre SIGAMM” machine, hosted by the Observatoire de la Côte d’Azur. We thank P. Pinilla for discussions of the diffusion of dust particles. We thank the referee John Chambers for the comments that helped to improve this work.

Appendix A Azimuthal disc structure

We present here the 2D surface density structure of a 25 M_E planet embedded in a disc with α = 0.001 and H∕r = 0.05 (top panel of Fig. A.1) as wellas the η value in the 2D configuration. The data of these simulations correspond to the purple line in Fig. 1, indicating that the pebble-isolation mass has been reached when calculating the azimuthally averaged η profile. However, as can be seen in the bottom panel of Fig. A.1, a negative η value and with it a blocking of inward-drifting pebbles can be achieved at all azimuthal values when the pebble-isolation mass is reached. In principle, there is a small region in parameter space that might allow pebbles to “tunnel” through the pressure barrier, but the planet and with it the spiral waves orbit the central star at a frequency of Ω _P , whereas the pebbles orbit with a frequency Ω(r), which results in the trapping of the pebbles inside the pressure bump.

Fig. A.1 Surface density (top) and η value (bottom) for discs with α = 0.001, H∕r = 0.05 with an embedded 25 ME planet. The planet is located at r = 1 and ϕ = −1. A planet with this mass opens a partial gap in the disc density and has reached the pebble-isolation mass (Fig. 1). The η value reaches negative values outside of the planetary orbit over the whole azimuthal range, indicating that taking the azimuthally average η value to compute the pebble-isolation mass is correct.

Appendix B Comparison to 2D simulations

The pressure in isothermal 3D simulations is related to the gas volume density , while for isothermal 2D simulations, the pressure is related to the gas surface density with . However, the gas surface density and the gas volume density are related through (B.1)

where H_g is the gas scale height. For power-law discs in 3D with Σ_g ∝ r^s and with H∕r = const., this implies ρ_g ∝ r^s−1 and . In 2D discs, the surface density needs to have s > 1 to invert the pressure, while in 3D s > 2 is needed, making it harder to open a pressure bump in the protoplanetary disc.

In Fig. B.1 we show the η value as function of orbital distance for 2D discs with α = 0.001, H∕r = 0.05, and different planetary masses, the same as in Fig. 1. Clearly, a much lower planetary mass is needed to generate a pressure bump outside of the planetary orbit. In 2D simulations, planetary masses of about a factor of 1.5 less are needed to open a pressure bump in the protoplanetary disc than in 3D simulations.

Fig. B.1 η parameter as function of orbital distance from the planet for α = 0.001, H∕r = 0.05 and different planetary masses in 2D discs. The location of the planet is always at r = 1. A negative η parameter indicates the pressure bump. Here a mass of about ≈10 Earth masses is needed to generate the pressure bump, which is about a factor of 1.5–2 lower than in the 3D case.

References

All Figures

Fig. 1 Pressure gradient parameter η as a functionof orbital distance from the planet for α = 0.001, H∕r = 0.05, and differentplanetary masses. The location of the planet is fixed at r = 1. A negative η parameter indicates the formation of a pressure bump and with it super-Keplerian gas velocities that stop the inwards drift of pebbles. Here a mass of about 25 Earth masses is needed to generate the pressure bump outside of the planetary orbit. Inside of the planetary orbit, η is always positive.

In the text

	Fig. 3 Pebble-isolation mass as a function of ∂lnP∕∂lnr of unperturbed discs with different surface density gradients. Steeper surface density slopes s result in more negative ∂lnP∕∂lnr values. All simulations have been performed for planets in discs with H∕r = 0.05, f = 0 and α = 0.001.
In the text

	Fig. 4 Pressure gradient parameter η as a function of orbital distance in a disc set-up with heating and cooling for two different planetary masses. The planet is placed at r = 1. The pebble-isolation mass is reached at ≈28.3 ME in the simulations, in good agreement with Eq. (10).
In the text

	Fig. 5 Pebble-isolation mass as a functionof the pressure gradient ∂lnP∕∂lnr and α in a disc with a constant H∕r = 0.05. The two black lines mark 20 and 50 ME. Clearly, higher values of viscosity and ∂lnP∕∂lnr result in significantly higher pebble-isolation masses.
In the text

Fig. 6 Midplane gas velocities as a function of orbital distance in discs with H∕r = 0.05, α = 0.001, and Σ ∝ r−0.5. Top plot: radial velocity in units of the sound speed at r = 1.0. A negative velocity indicates an inward flow of the gas. The red line indicates the radial velocity of the unperturbed disc. Bottom plot: azimuthal velocity in units of the Keplerian velocity. If vθ is larger than 1, the gas orbits super-Keplerian, indicating the pressure bump in the disc.

In the text

Fig. 7 Trajectories of pebbles with τf = 1.0 that started at (1.5;0) moving in the gas velocity field generated by planets with MP = 10 ME and MP = 25 ME . We plot the pebble trajectories in the rotating frame, and the planet position is marked by the black dot. The 10 ME planet does not generate a pressure bump, so the pebble drifts through all the way towards the inner boundary of the computational domain. The 25 ME planet, on the other hand, generates a pressure bump outside of its orbit and thus blocks the flow of pebbles.

In the text

Fig. 8 Evolution of the orbital distance of pebbles as they drift through a disc with an embedded 25 ME planet (top) and an embedded 30 ME planet (bottom). Even though the planet has generated a pressure bump outside of its orbit (maximum at r = 1.123, see Fig. 1), small pebbles with τf < 0.005 can drift through the pressure bump for the 25 ME planet, while the Stokes number of the particles has to be an order of magnitude smaller to drift through the pressure bump generated by the 30 ME planet.

In the text

Fig. 9 Evolution of the orbital distance of pebbles as they drift through a disc with an embedded 25 ME planet, where the radial velocity is fixed to − 1 × 10−5cs,0 (top) and − 3 × 10−5cs,0 (bottom). As the drift speed of the particles depends on the gas velocities, particles of different sizes cannot be blocked or drift through compared to the nominal case shown in Fig. 8. In particular, a slower gas flow allows for a more efficient particle trapping in the pressure bump generated by the planet.

In the text

	Fig. 10 Velocity perturbation Δv as a function of planet mass at the location of the pressure bump. A negative Δ v value indicates that the gas speed is super-Keplerian. The blue circles mark the lowest Stokes number of particles that can be stopped at the pressure bump assuming a radial gas velocity of 1 cm s−1 and following Eq. (16).
In the text

Fig. 11 Critical Stokes number of particles that are affected by turbulent diffusion as a function of absolute Δ v∕cs in the pressure bump and α. Particles with Stokes numbers lower than the critical value easily diffuse and can drift through a pressure bump. The black lines denote particles with Stokes numbers of 10, 1, 0.1, 0.01, and 0.001 from left to right. The white dots correspond to different planetary masses, where we have taken the maximum value of the pressure bump generated by the planet for simulations with different α. The planetary masses are marked with the white numbers next to the dots. The horizontal dashed black line markes our nominal α = 0.001.

In the text

	Fig. 12 Δv∕cs as function of orbital distance in discs with α = 0.001 and H∕r = 0.05 and embedded planets of different masses. Note that Δv∕cs = η∕(H∕r). The Δ v∕cs values shown here correspond to the white dots at α = 0.001 (dashed line) in Fig. 11.
In the text

Fig. 13 Stokes number of the pebbles in the simulations presented in Figs. 15 and 16 at a given time t and orbital distance r in the disc. The black lines mark Stokes numbers of 0.05 to 0.4. Note that each planet growth trajectory stars at a given point in r-t and the planet the moves to higher time values, indicating that the Stokes number of the pebbles accreted by the planet increases. Additionally, the planet migrates in the disc to smaller orbital distances, which also increase the Stokes number of the accreted particles. The over-plotted blue lines correspond to the growth tracks shown in Fig. 14, where the solid line corresponds to solid (pebble) accretion, the dot marks the pebble-isolation mass, and the dashed line corresponds to the gas accretion phase. The planetary growth track moves upwards to increasing time.

In the text

Fig. 14 Growth tracks of planets starting at several different initial positions in a disc that is already 1.5 Myr old for pebble-isolation masses given by Lambrechts et al. (2014) in Eq. (1) and for the new-found pebble isolation mass Miso in Eq. (26). The initial growth phase is the same for both isolation masses, so the growth tracks diverge only when the pebble-isolation mass (Eq. (1)) is reached, marked by the dots. The solid lines indicate pebble accretion, while the dashed lines mark gas accretion.

In the text

Fig. 16 Final core masses of planets Mc as a function of formation distance r0 and formation time t0 in the disc. The white lines correspond to core masses of 10, 15, and 20 ME (top to bottom) and are marked by white numbers. The black and blue lines have the same meaning as in Fig. 15. The top plot corresponds to the pebble-isolation mass stated in Lambrechts et al. (2014), while the bottom plot corresponds to the pebble-isolation mass in this work (Eq. (10)). Clearly, Eq. (10) allows for higher core masses, which is crucial for the formation of ice giants in our own solar system.

In the text

Fig. 15 Final masses of planets (total mass MP = Mc + Menv) as a function of formation distance r0 and formation time t0 in the disc. Planets that are below the dark blue line have reached pebble-isolation mass and can accrete gas. All planets that are below the white line have Mc < Menv, indicating that they have undergone runaway gas accretion. The black lines indicate the final orbital distance rf of the planet, namely 0.1, 0.5, 1.0, 5.0, 10.0, and 20.0 AU. The top plot corresponds to the pebble-isolation mass found in Lambrechts et al. (2014), while the bottom plot corresponds to the pebble-isolation mass in this work (Eq. (10)). The higher pebble-isolation mass suppresses the formation of gas giants in the very outer parts of the disc at late times.

In the text

	Fig. 2 Pebble-isolation mass as a function of α and for different aspect ratios H∕r. The pebble-isolation mass is fitted through different α values with asimple fit that also scales with .
In the text

Fig. 17 Total inward velocity of a dust particle (red) as a function of Stokes number τf and the corresponding gas velocity for a disc with α = 0.001 (magenta). Over-plotted are also the radial speeds of a 25 and a 50 ME planet in the same discs, assuming pure type I planetary migration. Even for the 50 ME planets, particles with τf > 0.001 drift faster than the planet and would thus accumulate at the pressure bump outside of the planetary orbit.

In the text

Fig. A.1 Surface density (top) and η value (bottom) for discs with α = 0.001, H∕r = 0.05 with an embedded 25 ME planet. The planet is located at r = 1 and ϕ = −1. A planet with this mass opens a partial gap in the disc density and has reached the pebble-isolation mass (Fig. 1). The η value reaches negative values outside of the planetary orbit over the whole azimuthal range, indicating that taking the azimuthally average η value to compute the pebble-isolation mass is correct.

In the text

	Fig. B.1 η parameter as function of orbital distance from the planet for α = 0.001, H∕r = 0.05 and different planetary masses in 2D discs. The location of the planet is always at r = 1. A negative η parameter indicates the pressure bump. Here a mass of about ≈10 Earth masses is needed to generate the pressure bump, which is about a factor of 1.5–2 lower than in the 3D case.
In the text