© ESO 2018

1 Introduction

Pebbles – solids in the mm-cm size range – are important ingredients for planet formation in protoplanetary disks as they can efficiently speed up the formation of planetary cores (e.g. Ormel & Klahr 2010; Johansen & Lacerda 2010; Lambrechts & Johansen 2012). Solids couple more or less to the disk gas depending on their size, internal density, and on gas density and temperature. The accretion of well-coupled solids on a planetary core is inefficient as they are fully dragged by the gas across the orbit of the planet. The accretion of decoupled solids is also inefficient because of their high velocity relative to the planet. But marginally coupled solids can be efficiently accreted as the drag force around the planetary core notably reduces their velocity relative to the planet (Ormel & Klahr 2010). This is the case for solids in the pebble size range in typical regions of planet formation.

Another interesting aspect of pebbles is their potential observational signatures in the (sub)mm continuum emission of protoplanetary disks, in particular when they are trapped at pressure maxima. These pebbles can form bright and dark emission rings when there is a radial pressure bump (e.g. Andrews et al. 2016) or lopsided emission structures if there is an azimuthal pressure bump such as a vortex in the disk gas (e.g. Fuente et al. 2017). Observing the cold dust emission of disks can give constraints about the physical process behind the dust trapping in these pressure bumps such as the possible existence of a planet (Regály et al. 2012; Akiyama et al. 2016).

The pebble isolation mass (PIM) is defined as the minimum planetary mass that prevents pebbles from drifting to and being accreted by a planet in a protoplanetary disk. As the planet grows, the spiral wakes that it generates in the protoplanetary disk progressively become shocks closer to the planet. This process leads to the opening of a partial gap in the gas disk about the orbital radius of the planet (Rafikov 2002). A pressure maximum builds up at the outer edge of the gap of the planet, where pebbles should be trapped (e.g. Haghighipour & Boss 2003; Pinilla et al. 2012). This trapping may have two major implications: the growth of the planet by pebble accretion ceases and the trapped pebbles can produce ring-like structures in the (sub-)mm continuum emission. Said differently, the PIM can be seen as the maximum mass of a planetary core for in-situ pebble-mediated planet formation, and as the minimum planetary mass that may produce observable ring-like features in observations of the radio continuum.

An accurate determination of the PIM, particularly its functional dependence on disk parameters, is especially important for global models of planet formation and evolution. Lambrechts et al. (2014) investigated the dependency of the PIM on the scale height of the disk and location of the planet, but not on turbulent viscosity of the disk. However, the formation of even a partial gap in the disk gas around the orbit of the planet is sensitive to the turbulent viscosity of the disk, as shown for instance in the gap-opening criterion formulated in Crida et al. (2006). Furthermore, the presence of a pressure maximum is a necessary condition for trapping solids that are marginally coupled to the gas beyond the orbit of the planet, but it is not sufficient. Turbulent diffusion also acts on solids and may kick them out of the equilibrium location defined by the pressure maximum. This depends on how strong turbulent diffusion is compared to the radial drift due to gas drag. Because turbulent diffusion depends both on the turbulent viscosity of the disk and the size of the solids via their stopping time (Youdin & Lithwick 2007), the PIM should therefore also depend on solids size. We note that very small particles, which are well coupled to the gas, are not trapped at the pressure maximum (e.g. Zhu et al. 2012) unless perhaps the planet is so massive that it keeps the gas outside the orbit of the planet from flowing across the gap.

The aim of this study is to examine how the turbulent viscosity of the disk impacts the PIM with a combination of hydrodynamical simulations and semi-analytical calculations. The manuscript is organized as follows: in Sect. 2, we describe our physical model and numerical set-up, and we present the results of simulations. Our gas plus dust simulations show that the inclusion of dust turbulent diffusion can increase the PIM by nearly an order of magnitude at high gas turbulent viscosities. In Sect. 3, we derive a semi-analytical expression for the PIMthat is valid for a broad range of disk scale heights, viscosities, and pebble sizes, and we discuss its accuracy. This expression can easily be used in population synthesis models of planet formation and evolution, or for the interpretation of radio disk observations. We finally summarize our results in Sect. 4.

2 Hydrodynamical simulations

The PIM is a priori a function of the disk’s scale height H, its turbulentkinematic viscosity ν, and the particles stopping time t_s. These quantities can be expressed as dimensionless parameters that are the disk’s aspect ratio h, α turbulent viscosity, and Stokes number St via $$h=\frac{H}{r},\alpha =\frac{\nu }{c_{\text{s}}H},St=\frac{t_{\text{s}}}{t_{\text{eddy}}},$$(1)

where r denotes the distance from the central star and c_s is the sound speed. The Stokes number in the disk midplane is defined as the ratio between the particles stopping time t_s (the time particles need to adjust their motion to the gas because of drag forces), and the eddy turnover time t_eddy, that is the correlation timescale of turbulent fluctuations, which is typically of the order of the orbital timescale (see e.g. Youdin & Lithwick 2007). For particle sizes much smaller than the molecular mean-free path (Epstein regime), which is the case in our semi-analytic calculation and our 2D simulations, the stopping time is t_s = πρ_ss∕2ΣΩ with ρ_s the particles internal density, s their physical radius, Σ and Ω the surface density and the angular velocity of the gas at the location of the particles (Weidenschilling 1977). Adopting t_eddy = Ω⁻¹ results in St = πsρ_s∕2Σ.

To examine how disk turbulence impacts the PIM, we adopted a two-step strategy. In the first step, we carried out series of 2D gas hydrodynamical simulations of disk-planet interactions for a range of values for h and α. More specifically, for each pairof {h, α}, we found the minimum planet mass for which the gas azimuthal velocity (v_φ) can reach or slightly exceed the Keplerian velocity (v_K) beyond the orbital radius of the planet due to the local deposition of angular momentum by the outer wake of the planet. This is the necessary condition for pebbles to be trapped beyond the orbit of the planet. We provide a fitting formula for this PIM as a function of h and α. In this step, the turbulent viscosity modifies the gas velocities through the gas momentum equation and results in different PIM for different viscosities. That means that in this step we investigated the indirect effect of turbulent viscosity on particle trapping. In the second step, we carried out a few 2D gas plus dust hydrodynamical simulations of disk-planet interactions with a range of dust sizes to investigate how much the turbulent diffusion of dust alters the PIM as inferred from gas-only simulations. In this step, we examined the direct influence of viscous turbulent diffusion on the pebble trapping in addition to the indirect effect that is measured in the first step. These gas plus dust simulations are carried out for a specific value of h, but for a range of α values. After describing the physical model and numerical set-up of our simulations in Sect. 2.1, we present the results of the gas simulations in Sect. 2.2 and those of the gas plus dust simulations in Sect. 2.3.

