# How much does turbulence change the pebble isolation mass for planet formation?

S. Ataiee, C. Baruteau, Y. Alibert, W. Benz

► To cite this version:

S. Ataiee, C. Baruteau, Y. Alibert, W. Benz. How much does turbulence change the pebble isolation mass for planet formation?. Astronomy & Astrophysics - A&A, 2018, 615, pp.A110. 10.1051/0004-6361/201732026. [hal-02369647](https://hal.science/hal-02369647v1)

HAL Id: hal-02369647

https://hal.science/hal-02369647v1

Submitted on 19 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# How much does turbulence change the pebble isolation mass for planet formation?

S. Ataiee<fnref n="1" />, C. Baruteau<fnref n="2" />, Y. Alibert<fnref n="1" />, and W. Benz<fnref n="1" />

Received 2 October 2017 / Accepted 24 March 2018

## ABSTRACT

*Context.* When a planet becomes massive enough, it gradually carves a partial gap around its orbit in the protoplanetary disk. A pressure maximum can be formed outside the gap where solids that are loosely coupled to the gas, typically in the pebble size range, can be trapped. The minimum planet mass for building such a trap, which is called the pebble isolation mass (PIM), is important for two reasons: it marks the end of planetary growth by pebble accretion, and the trapped dust forms a ring that may be observed with millimetre observations.

*Aims.* We study the effect of disk turbulence on the PIM and find its dependence on the gas turbulent viscosity, aspect ratio, and particles Stokes number.

*Methods.* By means of 2D gas hydrodynamical simulations, we found the minimum planet mass to form a radial pressure maximum beyond the orbit of the planet, which is the necessary condition to trap pebbles. We then carried out 2D gas plus dust hydrodynamical simulations to examine how dust turbulent diffusion impacts particles trapping at the pressure maximum. We finally provide a semi-analytical calculation of the PIM based on comparing the radial drift velocity of solids and the root mean square turbulent velocity fluctuations around the pressure maximum.

*Results.* From our results of gas simulations, we provide an expression for the PIM vs. disk aspect ratio and turbulent viscosity. Our gas plus dust simulations show that the effective PIM can be nearly an order of magnitude larger in high-viscosity disks because turbulence diffuse particles out of the pressure maximum. This is quantified by our semi-analytical calculation, which gives an explicit dependence of the PIM with Stokes number of particles.

*Conclusions.* Disk turbulence can significantly alter the PIM, depending on the level of turbulence in regions of planet formation.

**Key words.** protoplanetary disks – planet-disk interactions – submillimeter: planetary systems

## 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

<footnote n="1">¹ University of Bern, Physics Institute, Space Research and Planetary Sciences, Sidlerstrasse 5, 3012 Bern, Switzerland
e-mail: sareh.ataiee@space.unibe.ch</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS, UPS, Toulouse, France</footnote>
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 PIM that 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 turbulent kinematic viscosity $\nu$, and the particles stopping time $t_s$. These quantities can be expressed as dimensionless parameters that are the disk's aspect ratio $h$, $\alpha$ turbulent viscosity, and Stokes number $St$ via

$$h = \frac{H}{r}, \quad \alpha = \frac{\nu}{c_s H}, \quad St = \frac{t_s}{t_{\text{eddy}}}, \quad (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_{\text{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 = \pi \rho_s s / (2\Sigma\Omega)$ with $\rho_s$ the particles internal density, $s$ their physical radius, $\Sigma$ and $\Omega$ the surface density and the angular velocity of the gas at the location of the particles (Weidenschilling 1977). Adopting $t_{\text{eddy}} = \Omega^{-1}$ results in $St = \pi s \rho_s / (2\Sigma)$.

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 $\alpha$. More specifically, for each pair of $\{h, \alpha\}$, we found the minimum planet mass for which the gas azimuthal velocity ($v_\varphi$) 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 $\alpha$. 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 $\alpha$ 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, \varphi\}$. 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\pi$. The initial surface density profile of the gas is $\Sigma_0(r) = 10^{-4}(M_*/a^2) \times (r/a)^{-1/2}$ with $M_*$ the mass of the central star. This corresponds to a gas surface density $\sim 900 \text{ g cm}^{-2}$ at 1 au for a solar-mass star. The Toomre-Q parameter never falls below $\approx 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 \propto r^{-1}$. The aspect ratio $h$ is varied from 0.03 to 0.06 by steps of $5 \times 10^{-3}$. The $\alpha$ turbulent viscosities taken in this work are $5 \times 10^{-4}$, $10^{-3}$, $2 \times 10^{-3}$, $5 \times 10^{-3}$, and $10^{-2}$. This range of $h$ and $\alpha$ values are meant to reflect the range of temperatures and rates of turbulent transport expected in the inner regions ($r \lesssim 20 \text{ 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 $\langle \delta r \rangle = D_d \, dt \, \partial_r \Sigma / \Sigma$ and standard deviation $\sigma_r = \sqrt{2D_d \, dt}$, and the azimuthal kicks have mean $\langle \delta \varphi \rangle = D_d \, dt \, \partial_\varphi \Sigma / (r^2 \Sigma)$ and standard deviation $\sigma_\varphi = \sigma_r / r$. For the dust's turbulent diffusion $D_d$, we use $D_d = v(1 + 4St^2) / (1 + St^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) \propto s^{-1}$ for sizes $s$ between 1 mm and 300 mm. The particles internal density is set to $2 \text{ g cm}^{-3}$. These particle sizes and internal density correspond to initial Stokes numbers $St \in [3 \times 10^{-4} - 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 \approx 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 \in [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 $\epsilon = 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_\star$), the length unit is the planet's (fixed) semi-major axis ($a$), and the time unit is $\Omega_K^{-1}(a)$. We denote by $q$ the planet-to-star mass ratio and by $q_\oplus$ the Earth-to-Sun mass ratio ($q_\oplus \approx 3 \times 10^{-6}$). The simulations with only gas, which are presented in Sect. 2.2, are run over $t_{\text{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 $\alpha$ turbulent viscosity. As already 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, \alpha$}; 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_\varphi = v_K \sqrt{1 - \eta}$ with $\eta = -h^2 \partial \log P / \partial \log r$, another way of looking for $q_g$ is to find

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Illustration of the determination of the PIM via gas-only hydrodynamical simulations, where we find the minimum planet-to-star mass ratio ($q_g$) for the formation of a pressure maximum beyond the orbit of the planet. The quantity $\eta = -h^2 \partial \log P / \partial \log r$ (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 {$q_g, h$} at $\alpha = 5 \times 10^{-4}$. The shaded black area shows where $\eta \approx 0$ for the practical determination of $q_g$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Calculation of the PIM via gas-only hydrodynamical simulations ($q_g$, PIM normalized to the star mass). Results are shown for the range of disk aspect ratios $h$ and $\alpha$ turbulent viscosities considered in this work. The blue dashed curve shows the best fit for $(q_g/h^3)^2$ vs. $\alpha$, given by Eq. (2).</figcaption></figure>

the planet-to-star mass ratio from which $\eta$ cancels out beyond the orbit of the planet. Figure 1 shows the azimuthally-averaged radial profile of $\eta$ for $\alpha = 5 \times 10^{-4}$, 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, $\eta > 0$ everywhere in the disk except at the location of the radial pressure maximum, where $\eta \approx 0$.

The values of $q_g$ obtained for all of our models are shown in Fig. 2. The quantity $(q_g/h^3)^2$ is displayed as a function of $\alpha$ to highlight a linear relationship between both quantities. We find the best agreement for the following expression:

$$ \left(\frac{q_g}{h^3}\right)^2 \approx 82.33\alpha + 0.03. \qquad (2) $$
<figure><img src="image_4.png" /><figcaption>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).</figcaption></figure>

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 \qquad (3) $$

is given by Eq. (10) in Baruteau et al. (2014), which is written as

$$ q_{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} \qquad (4) $$

We find that $q_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.

<figure><img src="image_5.png" /><figcaption>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 = q_g$ (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.</figcaption></figure>

## 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 \ge 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 \ll 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 d subscript highlights that it is
<figure><img src="image_6.png" /><figcaption>**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 = q_g$, 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 = q_d$, 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/N_0$, which is the ratio of particles number with size >10 mm that remain outside of the orbital radius of the planet.</figcaption></figure>

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 remain beyond the orbital radius of the planet over time. We denote by $N_0$ the number of such particles at the beginning of the gas + dust simulations (at time $t = t_{run}$). If the ratio $N/N_0$ 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_0$. 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 $\alpha = 10^{-2}$, 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.

### 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 the orbit of the planet if the radial drift velocity ($v_{drift}$) of the particles inside of the pressure bump, which is positive, becomes larger than the root mean square turbulent velocity fluctuations ($v_{turb}$). 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_{\text{turb}}$, and their Eqs. (57) and (59) for that of $v_{\text{drift}}$. These equations are written as

$$v_{\text{turb}} = \sqrt{\frac{1 + 4St^2}{(1 + St^2)^2}} \frac{\alpha h v_K}{h}, \qquad (5)$$

$$v_{\text{drift}} = \frac{St}{1 + St^2} \frac{\partial \log P}{\partial \log r} h^2 v_K. \qquad (6)$$

Equating these two velocities gives

$$\left. \frac{\partial \log P}{\partial \log r} \right|_{v_{\text{turb}}=v_{\text{drift}}} = \frac{\sqrt{\alpha}}{h} \sqrt{\frac{1}{St^2} + 4}. \qquad (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 \gg 1$, $\partial \log P / \partial \log r$ tends to $2\sqrt{\alpha}/h$, independently of $St$, while it tends to $\sqrt{\alpha}/(St h)$ for $St \ll 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 $\Sigma(r) = \Sigma_0(r) C_\Sigma(r)$, where $\Sigma_0$ denotes the initial (unperturbed) density profile, and $C_\Sigma$ is a function of $q$, $h$, $\alpha$, and $r$ through an integral function that depends on $s = -\partial \log \Sigma_0 / \partial \log r$ and $f = \partial \log h / \partial \log r$. Assuming a steady-state viscous disk model, where the mass accretion rate $\dot{M} = 3\pi v \Sigma_0$ is uniform, we have $s = 1/2 + 2f$ (since $\alpha$ is taken uniform), which implies that

$$\frac{\partial \log P}{\partial \log r} = -\frac{3}{2} + \frac{\partial \log C_\Sigma}{\partial \log r}. \qquad (8)$$

Anticipating that $C_\Sigma(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_\Sigma(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}{St^2} + 4}. \qquad (9)$$

For each pair {$h, \alpha$}, 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 ($\partial \log C_\Sigma / \partial \log r \approx 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 $\ge 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_{g,A}}{h^3}\right)^2 \approx 37.7\alpha + 0.01. \qquad (10)$$

<figure><img src="image_7.png" /><figcaption>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)).</figcaption></figure>

By comparison with Eq. (2), we see that for each pair {$h, \alpha$}, $q_{g,A}$ is smaller than its value in our hydrodynamical simulations by a factor $\approx 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_{d,A}}{q_{g,A}} \approx 1 + 0.20 \left( \frac{\sqrt{\alpha}}{h} \sqrt{\frac{1}{St^2} + 4} \right)^{0.7}, \qquad (11)$$

which, from Eq. (7), can be recast more simply as

$$\frac{q_{d,A}}{q_{g,A}} \approx 1 + 0.20 \left( \frac{\partial \log P}{\partial \log r} \right)_{v_{\text{turb}}=v_{\text{drift}}}^{0.7}. \qquad (12)$$

As shown by the lower panel in Fig. 7, the agreement between Eq. (11) and our results of semi-analytical calculations is within $\approx 10\%$. In particular, we note that for a given pair {$h, \alpha$}, $q_{d,A}/q_{g,A}$ decreases with increasing $St$ for $St \ll 1$, and progressively becomes independent of $St$ for $St \ge 1$, in agreement with the qualitative behaviour of $\partial \log P / \partial \log r$ 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 \ge 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 $\alpha = 10^{-2}$, $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)
<figure><img src="image_8.png" /><figcaption>**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 ($q_{d,A}$) and that without ($q_{g,A}$), as a function of the particles Stokes number. Both $q_{g,A}$ and $q_{d,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 ($q_{d,a}$) and that obtained with our fitting formula ($q_{d,A fit}$, Eq. (11)).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Ratio of the normalized PIM for $St \ge 0.05$ particles obtained in our gas plus dust hydrodynamical simulations ($q_d$) and by the semi-analytical calculation given by Eq. (11) ($q_{d,A}$). Results are for a disk aspect ratio $h = 0.04$ and for all our $\alpha$ values.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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).</figcaption></figure>

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 $\alpha \in [5 \times 10^{-4} - 2 \times 10^{-3}]$, there is very good agreement in both the maximum value of $\partial \log P / \partial \log r$ and its radial location, which explains why $q_d/q_{d,A} \approx 1$ for these values of $\alpha$. For $\alpha = 5 \times 10^{-3}$ and $10^{-2}$, however, the agreement is not so good, the maximum in $\partial \log P / \partial \log r$ 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_{\text{drift}} \approx v_{\text{turb}}$ than the PIMs obtained in the simulations, which explains why $q_d/q_{d,A} > 1$ for these values of $\alpha$. 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_\Sigma$ 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 $\approx 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.

## 4. Summary

In this work we have examined how disk turbulence affects the PIM. Our findings can be summarized as follows:
<figure><img src="image_11.png" /><figcaption>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 $\alpha = 10^{-3}$ and for particles Stokes number $St = 0.1$.</figcaption></figure>

1. First, we carried out gas hydrodynamical simulations of planet-disk interactions to find the minimum planet mass to build up a radial pressure bump outside the orbit of the planet, as a function of disk aspect ratio $h$ and $\alpha$-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 following expression: $(q_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_{g,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 velocity remains 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 $\ge St$ gives rise to 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_{d,A} \approx h^3 \underbrace{\sqrt{37.3\alpha + 0.01}}_{q_{g,A}} \times \left[ 1 + 0.2 \left( \frac{\sqrt{\alpha}}{h} + \sqrt{\frac{1}{St^2} + 4} \right)^{0.7} \right]. \quad (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 $\sim 3$ for $\alpha = 10^{-2}$).

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 ($\epsilon = 0.6H$ instead of $0.4H$), found that for $\alpha = 6 \times 10^{-3}$ and $h = 0.056$ the PIM is about 50 Earth masses, while our Eq. (2) gives $\sim 42$ Earth masses. Another comparison can be made with Lambrechts et al. (2014), whose 3D simulations resulted in a PIM of $\sim 20$ Earth masses for $\alpha = 6 \times 10^{-3}$ and $h = 0.05$, while we find $\sim 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 observations would 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 constructive report that helped clarify this paper. We would also like to thank Paola Pinilla for her useful comments.

## References

Akiyama, E., Hasegawa, Y., Hayashi, M., & Iguchi, S. 2016, ApJ, 818, 158

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40

Baruteau, C., & Masset, F. 2008a, ApJ, 672, 1054

Baruteau, C., & Masset, F. 2008b, ApJ, 678, 483

Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927

Baruteau, C., Crida, A., Paardekooper, S.-J., et al. 2014, Protostars Planets VI, 667

Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, Space Sci. Rev., 205, 77

Bitsch, B., Morbidelli, A., Lega, E., & Crida, A. 2014, A&A, 564, A135

Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30

Charnoz, S., Fouchet, L., Aleon, J., & Moreira, M. 2011, ApJ, 737, 33

Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587

de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529
S. Ataiee et al.: How much does turbulence change the pebble isolation mass for planet formation?

Duffell, P. C. 2015, ApJ, 807, L11

Duffell, P. C., & MacFadyen, A. I. 2013, ApJ, 769, 41

Fuente, A., Baruteau, C., Neri, R., et al. 2017, ApJ, 846, L3

Haghighipour, N. & Boss, A. P. 2003, ApJ, 583, 996

Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353

Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107

Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35

Masset, F. 2000, A&AS, 141, 165

Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18

Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43

Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81

Rafikov, R. R. 2002, ApJ, 572, 566

Regály, Z., Juhász, A., Sándor, Z., & Dullemond, C. P. 2012, MNRAS, 419, 1701

Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588

Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6
<figure><img src="image_12.png" /><figcaption>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 $\Sigma \propto r^{-0.5}$ and uniform aspect ratios.</figcaption></figure>

## 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 $\alpha$ turbulent viscosity, and the initial pressure gradient $\partial \log P_0 / \partial \log r$ 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 \ge 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.

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_{g3d} = 25q_{\oplus}f_{\text{fit}}, \qquad (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], \qquad (A.2)$$

with $\alpha_3 = 10^{-3}$. In calculating $\partial \log P / \partial \log r$ 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 3D simulations of Bitsch et al. (2018). Both of these show overall good agreement and similar behaviours when varying $h$ and $\alpha$. 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.

<figure><img src="image_13.png" /><figcaption>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 $\Sigma \propto r^{-0.5}$, uniform aspect ratios and two values of the turbulent viscosity of the disk: $\alpha = 10^{-2}$ (upper panel) and $\alpha = 5 \times 10^{-4}$ (lower panel). The red stars indicate the PIM obtained from our 2D gas plus dust simulations (for $h = 0.04$).</figcaption></figure>

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 $\alpha = 10^{-3}$, their Fig. B1 shows a PIM of about $17M_\oplus$, while we obtain $14M_\oplus$ 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 $\alpha = 10^{-3}$), and running times that are too short lead to overestimating the PIM.

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_{\text{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}}, \qquad (A.3)$$

which coincides with Eq. (6) in the limit $St \ll 1$. For the root mean square turbulent velocity fluctuations, instead of Eq. (5), Bitsch et al. (2018) have used

$$v_{\text{turb}} \sim D_d/H_d, \qquad (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 $v$. We note that our expression for $D_d$ tends to $v$ only in the limit $St \ll 1$. Also, we stress that in the same limit $St \ll 1$, Eq. (A.4) gives $v_{\text{turb}} \approx \alpha c_s$, while our Eq. (5) gives $v_{\text{turb}} \approx \sqrt{\alpha c_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_{\text{turb}} = v_{\text{drift}}$ is

$$\left. \frac{\partial \log P}{\partial \log r} \right|_{v_{\text{turb}}=v_{\text{drift}}} = \frac{\alpha}{h} \frac{1}{S_t}, \qquad (A.5)$$

which is a factor $\alpha^{-1/2}$ smaller than our expression given by Eq. (7) in the limit $St \ll 1$, due to the different expressions for $v_{\text{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_{d3d}}{q_{g3d}} \approx 1 + 4.2 \frac{\alpha}{S_t}, \qquad (A.6)$$

which, using Eq. (A.5), can be recast as

$$\frac{q_{d3d}}{q_{g3d}} \approx 1 + 0.21 \left( \frac{h}{0.05} \right) \left. \frac{\partial \log P}{\partial \log r} \right|_{v_{\text{turb}}=v_{\text{drift}}}. \qquad (A.7)$$

Interestingly, Eq. (A.7) formally resembles our expression given by Eq. (12), except mainly for the scaling with the pressure gradient (which in our case comes to the 0.7 power). However, because we have different expressions for the pressure gradient $\partial \log P / \partial \log r$, 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 $\alpha = 10^{-2}$ (upper panel) and $\alpha = 5 \times 10^{-4}$ (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 $\alpha = 10^{-2}$. However, for $\alpha = 5 \times 10^{-4}$, 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 on the 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 $\partial \log T / \partial \log r$. 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.