2.1 Physical model and numerical set-up

We carried out our gas and gas plus dust hydrodynamical simulations with the code Dusty FARGO-ADSG. It is an extended version of the public code FARGO-ADSG (Masset 2000; Baruteau & Masset 2008a,b), which includes dust (Baruteau & Zhu 2016).

Gas – The gas continuity and momentum equations are solved in a 2D grid with polar coordinates {r, φ}. We use 666 grid cells uniformly spaced in radius between 0.5a and 2.5a (a denotes the semi-major axis of the planet), and 600 cells evenly spaced in azimuth between 0 and 2π. The initial surface density profile of the gas is ${\Sigma }_0(r)={10}^{-4}(M_{\star }/a^2)\times {(r/a)}^{-1/2}$ with M_⋆ the mass of the central star. This corresponds to a gas surface density ~900 g cm⁻² at 1 au for a solar-mass star. The Toomre-Q parameter never falls below ≈25 throughout our series of simulations, and for this reason gas self-gravity is discarded. Furthermore, for simplicity we adopt a locally isothermal equation of state with the gas pressure given by $P=\Sigma c_s^2$, where the sound speed (or temperature) are fixed in time. The disk aspect ratio, h = c_s∕v_K, is assumed to be uniform, which translates to a temperature profile T ∝ r⁻¹. The aspect ratio h is varied from 0.03 to 0.06 by steps of 5 × 10⁻³. The α turbulent viscosities taken in this work are 5 × 10⁻⁴, 10⁻³, 2 × 10⁻³, 5 × 10⁻³, and 10⁻². This range of h and α values are meant to reflect the range of temperatures and rates of turbulent transport expected in the inner regions (r ≲ 20 AU) of protoplanetary disks (e.g. Simon et al. 2015).

Dust – When dust is included in the simulations, it is modelled as Lagrangian test particles that feel the gravity of the star and planet, as well as the gas drag. We do not consider the dust back-reaction on the gas in this study. This effect might be important in long-term studies as we briefly discuss in Sect. 4. Turbulence is modelled by applying stochastic kicks to the particles position following the method in Charnoz et al. (2011), in which the spatial dependence of the turbulent diffusion of the dust is discarded, however. Specifically, the radial kicks have mean ⟨δr⟩ = D_d dt ∂_rΣ∕Σ and standard deviation ${\sigma }_r=\sqrt{2D_{\text{d}}\text{d}t}$, and the azimuthal kicks have mean ⟨δφ⟩ = D_d dt ∂_φΣ∕(r²Σ) and standard deviation σ_φ = σ_r∕r. For the dust’s turbulent diffusion D_d, we use $D_{\text{d}}=\nu (1+4S{t}^2)/{(1+S{t}^2)}^2$ from Youdin & Lithwick (2007). When the particles stopping time becomes shorter than the hydrodynamical timestep set by the gas CFL condition in the simulations, the so-called short friction time approximation is used to update particles velocities (Johansen & Klahr 2005). In this work, we use 50 000 dust particles with a size distribution n(s) ∝ s⁻¹ for sizes s between 1 mm and 300 mm. The particles internal density is set to 2 g cm⁻³. These particle sizes and internal density correspond to initial Stokes numbers St ∈ [3 × 10⁻⁴ − 0.2]. The Stokes number of the particles changes as they drift or diffuse in the disk inversely proportionally to the gas surface density at the particles location. The size range of the particles is chosen such that, first, there are enough particles with St > 0.05 that are trapped at the edge of the gap of the planet and, second, there are still particles with St < 0.05 that are not trapped at the gap edge and drift to the inner disk. The St ≈ 0.05 threshold is meant to be a slightly lower bound for the typical range of Stokes numbers for which pebble accretion is most efficient, which is St ∈ [0.1 − 1] (see e.g. Baruteau et al. 2016). Said differently, since this study deals with the PIM, we do not investigate the trapping of particles with a smaller Stokes number, as they should not contribute significantly to pebble accretion (e.g. Lambrechts & Johansen 2014). However, one can estimate the PIM for smaller Stokes numbers from the semi-analytical formula presented in Sect. 3.

Planet – The planet is held on a fixed circular orbit in all our simulations, and its gravitational potential is smoothed over a softening parameter ϵ = 0.4H(a). Its mass is gradually increased over the first five orbits in the simulations to let the gas adjust to the insertion of the planet in the computational grid. To avoid reflections of the planet wakes near to the radial edges of the grid, so-called wave-killing zones are used in which gas fields are damped towards their initial radial profile (de Val-Borro et al. 2006). All calculations are carried out in a frame co-rotating with the planet.

Units and running times – Results of simulations are expressed in the following units: the mass unit is the mass of the central star (M_⋆), the length unit is the planet’s (fixed) semi-major axis (a), and the time unit is ${\Omega }_{\text{K}}^{-1}(a)$. We denote by q the planet-to-star mass ratio and by q_⊕ the Earth-to-Sun mass ratio (q_⊕≈ 3 × 10⁻⁶). The simulations with only gas, which are presented in Sect. 2.2, are run over t_run = 5000, 4000, 3000, 2000 or 1000 planet orbits for our models with lowest to highest viscosity. These running times are chosen such that the disk profiles reach a steady state. Simulations with gas plus dust, which are presented in Sect. 2.3, are restarts of the gas-only simulations, where dust particles are inserted between r = 1.4a and 1.41a uniformly in radius and azimuth.

2.2 Results of gas-only simulations: minimum planet mass to form a local pressure bump beyond the planet’s orbit

Our gas hydrodynamical simulations were carried out to find the minimum planet mass for which a radial pressure maximum forms beyond the orbit of the planet as a function of disk aspect ratio h and α turbulent viscosity. Asalready mentioned at the beginning of Sect. 2, this was carried out by looking for the minimum planet-to-star mass ratio (q_g) for which the gas azimuthal velocity slightly exceeds the Keplerian velocity beyond the planet’s orbit, for each pair of {h, α}; the subscript g in q_g indicates that this is the PIM inferred from the gas-only simulations without considering the effect of turbulent diffusion of the dust. Equivalently, since $v_{\phi }=v_{\text{K}}\sqrt{1-\eta }$ with η = −h²∂logP∕∂logr, another way of looking for q_g is to find the planet-to-star mass ratio from which η cancels out beyond the orbit of the planet. Figure 1 shows the azimuthally-averaged radial profile of η for α = 5 × 10⁻⁴, obtained when the planet mass reaches the PIM for different values of the disk’s aspect ratio h. As can be seen in the figure, η > 0 everywhere in the disk except at the location of the radial pressure maximum, where η ≈ 0.

The values of q_g obtained for all of our models are shown in Fig. 2. The quantity ${(q_{\text{g}}/h^3)}^2$ is displayed as a function of α to highlight a linear relationship between both quantities. We find the best agreement for the following expression: $${\left(\frac{q_{\text{g}}}{h^3}\right)}^2\approx 82.33\alpha +\mathrm{0.03.}$$(2)

Although this relationship shows that the PIM has a much stronger dependency on h than α, the gas turbulent viscosity can notably affect the PIM since α can take a large range of values in regions of planet formation. For example, in a disk with h = 0.04 near to the location of the planet, the PIM increases from ~6 to ~20 Earth masses when α increases from 5 × 10⁻⁴ to 10⁻².

Equation (2) can be regarded as the condition for the opening of a partial gap in the gas around the orbit of the planet. This certainly differs from the gap-opening criterion formulated in Crida et al. (2006), which assumes planets carve a gap with a gas density drop of about 90%. The minimum planet-to-star mass ratio (q_gap) that satisfies the gap-opening criterion of Crida et al. (2006), $$\frac{3}{4}\frac{h}{{(q/3)}^{1/3}}+\frac{50\alpha h^2}{q}=1$$(3)

is given by Eq. (10) in Baruteau et al. (2014), which is written as $$q_{\text{gap}}=100\alpha h^2{\left[{\left(\sqrt{1+\frac{3h}{800\alpha }}+1\right)}^{1/3}-{\left(\sqrt{1+\frac{3h}{800\alpha }}-1\right)}^{1/3}\right]}^{-3}.$$(4)

We find thatq_g∕q_gap varies from about 3 to 12% for our range values of h and α. This point is further accentuated in Fig. 3, which shows that the gas surface density depletion in our models is between 10 and 20%. From this we can take away that a planet that builds up a 20% drop in the gas surface density about its orbit has reached the PIM.

We point out that Eq. (2) implies small PIMs for low-viscosity disks. According to Eq. (2), for α = 0, the PIM ranges from ~1.5 Earth masses for h = 0.03 to ~12.5 Earth masses for h = 0.06. This is consistent with the results of Duffell & MacFadyen (2013) that basically every non-migrating planet is able to carve a partial gap or dip about its orbit in an inviscid disk.

Fig. 1 Illustration of the determination of the PIM via gas-only hydrodynamical simulations, where we find the minimum planet-to-star mass ratio (qg) for the formation of a pressure maximum beyond the orbit of the planet. The quantity η = −h2 ∂logP∕∂logr (averaged over azimuth) is shown as a function of orbital radius r (normalized to the planet’s semi-major axis a). Each curve is for a pair {qg, h} at α = 5 × 10−4. The shaded black area shows where η ≈ 0 for the practical determination of qg.

Fig. 2 Calculation of the PIM via gas-only hydrodynamical simulations (qg, PIM normalized to the star mass). Results are shown for the range of disk aspect ratios h and α turbulent viscosities considered in this work. The blue dashed curve shows the best fit for ${(q_{\text{g}}/h^3)}^2$ vs. α, given by Eq. (2).

Fig. 3 Gas surface density profile in all the gas-only hydrodynamical simulations where the planet takes its PIM. Line transparency varies with the α turbulent viscosity, and colour varies with the disk aspect ratio (see legend).

Fig. 4 Results of gas plus dust hydrodynamical simulations without dust turbulent diffusion, which show that a planet with a planet-to-star mass ratio q = qg (the normalized PIM that is inferred from our gas-only simulations) does trap solid particles with St > 0.05 at the pressure maximum beyond its orbit.

2.3 Results of gas plus dust simulations: importance of turbulent diffusion on dust trapping at the pressure bump

In the previous section, we used gas-only hydrodynamical simulations to obtain the minimum planet-to-star mass ratio (q_g) for a planet to form a pressure maximum beyond its orbit, which is a necessary condition for trapping pebbles. To ensure that planets with q = q_g trap pebbles in the absence of dust turbulent diffusion, we performed several gas plus dust simulations without dust turbulent diffusion for h = 0.04 and all our α values. Figure 4 demonstrates the trapping of particles with St ≳ 0.05 at the pressure maximum beyond the orbit of the planet, as expected. We note that the trapping location increases with decreasing the alpha turbulent viscosity, which is due to planets carving wider gaps in less viscous disks (e.g. Crida et al. 2006). The particles with St ≪ 0.05 are well coupled to the gas and have not had enough time to drift significantly over the simulation time.

While the formation of a radial pressure bump is a necessary condition to trap pebbles beyond the orbit of the planet, it is not sufficient. We also needed to make sure that the turbulent diffusion of the dust does not diffuse particles out of the pressure bump. Because of our large parameter space, finding the PIM via gas plus dust hydrodynamical simulations for every single set of {h, α, St} would be computationally expensive. We thus chose to focus on one particular aspect ratio, h = 0.04, and find the minimum planet-to-star mass ratio that effectively traps dust particles with St > 0.05 at the pressure bump. Those particles have a size s > 10 mm for all our α turbulent viscosities. As already argued in Sect. 2.1, the St = 0.05 threshold can be seen as a lower bound for the typical range of particle Stokes numbers for which pebble accretion is efficient.

For each value of α, we start from q = q_g (the PIM as inferred from the gas-only simulations presented in Sect. 2.2) and increase the planet mass until dust trapping for the aforementioned sizes is observed for the whole running time of our gas + dust simulations. The corresponding planet-to-star mass ratio is denoted by q_d, where the dsubscript highlights that it is the mass ratio of PIM to the central star as inferred from gas plus dust simulations. Dust trapping can be quantified by counting the number of particles (N) with size s > 10 mm that remainbeyond the orbital radius of the planet over time. We denote by N₀ the number of such particles at the beginning of the gas + dust simulations (at time t = t_run). If the ratio N∕N₀ stays close to unity for the whole running time of the gas plus dust simulations, the planet has then reached its actual PIM.

The results of our gas plus dust simulations are summarized in Fig. 5. The first two rows of panels show the location of the dust particles overplotted on the relative perturbation of the gas surface density. The first row indicates q = q_g and the second row indicates q = q_d. The alpha turbulent viscosity decreases from left to right throughout the panels. The comparison between these first two rows of panels makes it clear that for q = q_g, most of the large particles have diffused out of the pressure bump and have mainly drifted inwards towards the star. This is further highlighted in the third row of panels, which represents the time evolution of N∕N₀. Clearly for the smallest viscosities, q_d and q_g are in close agreement, but the discrepancy significantly increases with increasing viscosity. For instance, for α =10⁻², the PIM obtained from the gas plus dust simulations is about 8.5 times larger than when obtained from the gas-only simulations. It is actually about half a Jupiter mass in the former case, and about 20 Earth masses in the latter case for a solar mass central star.

We note that the above values of q_d are for the trapping of St > 0.05 particles at the pressure bump. A different (in particular smaller) threshold value should result in a different value for q_d. Instead of finding the Stokes- or size-dependence of q_d via gas plus dust simulations, we now turn to the use of analytical calculations, which are presented in the next section.

Fig. 5 Results of the gas plus dust hydrodynamical simulations with dust turbulent diffusion for h = 0.04. The alpha turbulent viscosity decreases from left to right. First and second rows show the relative perturbation of the gas surface density (in black and white) with dots showing the location of dust particles (the colour of the dots indicates particles size). The first row shows results for q = qg, that is when the mass of the planet equals the PIM inferred from the gas-only simulations (q⊕ is the Earth to Sun mass ratio). The second row represents q = qd, that is when the mass of the planet equals the PIM obtained from the gas plus dust simulations. In each panel the white dashed circle indicates the orbital radius of the planet. The third row of panels represent the time evolution of N∕N0, which is the ratio of particles number with size >10 mm that remain outside of the orbital radius of the planet.

3 Semi-analytical formula for the pebble isolation mass

In the previous section we showed that dust turbulent diffusion could significantly alter the PIM that is inferred from the necessary condition to form a radial pressure bump beyond the orbit of the planet. As dust turbulent diffusion depends on the Stokes number (in addition to the gas turbulent viscosity), the PIM should explicitly depend on the Stokes number (in addition to the gas turbulent viscosity and aspect ratio). To examine how the PIM depends on the Stokes number, we turn to analytical and semi-analytical calculations. The idea is that particles get trapped at the pressure bump outside th orbit of the planet if the radial drift velocity (v_{d rift}) of the particles inside of the pressure bump, which is positive, becomes larger than the root mean square turbulent velocity fluctuations (v_{t urb}). When such a condition is satisfied, the particles that are kicked into the gap by the turbulence drift back towards the pressure maximum. We use Eq. (37) of Youdin & Lithwick (2007) to calculate the radial profile of v_{t urb}, and their Eqs. (57) and (59) for that of v_drift. These equations are written as $$\begin{array}{lll}{v}_{\text{t}urb}\hfill & =\hfill & \sqrt{\frac{1+4S{t}^2}{{(1+S{t}^2)}^2}\alpha }h{v}_{\text{K}},\hfill \\ v_{\text{d}rift}\hfill & =\hfill & \frac{St}{1+S{t}^2}\frac{\partial \log P}{\partial \log r}{h}^2{v}_{\text{K}}.\hfill \end{array}$$

Equating these two velocities gives $$\frac{\partial \log P}{\partial \log r}{|}_{v_{\text{t}urb}=v_{\text{d}rift}}=\frac{\sqrt{\alpha }}{h}\sqrt{\frac{1}{S{t}^2}+4}.$$(7)

This is the minimum (positive) pressure gradient required in the gap to keep the particles with Stokes number St from diffusing inwards through the gap. Equation (7) shows that for St ≫ 1, ∂ logP∕∂logr tends to $2\sqrt{\alpha }/h$, independently of St, while it tends to $\sqrt{\alpha }/(Sth)$ for St ≪ 1. We thus expect the PIM to change behaviour when varying St around St = 1, as shown below.

To find the planet mass that produces the pressure gradient satisfying Eq. (7) in the gap, we need to know how the planet modifies the gas pressure profile, or actually the gas density profile, in a locally isothermal disk model (where the planet does not alter the gas temperature). We use the density profile given by Eqs. (12)–(16) of Duffell (2015) which is, overall, in fairly good agreement with the results of disk-planet simulations for partial gap-opening planets; we return to this point later in this section. This density profile can be written as Σ (r) = Σ₀(r) C_Σ(r), where Σ ₀ denotes the initial (unperturbed) density profile, and C_Σ is a function of q, h, α, and r through an integral function that depends on s = −∂logΣ₀∕∂logr and f = ∂logh∕∂logr. Assuming a steady-state viscous disk model, where the mass accretion rate Ṁ = 3πνΣ₀ is uniform, we have s = 1∕2 + 2f (since α is taken uniform), which implies that $$\frac{\partial \log P}{\partial \log r}=-\frac{3}{2}+\frac{\partial \log C_{\Sigma }}{\partial \log r}.$$(8)

Anticipatingthat C_Σ(r) and the PIM have a weak dependence on the disk’s flaring index f, as is shown below, we take f = 0 and s = 1∕2, as in our hydrodynamical simulations. Because C_Σ(r) does not have an easily tractable form, we use an iterative method to find the right planet-to-star mass ratio q such that Eqs. (7) and (8) are satisfied, that is such that $$\frac{\partial \log C_{\Sigma }}{\partial \log r}=\frac{3}{2}+\frac{\sqrt{\alpha }}{h}\sqrt{\frac{1}{S{t}^2}+4}.$$(9)

For each pair {h, α}, we first find the planet-to-star mass ratio q_g,A such that a radial pressure maximum builds up beyond the orbit of the planet (∂ logC_Σ∕∂logr ≈ 3∕2 from Eq. (8)). Then, for each Stokes number St, we increase the planet-to-star mass ratio from q = q_g,A until Eq. (9) is satisfied in the gap. The planet-to-star mass ratio that we obtain, which we denote by q_d,A, corresponds to the PIM-to-star mass ratio for the effective trapping of particles with Stokes number ≥ St beyond the planet gap. The subscript A in q_g,A and q_d,A indicates that these quantities are calculated using our semi-analytical method.

The calculation of q_g,A with the analytic density profile of Duffell (2015) is well approximated by $${\left(\frac{q_{\text{g},\text{A}}}{h^3}\right)}^2\approx 37.7\alpha +\mathrm{0.01.}$$(10)

By comparison with Eq. (2), we see that for each pair {h, α}, q_g,A is smaller than its value in our hydrodynamical simulations by a factor ≈1.5. This comparison is illustrated in Fig. 6.

Figure 7 shows the results of our semi-analytical calculation of the ratio q_d,A∕q_g,A. For the range of Stokes numbers, aspect ratios, and alpha turbulent viscosities that we explored, we obtain good agreement with the following expression: $$\frac{q_{\text{d},\text{A}}}{q_{\text{g},\text{A}}}\approx 1+0.20{\left(\frac{\sqrt{\alpha }}{h}\sqrt{\frac{1}{S{t}^2}+4}\right)}^{0.7},$$(11)

which, from Eq. (7), can be recast more simply as $$\frac{q_{\text{d},\text{A}}}{q_{\text{g},\text{A}}}\approx 1+0.20{\left(\frac{\partial \log P}{\partial \log r}\right)}_{v_{\text{t}urb}=v_{\text{d}rift}}^{0.7}.$$(12)

As shown by the lower panel in Fig. 7, the agreement between Eq. (11) and our results of semi-analytical calculations is within ≈10%. In particular, we note that for a given pair {h, α}, q_d,A∕q_g,A decreases with increasing St for St ≪ 1, and progressively becomes independent of St for St ≳ 1, in agreement with the qualitative behaviour of ∂logP∕∂logr given by Eq. (7) and discussed above. We finally stress that it is not entirely surprising that q_d,A∕q_g,A features the pressure gradient given by Eq. (7), since it is the parameter that sets the effective trapping of particles with a given Stokes number in the presence of dust turbulent diffusion.

Figure 8 compares the PIM inferred from the gas plus dust hydrodynamical simulations presented in Sect. 2.3 with the PIM obtained with the above semi-analytical calculation, and expressed via Eq. (11); we denote by q_d its value normalized to the mass of the central star. Results are for h = 0.04 and St ≥ 0.05 particles. We see that q_d and q_d,A are in good agreement overall, although q_d tends to be larger than q_d,A at large viscosities. For example, for α = 10⁻², q_d is about 2.9 times larger than q_d,A.

The difference between q_d and q_d,A can be due to the discrepancy between the analytical gap density profile of Duffell (2015) and that in our hydrodynamical simulations. To investigate this, we show in Fig. 9 the radial profiles of the gas pressure (upper panel) and the logarithmic radial pressure gradient (lower panel) obtained from our simulations (solid curves) and from (Duffell (2015); dashed curves). Curves are shown for various alpha turbulent viscosities, but for the same planet mass and disk aspect ratio. For α ∈ [5 × 10⁻⁴ − 2 × 10⁻³], there is very good agreement in both the maximum value of ∂logP∕∂logr and its radial location, which explains why q_d∕q_d,A ≈ 1 for these values of α. For α = 5 × 10⁻³ and 10⁻², however, the agreement is not so good, the maximum in ∂logP∕∂logr being much larger and located closer to the planet in the analytical profile than in the simulations.Consequently, a smaller planet mass is needed to satisfy the condition v_drift ≃ v_turb than the PIMs obtained in the simulations, which explains why q_d∕q_d,A > 1 for these values of α. Based on this comparison, we thus suggest that Eq. (11) should be considered as a lower limit for the PIM in high-viscosity disks, and a good estimate for the PIM in low-viscosity disks.

The expressions of q_g,A and q_d,A in Eqs. (10) and (11) are for a steady-state viscous disk model with flaring index f = 0. As stated above, the background pressure profile scales as r^−3∕2 for a steady-state viscous disk model, whatever the specific choice of density and temperature profiles. However, the perturbed gap density profile (C_Σ in Eq. (8)) depends on the gradients of the background density and temperature profiles in a way that does not feature the background pressure gradient. The PIM should therefore also depend on the background temperature gradient (or, equivalently, on the surface density gradient, since both quantities are related in a steady-state viscous disk model). We checked with our semi-analytical calculation that this dependence remains weak. This is illustrated in Fig. 10, which shows that q_d,A slightly decreases with increasing flaring index. This decrease is by ≈7% for f = 0.5, which corresponds to a uniform background temperature. We also verified this weak dependence by running gas plus dust simulations for a radiative disk model with background gas surface density in r^−15∕14 and temperature in r^−3∕7 (e.g. Bitsch et al. 2014), for which we find basically the same PIM as with our fiducial disk profiles.

Fig. 6 PIM (normalized to the Earth mass) without effects of dust turbulent diffusion. A comparison between our results of gas-only hydrodynamical simulations (solid curves, Eq. (2)) with the prediction of the analytical gap model (dashed curves, Eq. (10) of Duffell (2015);.

Fig. 7 Semi-analytical calculation showing the dependence of the PIM on dust turbulent diffusion, using the method in Sect. 3. Upper panel: ratio between the normalized PIM with dust turbulent diffusion (qd,A) and that without (qg,A), as a functionof the particles Stokes number. Both qg,A and qd,A are calculated using the analytical gap density profile of Duffell (2015). Lower panel: relative difference between the normalized PIM obtained with our semi-analytical calculation (qd,a) and that obtained with our fitting formula (qd,A fit, Eq. (11)).

Fig. 8 Ratio of the normalized PIM for St ≥ 0.05 particles obtained in our gas plus dust hydrodynamical simulations (qd) and by the semi-analytical calculation given by Eq. (11) (qd,A). Results are for a disk aspect ratio h = 0.04 and for all our α values.

Fig. 9 Radial profiles of the gas pressure (upper panel) and of the logarithmic radial pressure gradient (lower panel) from our hydrodynamical simulations (solid curves) and from the analytic gap model (dashed curves) of Duffell (2015).

Fig. 10 Dependence of the PIM on the disk’s flaring index f, as obtained from our semi-analytical calculation for a disk aspect ratio h = 0.04, turbulent viscosity α = 10−3 and for particles Stokes number St = 0.1.

4 Summary

In this work we have examined how disk turbulence affects the PIM. Our findings can be summarized as follows:

• 1.

  First, we carried out gas hydrodynamical simulations of planet-disk interactions to find the minimum planet mass tobuild up a radial pressure bump outside the orbit of the planet, as a function of disk aspect ratio h and α-turbulent viscosity.The formation of such a pressure bump is a necessary but not sufficient condition for trapping pebbles beyond the orbit of the planet. By denoting q_g the ratio between the PIM and the mass of the central star, our results of simulations are best reproduced by the followingexpression: ${(q_{\text{g}}/h^3)}^2\approx 82.33\alpha +0.03$. Planets that reach the PIM deplete the surface density in the disk gas around their orbit by 10 to 20%. Our semi-analytical calculation based on the gap density profile of Duffell (2015) gives a PIM-to-star mass ratio q_g,A such that ${(q_{\text{g},\text{A}}/h^3)}^2\approx 37.3\alpha +0.01$. The difference between q_g and q_g.A resides in a slightly steeper gap pressure profile in Duffell (2015) than in our hydrodynamical simulations.

• 2.

  In addition to the presence of a radial pressure maximum, particle trapping requires that the drift velocityremains larger than the root mean square turbulent velocity fluctuations in the gap. This sets a minimum pressure gradient in the gap, which is given by Eq. (7). Our semi-analytical calculation shows that with dust turbulent diffusion, the effective trapping of particles with Stokes number ≥St gives riseto a correction factor by which q_g,A should be multiplied to obtain the actual PIM-to-star mass ratio (q_d,A). This correction factor is given by Eq. (11). Our final expression for the PIM-to-star mass ratio is written

$$q_{\text{d},\text{A}}\approx \underset{q_{\text{g},\text{A}}}{\underbrace{h^3\sqrt{37.3\alpha +0.01}}}\times \left[1+0.2{\left(\frac{\sqrt{\alpha }}{h}\sqrt{\frac{1}{S{t}^2}+4}\right)}^{0.7}\right].$$(13)

• 3.

  We also carried out gas plus dust hydrodynamical simulations with a range of particle sizes and we find that the PIM given by Eq. (13) is in good agreement with our results of simulations for low-viscosity disks. This PIM takes, however, smaller values than in the simulations for high-viscosity disks (by up to a factor ~3 for α = 10⁻²).

On the simulation side, this work is based on 2D hydrodynamical simulations of disk-planet interactions for a non-migrating planet. Such simulations require the use of a softening length in the planet’s gravitational potential, which can be seen as a crude way to mimic the effects of a finite vertical thickness. Morbidelli & Nesvorny (2012), who performed 2D simulations with a softening length slightly larger than ours (ϵ = 0.6H instead of 0.4H), found that for α = 6 × 10⁻³ and h = 0.056 the PIM is about 50 Earth masses, while our Eq. (2) gives ~ 42 Earth masses. Another comparison can be made with Lambrechts et al. (2014), whose 3D simulations resulted in a PIM of ~20 Earth masses for α = 6 × 10⁻³ and h = 0.05, while we find ~30 Earth masses. This comparison indicates that our choice of softening length gives overall good agreement with 3D simulations. Furthermore, if the migration of the planet is considered, a possible flux of pebbles coming from the disk inside the orbit of the planet may increase the PIM, depending on how the migration timescale and the pebbles radial drift timescale compare. This particular point requires a dedicated study.

A study of the PIM based on 3D gas-only hydrodynamical simulations along with 2D integrations of particle trajectories in the disk midplane has been recently undergone by Bitsch et al. (2018). A detailed comparison of our results with those of Bitsch et al. (2018) is found in Appendix A.

On the analytical side, our semi-analytical calculation of the PIM discards dust-gas interactions in the spirals and the direct gravitational interaction of the planet on the dust particles. Both of these effects could alter the predicted PIMs.

The stability of the dust ring that forms by piling up and growth of the particles should also be investigated in future works. If the lifetime of the ring is short compared to the age of the disk, the observed ring in millimetre observationswould need more massive planets to sustain a stronger pressure bump. On the other hand, the streaming instability might lead to a periodic growth of the planet by frequent destruction of the outer pressure bump. Therefore, studying the long-term evolution of dust trapping at the pressure maximum of planets with masses about the PIM and in low-viscosity disks (such that the low diffusion allows dust-to-gas ratios to reach unity) might probably need to take into account effects of dust drag onto the gas. Dust fragmentation and growth, which also depend on the disk turbulent diffusion, could also alter the PIM. These effects should be examined in future studies.

Acknowledgements

This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. We are grateful to the referee for her/his constructivereport that helped clarify this paper. We would also like to thank Paola Pinilla for her useful comments.

Appendix A Comparison with Bitsch et al. (2018)

Very recently, Bitsch et al. (2018) have studied the PIM dependency on the disk aspect ratio h, the α turbulent viscosity, and the initial pressure gradient ∂logP₀∕∂logr using gas-only 3D hydrodynamical simulations. Upon integrating the trajectory of particles with different Stokes numbers in the midplane of the disk, these authors have shown that their PIM could indeed trap particles with St ≥ 0.01 at the outer edge of the gap of the planet in the absence of dust turbulent diffusion. Using analytical estimations, they have obtained a correction factor that should be applied to their PIM expression to account for the effects of dust turbulent diffusion. Because the PIM expression obtained in Bitsch et al. (2018) can be seen as the 3D counterpart to that of the present study, we compare here both PIM formulae.

Fig. A.1 PIM (in Earth masses) without effects of dust turbulent diffusion. A comparison between the results of our 2D hydrodynamical simulations (solid curves) and the 3D simulations (dashed curves) of Bitsch et al. (2018); for Σ ∝ r−0.5 and uniformaspect ratios.

To ease the comparison, we express the formula of Bitsch et al. (2018) with our notations. We denote by q_g3d their PIM-to-star mass ratio without dust turbulent diffusion, and by q_d3d that obtained when accounting for dust turbulent diffusion. Eqs. (10) and (11) of Bitsch et al. (2018) can be recast as $$q_{\text{g3d}}=25q_{\oplus }{f}_{\text{fit}},$$(A.1)

where $$f_{\text{fit}}={\left(\frac{h}{0.05}\right)}^3\left[0.34{\left(\frac{{\log }_{10}{\alpha }_3}{{\log }_{10}\alpha }\right)}^4+0.66\right]\left[\frac{5.5+s-f}{6}\right],$$(A.2)

with α₃ = 10⁻³. In calculating ∂ logP∕∂logr in Eq. (11) of Bitsch et al. (2018), we use $P=\rho c_s^2$ with $\rho ={\Sigma }_0/(\sqrt{2\pi }H)$.

In Fig. A.1, we compare q_g as obtained in our 2D simulations with q_g3d from the 3Dsimulations of Bitsch et al. (2018). Both of these show overall good agreement and similar behaviours when varying h and α. Our PIM values are a factor 1.5 to 2 times smaller than in the 3D simulations of Bitsch et al. (2018). As shown in Bitsch et al. (2018), this can be understood by a more difficult gap formation in a 3D disk model than in 2D. Bitsch et al. (2018) have also carried out a few 2D simulations, which are directly comparable to ours. For a background surface density profile in r^−1∕2, h = 0.05, and α =10⁻³, their Fig. B1 shows a PIM of about 17M_⊕, while we obtain 14M_⊕ in our 2D simulations. While this is an overall good agreement, the 20% relative difference can be due to different grid resolutions (we use about twice as many grid cells in the radial direction), or perhaps to different simulation times (which are not specified in Bitsch et al. 2018). As we have experienced during preliminary tests, the gap profile can evolve on long timescales (e.g. a few thousand planet orbits for α =10⁻³), and running times that are too short lead to overestimating the PIM.

Fig. A.2 PIM (in Earth masses) with effects of dust turbulent diffusion. A comparison between the results of our semi-analytical calculation (solid curves) and those of (Bitsch et al. 2018; dashed curves) for Σ ∝ r−0.5, uniform aspect ratios and two values of the turbulent viscosity of the disk: α = 10−2 (upper panel) and α = 5 × 10−4 (lower panel). The red stars indicate the PIM obtained from our 2D gas plus dust simulations (for h = 0.04).

To account for the effects of dust turbulent diffusion on the PIM, Bitsch et al. (2018) have used a similar strategy as ours, in the sense that they also look for the planet mass in their simulations from which the root mean square turbulent velocity fluctuations of the particles (v_turb) become equal to their radial drift velocity (v_drift) within the gap. For the radial drift velocity, they have used $$v_{\text{drift}}=St\frac{\partial \log P}{\partial \log r}{h}^2{v}_{\text{K}},$$(A.3)

which coincides with Eq. (6) in the limit St ≪ 1. For the root mean square turbulent velocity fluctuations, instead of Eq. (5), Bitsch et al. (2018) have used $$v_{\text{turb}}~D_{\text{d}}/H_{\text{d}},$$(A.4)

with H_d a radial length scale for particle trapping in the pressure maximum, which is taken equal to the gas pressure scale height H, and with the dust’s turbulent diffusion D_d taken equal to the gas kinematic viscosity ν. We note that our expression for D_d tends to ν only in the limit St ≪ 1. Also, we stress that in the same limit St ≪ 1, Eq. (A.4) gives v_turb ≈ αc_s, while our Eq. (5) gives $v_{\text{turb}}\approx \sqrt{\alpha }{c}_{\text{s}}$. Eq. (5) thus leads to (much) larger turbulent velocity fluctuations than Eq. (A.4). Using Eqs. (A.3) and (A.4), the pressure gradient in the gap such that v_turb = v_drift is $$\frac{\partial \log P}{\partial \log r}{|}_{v_{\text{t}urb}=v_{\text{d}rift}}=\frac{\alpha }{h}\frac{1}{St},$$(A.5)

which is a factor α^−1∕2 smaller than our expression given by Eq. (7) in the limit St ≪ 1, due to the different expressions for v_turb in that limit.

Using Eq. (A.5), Bitsch et al. (2018) have found that the PIM with dust turbulent diffusion is written $$\frac{q_{\text{d3d}}}{q_{\text{g3d}}}\approx 1+4.2\frac{\alpha }{St},$$(A.6)

which, using Eq. (A.5), can be recast as $$\frac{q_{\text{d3d}}}{q_{\text{g3d}}}\approx 1+0.21\left(\frac{h}{0.05}\right)\frac{\partial \log P}{\partial \log r}{|}_{v_{\text{t}urb}=v_{\text{d}rift}}.$$(A.7)

Interestingly, Eq. (A.7) formally resembles our expression given by Eq. (12), except mainly for the scaling with the pressuregradient (which in our case comes to the 0.7 power). However, because we have different expressions for the pressure gradient ∂logP∕∂logr, our PIM values should differ a priori.

Figure A.2 compares q_d3d (see Eqs. (A.6), (A.1), and (A.2)) with our expression for q_d (see Eq. (13)) for α = 10⁻² (upper panel) and α = 5 × 10⁻⁴ (lower panel). For the regime of Stokes numbers that we are interested in, our PIM with dust turbulent diffusion is overall in good agreement with Bitsch et al. (2018) for α = 10⁻². However, for α = 5 × 10⁻⁴, we predicted smaller PIMs than Bitsch et al. (2018), for which the PIM displays no dependency on Stokes number for that viscosity. The differences are due to the combined effects of (i) different prescriptions for calculating the particle velocities (particularly the turbulent velocity fluctuations, as shown above), (ii) different dependencies of the PIM with the pressure gradient (compare Eqs. (A.7) and (12)), and also (iii) a larger q_g in Bitsch et al. (2018).

A last point of comparison with Bitsch et al. (2018) is the effect of varying the initial density and temperature gradients onthe PIM. While we have specialized to steady-state viscous disk models where both gradients are related, Bitsch et al. (2018) have varied both gradients independently. In Fig. 10, we show with our semi-analytical method that the PIM slightly decreases with increasing the background temperature gradient ∂ logT∕∂logr. Bitsch et al. (2018) have found with their simulations that the PIM slightly decreases with increasing the initial pressure gradient, which looks consistent with our findings, except that this dependence is found with increasing the initial density gradient for a fixed temperature profile (f = 0, see their Fig. 3), whereas no dependence of the PIM is found when varying the initial temperature gradient for a fixed surface density profile (s = 1∕2, see their Sect. 2.5). The reason for this is not totally clear, and may have to do with the steady-state assumption in our disk models. As stated before, starting with a non-steady-state disk model implies larger simulation times for a proper estimate of the PIM.

References

All Figures

Fig. 1 Illustration of the determination of the PIM via gas-only hydrodynamical simulations, where we find the minimum planet-to-star mass ratio (qg) for the formation of a pressure maximum beyond the orbit of the planet. The quantity η = −h2 ∂logP∕∂logr (averaged over azimuth) is shown as a function of orbital radius r (normalized to the planet’s semi-major axis a). Each curve is for a pair {qg, h} at α = 5 × 10−4. The shaded black area shows where η ≈ 0 for the practical determination of qg.

In the text

	Fig. 2 Calculation of the PIM via gas-only hydrodynamical simulations (qg, PIM normalized to the star mass). Results are shown for the range of disk aspect ratios h and α turbulent viscosities considered in this work. The blue dashed curve shows the best fit for ${(q_{\text{g}}/h^3)}^2$ vs. α, given by Eq. (2).
In the text

	Fig. 3 Gas surface density profile in all the gas-only hydrodynamical simulations where the planet takes its PIM. Line transparency varies with the α turbulent viscosity, and colour varies with the disk aspect ratio (see legend).
In the text

	Fig. 4 Results of gas plus dust hydrodynamical simulations without dust turbulent diffusion, which show that a planet with a planet-to-star mass ratio q = qg (the normalized PIM that is inferred from our gas-only simulations) does trap solid particles with St > 0.05 at the pressure maximum beyond its orbit.
In the text

Fig. 5 Results of the gas plus dust hydrodynamical simulations with dust turbulent diffusion for h = 0.04. The alpha turbulent viscosity decreases from left to right. First and second rows show the relative perturbation of the gas surface density (in black and white) with dots showing the location of dust particles (the colour of the dots indicates particles size). The first row shows results for q = qg, that is when the mass of the planet equals the PIM inferred from the gas-only simulations (q⊕ is the Earth to Sun mass ratio). The second row represents q = qd, that is when the mass of the planet equals the PIM obtained from the gas plus dust simulations. In each panel the white dashed circle indicates the orbital radius of the planet. The third row of panels represent the time evolution of N∕N0, which is the ratio of particles number with size >10 mm that remain outside of the orbital radius of the planet.

In the text

	Fig. 6 PIM (normalized to the Earth mass) without effects of dust turbulent diffusion. A comparison between our results of gas-only hydrodynamical simulations (solid curves, Eq. (2)) with the prediction of the analytical gap model (dashed curves, Eq. (10) of Duffell (2015);.
In the text

Fig. 7 Semi-analytical calculation showing the dependence of the PIM on dust turbulent diffusion, using the method in Sect. 3. Upper panel: ratio between the normalized PIM with dust turbulent diffusion (qd,A) and that without (qg,A), as a functionof the particles Stokes number. Both qg,A and qd,A are calculated using the analytical gap density profile of Duffell (2015). Lower panel: relative difference between the normalized PIM obtained with our semi-analytical calculation (qd,a) and that obtained with our fitting formula (qd,A fit, Eq. (11)).

In the text

	Fig. 8 Ratio of the normalized PIM for St ≥ 0.05 particles obtained in our gas plus dust hydrodynamical simulations (qd) and by the semi-analytical calculation given by Eq. (11) (qd,A). Results are for a disk aspect ratio h = 0.04 and for all our α values.
In the text

	Fig. 9 Radial profiles of the gas pressure (upper panel) and of the logarithmic radial pressure gradient (lower panel) from our hydrodynamical simulations (solid curves) and from the analytic gap model (dashed curves) of Duffell (2015).
In the text

	Fig. 10 Dependence of the PIM on the disk’s flaring index f, as obtained from our semi-analytical calculation for a disk aspect ratio h = 0.04, turbulent viscosity α = 10−3 and for particles Stokes number St = 0.1.
In the text

	Fig. A.1 PIM (in Earth masses) without effects of dust turbulent diffusion. A comparison between the results of our 2D hydrodynamical simulations (solid curves) and the 3D simulations (dashed curves) of Bitsch et al. (2018); for Σ ∝ r−0.5 and uniformaspect ratios.
In the text

Fig. A.2 PIM (in Earth masses) with effects of dust turbulent diffusion. A comparison between the results of our semi-analytical calculation (solid curves) and those of (Bitsch et al. 2018; dashed curves) for Σ ∝ r−0.5, uniform aspect ratios and two values of the turbulent viscosity of the disk: α = 10−2 (upper panel) and α = 5 × 10−4 (lower panel). The red stars indicate the PIM obtained from our 2D gas plus dust simulations (for h = 0.04).

In the text