© ESO 2020

1 Introduction

More than two decades have passed since the first discovery of an exoplanet (Mayor & Queloz 1995) and the number of confirmed exoplanets has now reached around 4000 and continues to increase¹. As the number of detections of exoplanets increases, the period–mass and period–radius distribution has shown nonuniformity in the occurrence rate of exoplanets. About half of Sun-like stars harbor close-in super-Earths with orbital periods of less than 85 days, radii of 1–4 R_⊕ (Earth radii), and masses of 2–20 M_⊕ (Earth masses; Fressin et al. 2013; Weiss & Marcy 2014). About 1% of Sun-like stars host hot Jupiters (≲0.1 au), while around 10% stars harbor cold Jupiters (~1–10 au) with a peak of occurrence rate of between 2 and 3 au (Johnson et al. 2010; Fernandes et al. 2019). The occurrence rate of the planets with an orbital distance of ~30–300 au are as low as 1% at most (Bowler 2016). Debate over the origin of exoplanetary systems is still ongoing.

Planets are formed in protoplanetary disks. One of the major planet formation theories currently being considered is the planetesimal accretion model. In this model, the building blocks of the planets are kilometer-sized planetesimals. As the planetary embryos accrete planetesimals and become massive enough, their gravitational focusing becomes more significant, and then the massive embryos grow faster than other smaller ones. This phase is called the runaway growth process (Kokubo & Ida 1996). After the runaway growth phase, the largest embryos (protoplanets) grow in an oligarchic manner, while most planetesimals remain small: this stage is called the oligarchic growth process (Kokubo & Ida 1998). The disk gas accretes onto a core formed in this way, and the total mass of the planet exceeds critical core mass (~ 10 M_⊕). Runaway gas accretion is then triggered and the planet evolves into a gas giant (e.g., Mizuno 1980; Pollack et al. 1996; Ikoma et al. 2000). In the planetesimal accretion model however, it is difficult to form gas giants within the lifetime of the disk at distances larger than ~5 au, unless the disk was massive (e.g., Kobayashi et al. 2011).

Pebble accretion is a new model of planet formation (e.g., Ormel & Klahr 2010; Ormel & Kobayashi 2012; Lambrechts & Johansen 2012, 2014; Lambrechts et al. 2014; Guillot et al. 2014; Ida et al. 2016), which may overcome several problems remaining in the theory based on planetesimal accretion (e.g., Kokubo & Ida 2000). Accreting millimeter-centimetre(mm-cm)-sized particles (pebbles) onto the proto-cores can form cores massive enough to trigger gas giant formation within the lifetime of gas disks (Lambrechts & Johansen 2012; Levison et al. 2015; Tanaka & Tsukamoto 2019; Bitsch et al. 2019). In addition to the formation of the cores of gas giants, the pebble accretion scenario has been applied to the size distribution of terrestrial planets in the Solar System (Morbidelli et al. 2015; Drążkowska et al. 2016), water delivery to the terrestrial planets (Morbidelli et al. 2016; Sato et al. 2016; Ida et al. 2019), the preference for prograde spin of the major planetary bodies in the Solar System (Johansen & Lacerda 2010; Visser et al. 2020), the formation of close-in super-Earths (Chatterjee & Tan 2014, 2015; Moriarty & Fischer 2015; Lambrechts et al. 2019; Izidoro et al. 2019; Bitsch et al. 2019), and the mass distribution of the planets around cool dwarf stars (Ormel et al. 2017; Schoonenberg et al. 2019). Pebble accretion in inviscid disks may explain the dichotomy between theinner super-Earths and outer gas giants (Fung & Lee 2018). However, the final configurations of the planetary systems strongly depends on the pebble flux (Bitsch et al. 2019). Even with the pebble accretion model, the origin of the nonuniformity of the distribution of exoplanets remains obscure as if it was enveloped in fog.

A key parameter to control the outcome of the planet formation models is the accretion probability of pebbles or smaller dust onto a growing body, which has been evaluated mostly in the Keplerian or sub-Keplerian shear flow around a protoplanet (Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014; Sellentin et al. 2013; Guillot et al. 2014; Ida et al. 2016; Visser & Ormel 2016). The loss of angular momentum due to the gas drag force leads to efficient accretion of pebble-sized particles compared to planetesimals. The particle scale height, which depends on the turbulence and particle size, has also been found to be important as it determines whether pebbles accrete two- or three-dimensionally.

Recent studies have reported the detailed 3D structure of the gas flow induced by a planet embedded in a protoplanetary disk (Ormel et al. 2015b; Fung et al. 2015; Lambrechts & Lega 2017; Cimerman et al. 2017; Kurokawa & Tanigawa 2018; Kuwahara et al. 2019; Béthune & Rafikov 2019). The anterior–posterior horseshoe flows extending in the orbital direction of the planet have a characteristic vertical structure like a column. A substantial amount of gas from the disk enters the gravitational sphere of the planet (Bondi or Hill sphere) at high latitudes (inflow), and exits through the midplane region of the disk (outflow).

The induced gas flow may alter the accretion probability of pebbles. Kuwahara et al. (2019) showed that the speed of midplane outflow increases with the protoplanet mass and eventually exceeds the terminal speed of incoming pebbles. This suggests that the outflow may inhibit accretion, whereas the polar inflow may act as the accretion window when the accretion is 3D. Popovas et al. (2018, 2019) incorporated dust particles into their hydrodynamical simulations and found that small particles (10 μm–1 cm) move away from the planet in the horseshoe flow and avoid accretion onto Earth- and Mars-sized planets. Ormel (2013) also reported a reduction in the accretion of small particles in 2D flow for an approximately Mars-sized core. From the analytical argument, Rosenthal et al. (2018) and Rosenthal & Murray-Clay (2018) proposed that the gas flow around the bound atmosphere may set the “flow isolation mass”.

Though these studies suggest the importance of considering the gas flow induced by a protoplanet, a comprehensive study investigating a range of core masses and pebble sizes is missing. Here, we perform hydrodynamical simulations of the gas flow around an embedded planet and compute the pebble trajectories and their accretion probabilities.

The structure of this paper is as follows. In Sect. 2 we describe the numerical method. In Sect. 3 we show the results obtained from a series of simulations and a comparison between the analytical estimations. In Sect. 4 we discuss the implications for planet formation. We summarize in Sect. 5.

2 Methods

2.1 Model overview

Gas flow around an embedded planet is perturbed by the gravity of the planet, and then a 3D flow structure is formed around the planet. In this study, we refer to the above-mentioned perturbed gas flow as the planet-induced (PI) gas flow, in contrast to the unperturbed shear flow. To investigate the influence of planet-induced gas flow on pebble accretion, we performed 3D hydrodynamical simulations (Sect. 2.2) and then calculated the trajectories of pebbles in the gas flow (Sect. 2.3). Finally, we compute the accretion probability of pebbles (Sect. 2.4). Through all of our simulations, the lengths, times, velocities, and densities are normalized by the disk scale height H, the reciprocal of the orbital frequency Ω⁻¹, the sound speed c_s, and the unperturbed gas density at the location of the planet ρ_disk, respectively.In this dimensionless unit system, the dimensionless mass of the planet (called “thermal mass”; Fung et al. 2015) is expressed by the ratio of the Bondi radius of the planet, R_Bondi, to the disk scale height, $$m=\frac{R_{\text{Bondi}}}{H}=\frac{G{M}_{\text{pl}}}{c_{\text{s}}^3\Omega },$$(1)

where G is the gravitational constant and M_pl is the mass of the planet. The Hill radius is given by $R_{\text{Hill}}={(m/3)}^{1/3}$ in this unit. When we assume a solar-mass host star and a disk temperature profile $T=270{(a/1\text{au})}^{-1/2}$ K, which corresponds to the minimum mass solar nebula (MMSN) model (Weidenschilling 1977b; Hayashi et al. 1985), M_pl is described by $$M_{\text{pl}}\simeq 12m{\left(\frac{a}{1\text{au}}\right)}^{3/4}{M}_{\oplus },$$(2)

where a is the orbital radius (Kurokawa & Tanigawa 2018).

2.2 Three-dimensional hydrodynamical simulations

In this study, we performed nonisothermal 3D hydrodynamical simulations of the gas of the protoplanetary disk around a planet. Our simulations were performed in a spherical polar coordinate co-rotating with a planet with Athena++ (White et al. 2016; Stone et al., in prep.). Most of our methods of hydrodynamical simulations are the same as described in detail in Kurokawa & Tanigawa (2018), but include several differences. Here we summarize the differences with Kurokawa & Tanigawa (2018) as below.

We ignored the headwind of the gas for all of our simulations. The external force in the Euler equation does not contain the global pressure force due to the sub-Keplerian motion of the gas. This assumption is justified because we focus on the shear regime of pebble accretion. There are two regimes of pebble accretion: the headwind regime and shear regime (Lambrechts & Johansen 2012; Ormel et al. 2017)². In the headwind regime, the approach speed of the pebble to the planet is dominated by the headwind of the gas, $$v_{\text{hw}}=\eta v_{\text{K}},$$(3)

where $$\eta =-\frac{1}{2}{\left(\frac{c_{\text{s}}}{v_{\text{K}}}\right)}^2\frac{\text{d}\ln P}{\text{d}\ln a},$$(4)

which is a dimensionless quantity characterizing the pressure gradient of the disk gas, v_K = aΩ, which is the Kepler velocity, and P is the pressure of the gas. On the other hand, the approach speed is dominated by the Keplerian shear in the latter case. The transition from the headwind to the shear regime occurs when the mass of the planet exceeds the transition mass given by $$M_{\text{t}}=\sqrt{\frac{1}{3}}\frac{v_{\text{hw}}^3}{G\Omega }\simeq 2.0\times {10}^{-4}{M}_{\oplus }{\left(\frac{a}{1\text{au}}\right)}^{3/2}{\left(\frac{v_{\text{hw}}}{30\text{m}{\text{s}}^{\text{-1}}}\right)}^3,$$(5)

(Lambrechts & Johansen 2012; Johansen & Lambrechts 2017). From Eq. (2), Eq. (5) can be rewritten as $$m\simeq 1.7\times {10}^{-5}{\left(\frac{a}{\text{1}\text{au}}\right)}^{3/4}{\left(\frac{v_{\text{hw}}}{30\text{m}{\text{s}}^{-1}}\right)}^3.$$(6)

Since the assumed planetary masses in this study are larger than the transition mass (Eq. (6)), we focus on the shear regime. We note that the effect of headwind on pebble drift is considered when we compute the accretion probability of pebbles (Sect. 2.4).

Kurokawa & Tanigawa (2018) fixed the size of the inner boundary for all of their simulations, but we varied it according to the mass of the planet. Assuming the density of the embedded planet, ρ_pl = 5 g cm⁻³ leads to the physical radius of the planet, R_pl, as given by $$\begin{array}{ll}{R}_{\text{pl}}\hfill & ={\left(\frac{3M_{\text{pl}}}{4\pi {\rho }_{\text{pl}}}\right)}^{1/3}\hfill \\ \hfill & =\left(\frac{9M_{\ast }}{4\pi {\rho }_{\text{pl}}}\right)\frac{R_{\text{Hill}}}{a}\left(\because R_{\text{Hill}}=a{\left(\frac{M_{\text{pl}}}{3M_{\ast }}\right)}^{1/3}\right)\hfill \\ \hfill & \simeq 3\times {10}^{-3}{m}^{1/3}{(\frac{{\rho }_{\text{pl}}}{5\text{g}{\text{cm}}^{-3}})}^{-1/3}{(\frac{M_{\ast }}{1M_{\odot }})}^{1/3}{(\frac{a}{1\text{au}})}^{-1},\hfill \end{array}$$(7)

where M_* and M_⊙ are the stellar mass and the solar mass. Our hydrodynamical simulations are computationally expensive, but the computational time is reduced when we increase the size of the inner boundary. Therefore, we regard the size of r_inn as being determined by Eq. (7) with a = 0.1 au, though we consider pebble accretion not only at 0.1 au, but also at various orbital radii. We confirmed that the size of the inner boundary does not affect our results.

Our simulations utilized the β cooling model, where the temperature T relaxes toward the background temperature T₀ with the dimensionless timescale β (e.g., Gammie 2001). To determine an adequate β value for a certain planetary mass, we followed the discussion in Kurokawa & Tanigawa (2018). Considering the relaxation time β for the temperature perturbation whose wavelength is equal to the Bondi radius of the planet, we fixed β = 1 at m = 0.1 (Malygin et al. 2017) and assumed that β scales linearly with the square of m, namely, $\beta ={\left(m/0.1\right)}^2$ (see Kurokawa & Tanigawa 2018, for the discussion).

We list our parameter sets in Table 1. The range of the dimensionless planetary masses, m = 0.03–0.3, corresponds to planets of between three Mars masses and one super-Earth in size, M_pl = 0.36–3.6 M_⊕, orbiting a solar-mass star at 1 au. Since the rapid increase of the gravity of the planet in the unperturbed disk affects the results of the simulations (Ormel et al. 2015a), the gravity of the planet is gradually inserted into the disk at the injection time, t_inj = 0.5. The gravity term is given by (Ormel et al. 2015a,b; Kurokawa & Tanigawa 2018; Kuwahara et al. 2019) $$F_{\text{p}}=-\nabla \left(\frac{m}{\sqrt{r^2+r_{\text{s}}^2}}\right)\left\{1-\exp \left[-\frac{1}{2}{\left(\frac{t}{t_{\text{inj}}}\right)}^2\right]\right\},$$(8)

where r is the distance from the center of the planet and r_s is the smoothing length. We assumed r_s = 0.1m (Kurokawa & Tanigawa 2018). The numerical resolution is [logr,  θ,  ϕ] = [128, 64, 128]. We adopted a logarithmic grid for the radial coordinates, which has a higher resolution in the vicinity of the planet.

2.3 Three-dimensional orbital calculation of pebbles

We calculated the trajectories of pebbles influenced by the planet-induced gas flow in the frame co-rotating with the planet (Fig. 1). We performed three types of orbital calculations: one is performed in the unperturbed shear flow (the gas has an initial condition profile), and the others are performed with the hydro-simulations data. We explain the detail of the orbital calculations in Sect. 2.3.2. The origin of the system is located at the position of the embedded planet. In order to investigate the effect of planet-induced gas flow on the motion of pebbles, we performed parameter studies for three planet masses, m = 0.03, 0.1, and 0.3, and various pebble sizes using the numerical results obtained from the above-mentioned hydrodynamical simulations. We used the final state of the hydro-simulations data (t = t_end), where the flow field seems to have reached a steady state.

Fig. 1 Schematic picture of orbital calculation of pebbles. A planet is located at the origin of the co-rotational frame. The dashed line represents the outer boundary of the hydrodynamical simulations. The starting point of orbital calculation is beyond the outer boundary of hydrodynamical simulations. Its y-component is fixed at 40RHill. The initial velocity of the pebble is the same with the Keplerian shear velocity. The gas velocity is assumed to be the speed of the Keplerian shear both inside and outside of rout in the shear flow case in the Stokes regime (Shear case), but is switched to the gas velocity obtained from the hydrodynamical simulations within rout in the planet-induced flow case in the Epstein regime (PI-Epstein case), and in the Stokes regime (PI-Stokes case).

2.3.1 Equation of motion

The dimensionless equation of motion of a pebble is described by $$\frac{\text{d}v}{\text{d}t}=\left(\begin{array}{c}2v_y+3x\\ -2v_x\\ 0\end{array}\right)-\frac{m}{r^3}\left(\begin{array}{c}x\\ y\\ z\end{array}\right)+F_{\text{drag}},$$(9)

where v = (v_x, v_y, v_z) is the velocity of the pebble and (x, y, z) are the coordinates of the pebble (Ormel & Klahr 2010; Visser & Ormel 2016). The first and second terms on the right-hand side of Eq. (9) are the Coriolis and tidal forces and the two-body interaction force with the embedded planet. The third term is the gas drag force acting on the pebble expressed by $$F_{\text{drag}}=-\frac{v-v_{\text{g}}}{\text{St}},$$(10)

where v_g is the gas velocity and St is the dimensionless stopping time of a pebble, referred to as the Stokes number, St = t_stopΩ. We assumed St = 10⁻³–1.

We omitted the dimensionless vertical component of the tidal force, − ze_z in Eq. (9). Through the test integrations with this term, we found that most pebbles settled in the midplane region of the disk before they enter the calculation domain of hydrodynamical simulations when St ≳ 0.1. In reality, turbulent stirring would, on average, balance the vertical component of the tide. The sources of this turbulence are discussed in Sect. 4.3.1. Though our hydrodynamical simulations do not resolve the turbulence, we assumed that the turbulent and tidal forces in the vertical direction canceled each other out. Considering the effect of random motion due to the turbulence (Xu et al. 2017; Picogna et al. 2018) is beyond the scope of this study.

The gas drag force is divided into two regimes, namely the Epstein and the Stokes regimes, depending on the relationship between the size of the particle and the mean free path of the gas. The stopping time of the particle in each regime is described by $t_{\text{stop}}=\{\begin{array}{c}\frac{{\rho }_{•}s}{{\rho }_{\text{g}}{c}_{\text{s}}},\left(\text{Epstein regime}:s\le \frac{9}{4}\lambda \right)\\ \frac{4{\rho }_{•}{s}^2}{9{\rho }_g{c}_s\lambda },\left(\text{Stokes regime}:s\ge \frac{9}{4}\lambda \right),\end{array}$

where ρ_• is the internal density of the pebble, s is the radius of the pebble, ρ_g is the density of the gas, and λ is the mean free path of the gas, $\lambda =\mu m_{\text{H}}/{\rho }_{\text{g}}{\sigma }_{\text{mol}}=1.44\text{cm}{(a/\text{1}\text{au})}^{11/4}$ with μ, m_H, and σ_mol being the mean molecular weight, μ = 2.34, the mass of the proton, and the molecular collision cross section, σ_mol = 2 × 10⁻¹⁵ cm² (Chapman & Cowling 1970; Weidenschilling 1977a; Nakagawa et al. 1986). The gas density at the midplane is given by ${\rho }_{\text{g}}={\Sigma }_{\text{g}}/\sqrt{2\pi }H$, where Σ_g is the gas surface density, ${\Sigma }_{\text{g}}=1700\text{g}{\text{cm}}^{-2}{(a/\text{1}\text{au})}^{-3/2}$. The pebble size in the MMSN can be expressed as $s=\{\begin{array}{c}3.4\times {10}^2\text{cm}\text{St}{\left(\frac{{\rho }_{•}}{2\text{g}{\text{cm}}^{-3}}\right)}^{-1}{\left(\frac{a}{1\text{au}}\right)}^{-3/2},\left(\text{Epstein regime}\right)\\ 33.2\text{cm}\sqrt{\text{St}}{\left(\frac{{\rho }_{•}}{2\text{g}{\text{cm}}^{-3}}\right)}^{-1/2}{\left(\frac{a}{1\text{au}}\right)}^{5/8},\left(\text{Stokes regime}\right)\end{array}$

(Lambrechts & Johansen 2012). From Eqs. (11)–(14), the transition from the Stokes to the Epstein regime occurs when $\text{St}=9.5\times {10}^{-3}\left(\frac{{\rho }_{•}}{2\text{g}{\text{cm}}^{-3}}\right){\left(\frac{a}{\text{1}\text{au}}\right)}^{17/4}.$(15)

Since the mean free path of the gas is proportional to the reciprocal of the gas density, $\lambda \propto {\rho }_{\text{g}}^{-1}$, the stopping time in the Stokes regime is independent of the gas density. In our model, the Stokes number, St, in Eq. (9) is independent of the local gas density, ρ_g, in the Stokes regime, whereas St is proportional to ${\rho }_{\text{g}}^{-1}$ in the Epstein regime. The gas density increases significantly in the vicinity of the planet, then the effective Stokes number decreases as the gas density increases. The threshold value of St for the transition varies as a function of the distance from the star (Eq. (15)). Because we discuss planet formation at various orbital radii, both the Epstein and Stokes regimes are considered for a range of St = 10⁻³–10⁰.

2.3.2 Numerical settings: 3D orbital calculation

We integrated Eq. (9) using a fifth-order Runge-Kutta-Fehlberg method (Fehlberg 1969; Shampine et al. 1976; Es-Hagh 2005). In this scheme each time step is controlled automatically depending on the relative acceleration of the pebble to the planet. The maximum relative error tolerance was set to 10⁻⁸, which ensures numerical convergence (Ormel & Klahr 2010; Visser & Ormel 2016). We fixed the y-coordinate of the starting point of pebbles at |y_s| = 40R_Hill (Ida & Nakazawa 1989). The x- and z-coordinates of the starting point of pebbles, x_s and z_s, are the parameters. We assumed that the initial velocity of the pebble is v_p,∞ = (0, −3∕2x, 0). We set the lower limits of z_s as z_s = 0 and x_s as x_s = 0.01b_x as values sufficiently small to be the common criteria for any planetary mass and Stokes number, which have on the order of ~ 10⁻⁴ [H], where b_x is the maximum impact parameter of accreted pebbles for the Shear case (Appendix A). The upper limits of x_s and z_s are set to be sufficiently large values, which have on the order of a few disk scale heights. We integrated Eq. (9) for pebbles with their initial spatial intervals of 0.01b_x in the x and z directions. Since the starting point is beyond the calculation domain of our hydrodynamical simulations, we assumed that the gas velocity was set to be Keplerian shear, v_g = (0, −3∕2x, 0), in the outer region of the domain, r ≥ r_out.

We performed three types of orbital calculation: shear flow case in the Stokes regime (Eq. (12); hereafter Shear case), wherewe adopted unperturbed Keplerian shear flow where the gas density is uniform; planet-induced flow case in the Epstein regime (Eq. (11); hereafter PI-Epstein case); and the same in the Stokes regime (hereafter PI-Stokes case). The gas density increases significantly in the vicinity of the planet in the PI-cases. Since the Stokes number in the Epstein regime depends on the gas density, the effective Stokes number decreases as the gas density increases. For the latter two cases, we switchedthe gas flow from Keplerian to planet-induced gas flow obtained by hydrodynamical simulations at r = r_out. We interpolatedthe gas velocity using the bilinear interpolation method (see Appendix B).

In the case of the m003 run, we found that the horseshoe flow formed unexpected vortices, which influences the pebble trajectories. The origin of these vortices is unknown, but it is likely to be a numerical artifact due to the spherical polar coordinates centered at the planet, in which the resolution becomes too low to resolve the horseshoe flow far from the planet when the assumed planet mass is small. Only in this case do we use the limited part of the calculation domain, r ≤ 0.3, to avoid the effects of the vortices. We discuss the width of horseshoe flow in Sect. 4.2.1.

We assumed that the pebble accreted onto the planet when the distance between the pebble and the planet became smaller than the following critical values: pebbles accreted when r ≤ 2r_inn for the Stokes regime, or r ≤ 0.3R_Bondi for the Epstein regime. In the Epstein regime, pebbles sometimes stagnate in the dense region near the planet. This is because the radial velocity of the pebble is so small but takes a nonzero value in the isolated envelope. To finish the calculation we assumed that pebblesaccreted onto the planet when they entered the bound atmosphere.

2.4 Calculation of accretion probability of pebbles

2.4.1 Width of accretion window and accretion cross section

We defined the width of the accretion window as $$w_{\text{acc}}(z)=x_{\text{max}}(z)-x_{\text{min}}(z),$$(16)

where x_max(z) and x_min(z) are the maximum and the minimum values of the x-component of the starting point of accreted pebbles at a certain height. In the unperturbed shear flow, the width of the accretion window, w_acc(z), is identical to the maximum impact parameter of accreted pebbles, b_x, when St < 1 as x_min(z) = 0 (Eq. (A.3)) Using this definition, we defined the accretion cross section as $$A_{\text{acc}}=2{\displaystyle {\int }_{z_{\text{min}}}^{z_{\text{max}}}{\displaystyle {\int }_{x_{\text{min}}(z)}^{x_{\text{max}}(z)}\text{d}}}x\text{d}z,$$(17)

where z_max and z_min are the maximum and the minimum value of the z-component of the starting point of accreted pebbles. The factor of two in Eq. (17) comes from the symmetricalstructure of the planet-induced gas flow (see Sect. 3.2 and Fig. 2). We reduced the spatial intervals stepwise near the edge of the accretion window, and determined x_max(z) and x_min(z) with sufficient accuracy. We confirmed that obtained w_acc has at least three significant digits.

Fig. 2 Streamlines of 3D planet-induced gas flow around the planet with m = 0.1 at t = 150. The red, green, and blue solid lines are the recycling streamlines, the horseshoe streamlines, and the Keplerian shear streamlines, respectively. The sphere is the Bondi sphere of the planet.

2.4.2 Accretion probability

Here we consider the accretion efficiency of pebbles in terms of the accretion probability of pebbles. We define the accretion probability of pebbles as $$P_{\text{acc}}=\frac{{\dot{M}}_{\text{p}}}{{\dot{M}}_{\text{disk}}},$$(18)

where Ṁ_p is the accretion rate of pebbles ontoa protoplanet and Ṁ_disk is the radial inward mass flux of pebbles in the gas disk described by $${\dot{M}}_{\text{disk}}=2\pi a{\Sigma }_{\text{p}}\left|v_{\text{drift}}\right|,$$(19)

where Σ_p is the surface density of pebbles and v_drift is the radial drift velocity of pebbles given by Weidenschilling (1977b) as $$v_{\text{drift}}=-\frac{2\text{St}}{1+{\text{St}}^2}{v}_{\text{hw}}.$$(20)

The density distribution of pebbles is described by $${\rho }_{\text{p}}(z)=\frac{{\Sigma }_{\text{p}}}{\sqrt{2\pi }{H}_{\text{p}}}\exp \left[-\frac{1}{2}{\left(\frac{z}{H_{\text{p}}}\right)}^2\right],$$(21)

where H_p is the scale height of pebbles (Dubrulle et al. 1995; Cuzzi et al. 1993; Youdin & Lithwick 2007): $$H_{\text{p}}={\left(1+\frac{\text{St}}{\alpha }\frac{1+2\text{St}}{1+\text{St}}\right)}^{-1/2},$$(22)

where α is the turbulent parameter in the disk introduced by Shakura & Sunyaev (1973). The sources of the turbulence are discussed in Sect. 4.3.1. Our calculation of accretion probability assumed that pebbles have a vertical distribution given by Eq. (21). This approach neglects the effect of random motion of individual particles, which is discussed in Sect. 4.2.3. The accretion rate of pebbles, Ṁ_p, divided into two formulas: $$\begin{array}{ll}{\dot{M}}_{\text{p},\text{2D}}\hfill & =2{\displaystyle {\int }_{x_{\text{min}}(0)}^{x_{\text{max}}(0)}{\Sigma }_{\text{p}}}{v}_{\text{p},\infty }\text{d}x,\hfill \end{array}$$(23)

in the 2D case, and $$\begin{array}{ll}{\dot{M}}_{\text{p},\text{3D}}\hfill & =4{\displaystyle {\int }_{z_{\text{min}}}^{z_{\text{max}}}{\displaystyle {\int }_{x_{\text{min}}(z)}^{x_{\text{max}}(z)}{\rho }_{\text{p},\infty }}}(z)v_{\text{p},\infty }\text{d}x\text{d}z,\hfill \end{array}$$(24)

in the 3D case. In order to account for the accretion from both x > 0 and from x < 0, we multiply Eq. (23) by two and Eq. (24) by four, respectively. The accretion probabilities in both 2D and 3D do not depend on the orbital radius, a (see Appendix C).

The existence of massive pebble reservoirs (~10² M_⊕) in the outer regions of the young disks has been suggested by (sub)millimeter dust observations (Ricci et al. 2010; Andrews et al. 2013; Ansdell et al. 2016). The dust mass is correlated with the accretion rate onto young stars. Given a typical T Tauri star with a disk accretion rate 10⁻⁸ M_⊙ Myr⁻¹, where M_⊙ is the solar mass, the dust mass is ~10² M_⊕ (Manara et al. 2016). The pebble flux can be analytically estimated as ~ 100–300 M_⊕ Myr⁻¹ when we assume the disk lifetime is ~1–3 Myr (Lambrechts & Johansen 2014). We fixed the inward pebble mass flux as Ṁ_disk = 10²M_⊕ Myr⁻¹, which is consistent with the typical value of the pebble flux used in a previous study (Lambrechts et al. 2019). Since the pebble mass flux is fixed but the radial drift velocity is not, the dust-to-gas ratio, Σ_p ∕Σ_g, varies as a function of the Stokes number and the orbital radius, $${\Sigma }_{\text{p}}/{\Sigma }_{\text{g}}=1.22\times {10}^{-5}\times (\frac{1+{\text{St}}^2}{\text{St}}){(\frac{a}{1\text{au}})}^{1/2},$$(25)

where the pebble surface density is derived from Eq. (19) with the aspect ratio of the disk, $H/a\simeq 0.033{(a/1\text{au})}^{1/4}$, and the typical value of the headwind, $v_{\text{hw}}/c_{\text{s}}\simeq 0.05{(a/1\text{au})}^{1/4}$, in the MMSN model.

3 Results

3.1 Results overview

The main subject of this study is to clarify the influence of the planet-induced gas flow on pebble accretion. In Sect. 3.2, we show the characteristic 3D structure of the planet-induced gas flow field obtained by 3D hydrodynamical simulations. In Sect. 3.3, we show the results of orbital calculations both in 2D and in 3D. Section 3.4 shows the dependence of the accretion probability of pebbles on the planetary mass and the Stokes number.

3.2 Three-dimensional planet-induced gas flow

The gas flow around an embedded planet is perturbed by the planet, which forms the characteristic 3D structure of the flow field (Ormelet al. 2015b; Fung et al. 2015; Cimerman et al. 2017; Lambrechts & Lega 2017; Kurokawa & Tanigawa 2018; Kuwahara et al. 2019; Béthune & Rafikov 2019). The fundamental features of the flow field are as follows: (1) gas from the disk enters the Bondi or Hill sphere at high latitudes (inflow) and exits through the midplane region of the disk (outflow; the red lines of Fig. 2). (2) The horseshoe flow exists in the anterior-posterior direction of the orbital path of the planet (the green lines of Fig. 2). The horseshoe streamlines have a columnar structure in the vertical direction. (3) The Keplerian shear flow extends inside and outside the orbit of the planet (the blue lines of Fig. 2).

The above-mentioned three features are commonly found both in the isothermal and in the nonisothermal simulations. However, anisolated envelope emerges around the planet in the nonisothermal case. The inner part of the envelope, whose size is ~ 0.5R_Bondi, is isolated from the recycling flow of the gas. The envelope has a slightly higher temperature than that of the disk gas due to compression, but we do not consider the thermal effect (e.g., evaporation) on pebble accretion. The profile of the gas flow in terms of the velocity, the density, the temperature, and the 3D structure does not change after several tens of orbits, which means that the flow field has reached a steady state. A detailed description and the differences between isothermal and nonisothermal 3D gas flows have been summarized in Kurokawa & Tanigawa (2018) and Kuwahara et al. (2019). In the following sections, we investigate the influence of the nonisothermal planet-induced gas flow on pebble accretion.

3.3 Orbital calculations

3.3.1 Pebble accretion in 2D

We first focus on the 2D limit of pebble accretion in the Stokes regime, that is all of the pebbles settled in the midplane of the disk and the Stokes number of pebbles does not depend on the gas density. Figure 3 shows the trajectories of pebbles at the midplane of the disk. The width of the accretion window becomes narrower in the PI-Stokes case than in the Shear case, in particular for the smaller Stokes numbers, St ≲ 0.1 (Figs. 3e–g). In the Shear case, the pebbles that successfully accrete onto the planet originated in the vicinity of the planetary orbit, x = 0. In the PI-Stokes case, however, most pebbles coming from this region move away from the planet along the horseshoe flows. This trend is consistent with the results of a previous study (Popovas et al. 2018). The pebbles coming from a window between the horseshoe and the shear regions can accrete onto the planet. We found voids in the upper left and lower right regions of the planet when the pebbles are small (Figs. 3e–g). The outflow speed of the gas at the midplane region of the disk is given by $\sqrt{3}/2m{c}_{\text{s}}$ for the range of planetary mass considered, which influences small pebbles when $\text{St}\lesssim \sqrt{m}$ as predicted by Kuwahara et al. (2019). The midplane outflow and resulting voids in pebble trajectories are not found in 2D simulations (Ormel 2013, their Fig. 12), but appear in the 3D ones.

For the larger Stokes number, St = 1, we obtained similar results in both the Shear and PI-Stokes cases (Figs. 3d and h). When the Stokes number becomes large, that is, when the stopping time of the pebble increases, pebbles are loosely coupled with the gas. Though the disk gas is perturbed by the planet, its influences on the larger pebbles become weak. Nevertheless, the motion of large pebbles is partially affected by the planet-induced gas flow on a large scale. We found that the shape of the accretion cross section of pebbles in the planet-induced gas flow differs from that in the unperturbed gas flow (see Sect. 3.3.2 and Fig. 7). In the PI-Epstein case, the shape of trajectories of pebbles does not differ significantly from that in the PI-Stokes case when we plot the trajectories with the same intervals as in Fig. 3: 0.05 [H]. However, the width of the accretion window, the accretion cross section, and the accretion probability in the PI-Epstein case do not match with those in the PI-Stokes case (see Sects. 3.4.1 and 3.4.2).

The dependence on the planetary mass in the PI-Stokes case is shown in Fig. 4. In our previous study, we found that the planet chiefly perturbs the surrounding disk gas at the typical scale of the smaller of the Bondi and Hill radii, and the outflow speed increases with the planetary mass (Kuwahara et al. 2019). In both Figs. 4a and b, the perturbed region can be scaled by the Bondi radius, but their sizes are different by an order of magnitude. When the planetary mass becomes larger, the perturbed region and the outflow speed become larger. Therefore, one can see the significant difference in Figs. 4a and b shown with the same scale. Whereas the pebbles with St = 10⁻² are hardly influenced by the perturbation in the gas flow induced by the planet when m = 0.03, the same-sized pebbles are highly influenced when m = 0.3.

Fig. 3 Trajectories of pebbles in Shear case (top) and PI-Stokes case (bottom) with different Stokes numbers at midplane around an embedded planet with m = 0.1. We set zs = 0 for all cases. The red and blue solid lines correspond to the trajectories of pebbles which accreted and did not accrete onto the planet, respectively. The dashed and dotted circles show the Hill and the Bondi radius of the planet, respectively. The black dots at the center of each panel denote the position of the planet. The interval of pebbles at their initial locations is 0.05 [H].

3.3.2 Pebble accretion in 3D

Next we focus on the 3D behavior of pebble accretion in the Stokes regime. Figure 5 shows the 3D trajectories of pebbles with St = 10⁻² around an embedded planet with m = 0.1 in the PI-Stokes case. The behavior of pebble accretion does not differ significantly from that in the midplane (Fig. 3f). Since the horseshoe streamlines extend in the vertical direction keeping their configurations (Fig. 2), pebbles coming from the vicinity of the planetary orbit move away from the planet and do not accrete onto it. The pebbles coming from the narrow region between the horseshoe and the shear regions sharply descend toward the planet with terminal velocity when they approach the planet from the polar region. This downward motion is due to the gravity of the planet and the advection due to the gas inflow (Fig. 2). Disk gas enters the Bondi sphere from the polar region and from the altitude where the planet gravity dominates (z ~ R_Bondi). The inflow has the velocity component of the − z direction, which causes rapid changes in the motion of the pebbles. The outflow of the gas also deflects the trajectories of pebbles and creates the voids in the second and fourth quadrants of the x–y plane because the outflows extend in the vertical direction (Kuwahara et al. 2019). The vertical scale of the outflow is ~ R_Bondi in the nonisothermal case.

Fig. 4 Trajectories of pebbles in the PI-Stokes case with St = 10−2 at the midplane around the embedded planets with m = 0.03 (left panel) and 0.3 (right panel). We set zs = 0. The red and blue solid lines correspond to the trajectories of pebbles which accreted and did not accrete onto the planet, respectively. The dashed and dotted circles show the Hill and the Bondi radius of the planet, respectively. The sizes of the Hill radius are 0.22 [H] for m = 0.03 and 0.46 [H] for m = 0.3. The black dots at the center of each panel denote the position of the planet. The interval of pebbles at their initial locations is 0.05 [H].

Fig. 5 Trajectories of pebbles with St = 10−2 around an embedded planet with m = 0.1 in the PI-Stokes case. The height of the initial position of the pebbles is zs = 0.7RHill. The red and blue solid lines correspond to the trajectories of pebbles which accreted and did not accrete onto the planet, respectively. The black dot denotes the position of the planet. The sphere around the planet represents the Hill radius of the planet. We only plot the trajectories within the region where r < 10RHill. The interval of pebbles at their initial locations is 0.05 [H].

3.4 Accretion probability of pebbles

3.4.1 Width of accretion window and accretion cross section

Figure 6 shows the changes of the width of the accretion window in the midplane region as a function of the Stokes number for different planetary masses m = 0.03, 0.1, and 0.3. In the Shear case, the numerically calculated width of the accretion window is almost consistent with the analytical estimation of the width of the accretion window (Eq. (A.4)). The width of the accretion window in the unperturbed shear flow follows a power law with an index of one-third.

In the PI-Stokes case, since the gas flow effects are weaker for the larger pebbles, the widths of accretion windows almost agree with those in the Shear case when St ≳ 0.01 (for m = 0.03) or St ≳ 0.1 (for m = 0.1, 0.3). On the other hand, the widths of accretion windows deviate from those in the Shear case when the Stokes number is small, St ≲ 10⁻²–10⁻¹. In the smallest Stokes number, St = 10⁻³, the width of the accretion window in the PI-Stokes case decreases by one or two orders of magnitude compared to that in the Shear case.

The reduction of the widths of accretion windows at the midplane becomes more prominent in the PI-Epstein case. This is because in the Epstein regime the Stokes number is proportional to the reciprocal of gas density. Since the gas density is higher around the planet due to its gravity, the effective Stokes number decreases as the pebble approaches the planet. In particular, pebbles with St = 10⁻³ do not accrete onto the planet with the spatial intervals of 0.01b_x in the x-direction. This means that the width of the accretion window is smaller than w_acc(0) ≲ 10⁻⁴ for St = 10⁻³. We truncate the width of the accretion window in PI-Epstein case at St = 3 × 10⁻³ in Fig. 6.

We plotted the accretion cross section of pebbles both in the Shear case and in the PI-Stokes case for a planet with m = 0.1 (Fig. 7). In both cases, the accretion cross section decreases with the Stokes numbers.

Firstly, we consider the results for the smaller Stokes numbers, St ≤ 0.1 (Figs. 7a–c, e–g). In the Shear case, since the pebbles approach the planet almost linearly (Fig. 3), the projected accretion cross sections are located near the planetary orbit. The value of x_min(z) changes only slightly, but the x_max(z) value decreases with height. Since the gravity of the planet acting on the pebbles becomes weaker at high altitudes, x_max(z) should decrease with height.

In the PI-Stokes case, the accretion cross section shifts to the right as a whole, and its shape becomes narrower compared to that in the unperturbed shear flow (Figs. 7e–g). This is due to the horseshoe flow, which lies across the orbit of the planet and has a columnar structure in the vertical direction. A notable result can be seen in the zoomed view in Fig. 7e. We found an accretion window above the midplane region. The width of the accretion window broadens at high altitudes, z ~ 0.2–0.3 [H]. The emergence of this window is related to the vertical scale of the outflow of the gas. In the midplane region, the pebbles with St = 10⁻³ are inhibited from accreting onto the planet by the planet-induced outflow. However, the pebbles approaching from an altitude higher than the vertical extent of the outflow (z ~ R_Bondi) can accrete onto the planet. Accretion windows at high latitudes above the outflow region for small pebbles were also found for other planet masses.

For the larger Stokes number, St = 1.0, the accretion cross section shifts to the right in both the Shear case and the PI-Stokes case (Figs. 7d and h). When the Stokes number becomes larger, the three-body effect becomes more important than the effect of the gas. Thus the behavior of the pebble is similar in both the Shear and the PI-Stokes cases in the region close to the planet near the midplane of the disk. However, at high altitudes, z ≳ 0.6 [H], we found the opposite trend between Figs. 7d and h. In the PI-Stokes case, the tip of the accretion cross section turns to the right in contrast to the result in the Shear case (Fig. 7d).

Figure 8 shows the differences between the integrated accretion cross section in the Shear, PI-Stokes, and PI-Epstein cases, and the dependence on the planetary mass. In the Shear case, the numerically calculated accretion cross section matches the analytical estimation (Eq. (A.8)). The accretion cross section in the unperturbed shear flow follows a power law with an index of two-thirds.

In the PI-Stokes case, the accretion cross sections agree with those in the Shear case when St ≳ 0.01 (for m = 0.03) or St ≳ 0.3 (for m = 0.1 and 0.3). On theother hand, the accretion cross sections deviate from those in the Shear case when the Stokes number is small, St ≲ 3 × 10⁻³–10⁻¹. For the pebbles with St = 10⁻³, the accretion cross section in the PI-Stokes case decreases by an order of magnitude compared to that in the Shear case.

Similarly to the 2D case (Fig. 6), the reduction of the accretion cross section becomes more prominent in the PI-Epstein case (Fig. 8). We note that, unlike Fig. 6, we can compute the accretion cross section for the pebbles with St = 10⁻³ because w_acc(z) takes a nonzero value above the outflow region.

Fig. 6 Width of accretion window in the midplane, wacc(0), as a function of the Stokes number in the PI-Stokes case (solid lines), PI-Epstein case (dashed-lines), and the Shear case (dotted lines). The dashed-dotted lines correspond to the analytical estimation for the Shear case expressed by Eq. (A.4). Colors indicate the mass of the planet: m = 0.03 (red), m = 0.1 (yellow), and m = 0.3 (blue).

Fig. 7 Accretion cross section with the different Stokes numbers for the planet with m = 0.1 in the Shear case (top) and the PI-Stokes case (bottom). We assumed α = 10−3 and the dust-to-gas ratio is equal to 10−2. Color represents the density of the pebbles expressed by Eq. (21) normalized by the gas density. The two panels in panels a and e show the enlarged outlines of accretion cross sections. We note that the color contour is saturated for the ρp ≲ 10−3.

Fig. 8 Accretion cross section as a function of the Stokes number as a function of the Stokes number in the PI-Stokes case (solid lines), PI-Epstein case (dashed-lines), and the Shear case (dotted lines). The dashed-dotted lines correspond to the analytical estimation for the Shear case expressed by Eq. (A.8). Colors indicate the mass of the planet: m = 0.03 (red), m = 0.1 (yellow), and m = 0.3 (blue).

Fig. 9 Accretion rate (left vertical axis) and probability (right vertical axis) as a function of the Stokes number in the PI-Stokes case (solid lines), PI-Epstein case (dashed-lines), and the Shear case (dotted lines). The dashed-dotted lines correspond to the analytical estimation for the Shear case expressed by Eqs. (A.12) and (A.13). Top: 2D case. Bottom: 3D case (α = 10−3). Colors indicate the mass of the planet: m = 0.03 (red), m = 0.1 (yellow), and m = 0.3 (blue).

3.4.2 Accretion probability

Figure 9 shows the accretion rate and the accretion probability as a function of the Stokes number for the different planetary masses. In the Shear case, the numerical results are consistent with the analytical estimations both in 2D and 3D (Eqs. (A.12) and (A.13)). The accretion rate (probability) increases as the planetary mass increases, but decreases as the Stokes number increases in 2D (Liu & Ormel 2018). In 3D however, the accretion probability has a peak at St ~ 0.03 (Ormel 2017; Ormel & Liu 2018).

In the PI-Stokes case, the situation changed. In the 2D case, the accretion probabilities either match or are slightly larger than that in the Shear case when St ≳ 3 × 10⁻³–10⁻², but deviate from it when St become smaller than the preceding value (Fig. 9a). Suppression of pebble accretion becomes significant when the Stokes number is small in the 2D limit. However, the coincidence in accretion probabilities between unperturbed and perturbed flow cases for St ≳ 3 × 10⁻³–10⁻² is counterintuitive, in particular for m = 0.1 and 0.3 because the width of the accretion window decreases significantly in the planet-induced gas flow when St ≲ 10⁻¹ for m = 0.1 and 0.3 (Fig. 6). This phenomenon can be understood as follows: in the Shear case, the accretion window lies near the planetary orbit, where the relative velocity between pebbles and the planet is low because the velocity is given by v_p,∞ = −3∕2xe_y. In the PI-Stokes case, the accretion cross section is shifted to the right (Fig. 7), where the relative velocity is large. In the midplane region, the offset of the reduction of the accretion cross section by the increase of relative velocity leads to the accretion probability in the perturbed flow case being almost equal to that in the shear flow for St ≳ 3 × 10⁻³–10⁻².

In the 3D case where we assumed α = 10⁻³ as a nominal value, we found the accretion probability in the PI-Stokes case either matches or is slightly larger than that in the Shear case for the planet with m = 0.1 and 0.3 (Fig. 9b). Similarly to the 2D case, the reduction of the accretion cross section and the increase of relative velocity cancel each other out. Since the width of the accretion window of the pebbles with St = 10⁻³ increases at high altitudes (Fig. 7e), the accretion probability in the planet-induced gas flow never falls below that in the unperturbed shear flow.

The above explanations shall not apply to m = 0.03, where the accretion probability for the smaller pebbles in the PI-Stokes case is smaller than that in the Shear case. Because the horseshoe width is proportional to the square of the planetary mass (Masset & Benítez-Llambay 2016), the shifting rate of the accretion cross section decreases as the planetary mass decreases. As a consequence, when the planetary mass is small, m = 0.03, the accretion cross section does not shift to the right enough to cancel its reduction.

In the PI-Epstein case, the reduction of the width of the accretion window and the accretion cross section become more significant than those in the PI-Stokes case for the smaller pebbles (Figs. 6 and 8). The increase of the relative velocities of pebbles does not fully offset the significant reduction of the width of the accretion window and the accretion cross section. Therefore, the accretion probability in the PI-Epstein case becomes smaller than that in the Shear- and PI-Stokes cases regardless of assumed St and m. Since the pebbles with St = 10⁻³ do not accrete onto the planet in the Epstein regime in the 2D case, the accretion probabilities are truncated at St = 3 × 10⁻³ in Fig. 9a. In the last paragraph of Sect. 3.3.1 we expected w_acc(0) ≲ 10⁻⁴ for St = 10⁻³. Based on Eq. (A.12), we have $$\begin{array}{ll}{\dot{M}}_{\text{p},\text{2D}}\hfill & \propto \left(x_{\text{max}}^2(0)-x_{\text{min}}^2(0)\right),\hfill \\ \hfill & =\left(x_{\text{max}}(0)+x_{\text{min}}(0)\right)w_{\text{acc}},\hfill \end{array}$$(26)

where x_max(0) + x_min(0) is expected to have an order of ~10⁻¹. If the pebbles can accrete onto the planet in the midplane with narrower spatial intervals than those investigated in this study, from these estimates we expect that the accretion probability for the pebbles with St = 10⁻³ is P_acc ≲ 4 × 10⁻⁴. This estimation proved to be valid when m = 0.03 and 0.1 (see the following paragraph and Fig. 11).

Figures 10 and 11 show the dependence on the turbulent parameter. Each panel shows the accretion probability for a planet with m = 0.03, 0.1, and 0.3. The accretion probability in the Shear case increases as α decreases and approaches that in the 2D case. When m = 0.03 (Fig. 10a), the accretion probability in the PI-Stokes case also increases as α decreases, and remains smaller than that in the Shear case. However, for m = 0.1 and 0.3 (Figs. 10b and c), the accretion probability for the smaller pebbles (St ≲ 10⁻³–10⁻²) becomes smaller than that in the Shear case, only when the turbulent parameter falls below α ≲10⁻⁵.

In the PI-Epstein case, the accretion probability for a pebble with St = 10⁻³ takes a much smaller value, P_acc ≲ 4 × 10⁻⁴, in the weak turbulence for a planet with m = 0.03 and 0.1 as mentioned above (Figs. 11a and b). The accretion probability for a planet with m = 0.3 converges at P_acc ~ 2 × 10⁻² in the range of the turbulence parameter used in this study (Fig. 11c). In the PI-Epstein case, accretion probabilities never exceed that in the Shear case even for strong turbulence, α =10⁻².

Figure 12 shows the accretion probability as a function of both the planetary mass and the Stokes number for the various turbulence strengths, α. Under strong turbulence, the accretion probability is always below P_acc ≲ 0.2 (Figs. 12a–c). The accretion probability has a peak at St ~ 0.1. As the turbulence strength decreases, the maximum value of the accretion probability increases and the peak shifts to smaller St (Figs. 12d–i). Suppression of pebble accretion due to the planet-induced gas flow becomes prominent for smaller Stokes numbers and smaller turbulence parameters. The peak of accretion probability lies in the upper left (higher m and smaller St) region of Fig. 12g and hasP_acc ≳ 0.5, whereas the accretion probability remains below P_acc ≲ 0.05 in the corresponding region of Fig. 12i.

Fig. 10 Dependence of accretion rate (left vertical axis) and accretion probability (right vertical axis) on the turbulent parameter, α. Panels a–c: results of the planet with m = 0.03, 0.1, and 0.3 in the PI-Stokes case (solid lines) and Shear case (dotted lines).

Fig. 11 Same as Fig. 10, but the dashed lines show the accretion probability in the PI-Epstein case.

4 Discussion

4.1 Comparison to previous studies

Recently, Popovas et al. (2018) investigated pebble accretion onto planets with masses of between those of Mars and Earth (m = 0.008, 0.04, and 0.08) during planet-induced gas flow in the Epstein regime. These latter authors performed 3D nonisothermal simulations in the shearing box. After sufficient relaxation of the fluid motion and density distribution of the gas, they injected millions of tracer particles with different sizes (St = 3 × 10⁻⁵–3). The initial stratification of particles was the same as that of the gas. The authors measured the accretion rates of a subpopulation of particles that initially resided within the Hill sphere and reported that the accretion rate increases as the planetary mass increases, and that it scales linearly with pebble size regardless of the assumed planetary mass when St = 3 × 10⁻⁵–3 × 10⁻¹. In their subsequent study, they concluded that the convective motion within the envelope caused by the accretion heating does not significantly affect the accretion rate (Popovas et al. 2019). Similarly to their results, here we found that the accretion rates increase with planetary mass (Fig. 11). In our study however, the slope of the accretion rate changes significantly for the assumed planetary mass and turbulence parameter. The accretion rate does not scale linearly over the same pebble size range as in Popovas et al. (2018). The maximum accretion rate in Popovas et al. (2018) is Ṁ ~ 10² M_⊕ Myr⁻¹ when m = 0.08 and St ~ 1–3, but the achieved accretion rate in our study becomes smaller by one or two orders of magnitude even for a larger planet, m = 0.3 (Figs. 11b and c). The higher accretion rate obtained in Popovas et al. (2018) may be the result of the absence of turbulent stirring. In contrast to our study, Popovas et al. (2018) considered the vertical component of the tidal force, − zΩe_z, so that the particles settle to the midplane with the settling timescale, ${(\text{St}\Omega )}^{-1}$, which may lead to a high accretion rate for the larger pebbles. We discuss the effect of the turbulence on the accretion probability in Sect. 4.2.3.

In our previous study, we found that the speed of outflow of gas at the midplane region can be expressed by $u_{\text{out}}=\sqrt{3/2}m{c}_{\text{s}}$ (Kuwahara et al. 2019). Comparing the outflow speed to the terminal velocity of pebbles, we found that planetary mass with the potential to affect the accretion can be written as $m\gtrsim \sqrt{\text{St}}$. As the massof the planet increases, the outflow becomes fast and starts to prevent solid materials from accreting onto the core. As shown in Figs. 6–8, the reduction of the accretion cross section becomes significant for the pebbles with St ≲ m², which is consistent with the prediction in Kuwahara et al. (2019). However, the outflow does not completely prevent pebble accretion. Most pebbles approach the planet from the first and third quadrants of the x–y plane. In contrast, the outflow occurs dominantly in the second and fourth quadrants of the x–y plane.

Fig. 12 Accretion probability as a function of the planetary mass and the Stokes number for the turbulence parameters α = 10−6, 10−4, and 10−2 (bottom to top) in the Shear case (left column), the PI-Stokes case (middle column), and the PI-Epstein case (right column). The contours represent the accretion probabilities.

4.2 Topologies of planet-induced gas flow

4.2.1 Width of horseshoe flow

The shift of the accretion cross section is controlled by the width of the horseshoe region given by Masset & Benítez-Llambay (2016), $$w_{\text{HS}}=1.05\sqrt{m}.$$(27)

Equation (27) agrees with the horseshoe width of our hydrodynamical simulations for m01 and m03 runs. In m003 run, Eq. (27) overestimated the horseshoe width. This overestimation of the horseshoe width for a lower-mass planet (m ≤ 0.03) could be seen in the previous studies (Ormel et al. 2015b; Kurokawa & Tanigawa 2018; Kuwahara et al. 2019), but if w_HS continues to decrease as the planetary mass decreases, the accretion cross section only shifts very slightly. Thus, whether the accretion probability changes or not depends on whether the outflow barrier works for pebbles in the vicinity of the planet. However, the outflow speed decreases with the planetary mass, and it becomes difficult for the outflow to suppress pebble accretion when the planetary mass is small (Kuwahara et al. 2019). Assuming that the lower limit of the Stokes number is St = 10⁻³, according to our previous study, we would expect that only a planet larger than m ≥ 0.03 would change the accretion probability (Kuwahara et al. 2019). Pebble accretion in the planet-induced gas flow for the lower mass planets may not be much different from that in the unperturbed shear flow, as long as we focus on the shear regime of pebble accretion. When the planetary mass is much smaller than what is considered in this study, the effect of the headwind becomes more important (see the following section).

4.2.2 Effect of headwind

Pebble accretion probability depends on the topologies of the planet-induced gas flow. In this study, we did not considerthe sub-Keplerian motion of the gas as we focused on the shear regime. Given the nonzero η (Eq. (4)), the flow topology changes (Ormel 2013; Ormel et al. 2015a,b; Kurokawa & Tanigawa 2018). The flow topology is nearly axisymmetric with respect to x = 0 when we assume η = 0, but the symmetry breaks downwhen η≠0. The horseshoe streamlines shift to within the planetary orbit. Such asymmetric flow structure may affect the accretion probability of pebbles. The shift of the horseshoe may enhance the accretion of pebbles coming from the front of a planet’s orbit, whereas the relative velocity of accreted pebbles (Sect. 3.4) may be reduced. Therefore, the net effect of the headwind is not obvious. Considering the headwind regime is beyond the scope of this study, but it would be important for understanding pebble accretion onto a core smaller than the transition mass (Eq. (5)).

4.2.3 Effect of turbulence

Strong turbulent viscosity also changes the flow topology. Fung et al. (2015) performed 3D viscous hydrodynamical simulations, finding that the outward gas flow streamline vanished and then inflowing streams entering the Bondi sphere emerged between the horseshoe and shear streamlines in the first and third quadrants in the x–y plane in the midplane of the disk when α = 10⁻². This inflow may promote the accretion of pebbles and increase the accretion probability.

The effects of the turbulence on the accretion probability have been investigated in previous studies. Xu et al. (2017) investigated the effect of turbulent stirring on pebble accretion, performing three types of local 3D hydrodynamical simulations under the shearing-sheet approximation: the ideal magnetohydrodynamical (MHD) simulation (α = 4.5 × 10⁻²), nonideal MHD simulation (α = 6.8 × 10⁻⁴), and pure hydrodynamical simulation (nonturbulence). They measured the accretion rate of particles for an embedded planet with m = 3 × 10⁻³ and 3 × 10⁻² in the Epstein regime. In all cases, these latter authors found that pebble accretion occurs efficiently when St = 0.1–1, but the strong turbulence reduces the accretion probability when St ≲ 0.03 for low-mass planets (m = 3 × 10⁻³). In our study, the reduction of the accretion probability of small pebbles (St ≲0.01) can be seen even for m = 0.03 m = 0.03 (Fig. 11a), which corresponds to a high-mass planet in Xu et al. (2017). This may be the result of resolution in the vicinity of the planet. In Xu et al. (2017), the Bondi radius is resolved by approximately 1–3 cells, while it is resolved by at least 40 cells in the radial direction in our study. Since the typical scale of the perturbed region is the Bondi radius (Kuwahara et al. 2019), a higher resolution may be required to take into account the effects of planet-induced gas flow. However, we emphasize that in contrast to Xu et al. (2017) where the turbulent stirring was directly resolved, our model assumes that the turbulent diffusion balances the tide in the vertical direction. In addition, Xu et al. (2017) assumed η = 0.1 because they focus on the growth of the planet in the outer region of the disk, which may cause the differences in results. We discussed the effect of the headwind in the previous section.

Picogna et al. (2018) investigated the accretion of pebbles onto a planet embedded in a 3D globally isothermal disk. These latter authors modeled two types of disks: the turbulent disk where the turbulence is driven by hydrodynamic instability (vertical shear instability; VSI) and the laminar viscous disk where the particles experience stochastic kicks. The turbulence parameter corresponds to α = 1.737 × 10⁻³. They measured the accretion rate of particles in an equilibrium situation in the Epstein regime. The accretion probability of pebbles with St = 10⁻⁴–1 reaches ~ 1–10% and has a peak at St ~ 10⁻² for the planets at 5.2 au with 5 M_⊕ (m = 0.12) and 10 M_⊕ (m = 0.24) both in VSI turbulent and in laminar viscous disks. As shown in Figs. 11b and c, the accretion probabilities have peaks at St ~ 0.03 and fall in the range of ~1–10% when α ~ 10⁻³, which is consistent with Picogna et al. (2018). However, the accretion probability for a pebble with St = 10⁻³ still has the order of ~ 10% in the viscous disk in Picogna et al. (2018), while it has only ~1% in our study. In contrast to our study, Picogna et al. (2018) consider the random motion of the pebbles due to turbulence, which may lead to the high accretion probability of the smaller pebbles.

4.2.4 Effect of accretion heat

When the planet accretes solid materials, accretion heating is activated. Given the high accretion luminosity and the baroclinic fluid, spiral-like streamlines rise vertically above the growing planet and form an outflow column escaping the Hill sphere after the gas coming out from between the horseshoe and shear region in the first and third quadrants in the x–y plane of the midplane (Chrenko & Lambrechts 2019). In contrast to our results, such upward gas flow may suppress the accretion of pebbles coming from high latitudes.

4.3 Implications for the formation of planetary systems

We propose a formation scenario of planetary systems. Firstly we summarize two parameters, namely the turbulence parameter and the size distribution of the solid materials in a disk, and then give the global structure of the disk in terms of these parameters (Sects. 4.3.1 and 4.3.2). Finally we introduce our scenario based on our results (Sect. 4.3.3).

4.3.1 Turbulence in a protoplanetary disk

Pebble accretion depends on the turbulence strength in a disk. One of the possible origins of the disk turbulence is the magneto-rotational instability (MRI, Balbus & Hawley 1991). The MRI can develop in a region where the ionization degree of the gas is high and produces strong turbulent intensities of α ~ 10⁻³–10⁻². However, Ohmic diffusion would suppress the MRI, which creates an MRI-inactive region (dead-zone; Gammie 1996). The dead-zone lies between a few times 0.1 au and about a few tens of au in a disk, in which the α value should be smaller, α ≲ 10⁻⁴.

Even if the ionization degree is low, the turbulence is generated through hydrodynamic instabilities: vertical shear instability (VSI, Urpin & Brandenburg 1998; Urpin 2003; Nelson et al. 2013; Barker & Latter 2015; Lin & Youdin 2015), convective overstability (COV, Klahr & Bodenheimer 2003; Klahr & Hubbard 2014; Lyra 2014; Latter 2016), subcritical baroclinic instability (SBI, Klahr & Bodenheimer 2003; Petersen et al. 2007; Lesur & Papaloizou 2010), and zombie vortex instability (ZVI, Marcus et al. 2013, 2015). These instabilities develop in different regions of the disk. Assuming a MMSN disk model, the ZVI can develop in the inner region of the disk, ≲ 1 au, where the gas is more adiabatic, and produce α ~ 10⁻³; the COV can develop in the region 1–10 au and trigger a nonlinear state of vortex amplification coined SBI, which induces α ~ 10⁻⁴–10⁻²; the VSI can develop in the outer region of the disk, >10 au, where there tends to be isothermal and generate weak turbulence, α ~ 10⁻⁴ (Malygin et al. 2017, and references therein). Self-regulation may limit the turbulence strength produced by ZVI to as little as α ~ 10⁻⁵.

The turbulence strength beyond a few tens of au can be constrained observationally. To explain the signature of the millimeter-wave polarization and the ring structure of the HL Tau disk shown by Atacama Large Millimeter/submillimeter Array (ALMA) observations (e.g., ALMA Partnership 2015), the maximum dust size should be ~ 150 μm (St ~ 10⁻⁴–10⁻³), and the dust grains are settled to the midplane of the disk (Kataoka et al. 2016; Okuzumi & Tazaki 2019). This requires that the α value be low, α ≲10⁻⁴–10⁻³. Though we adopt the small dust size suggested from the polarization observations of ALMA in the following discussion, we note that it has been argued that a certain proportion of the dust has grown up to > 10 cm based on the Very Large Array (VLA) image at 1.3 cm (Greaves et al. 2008). Dullemond et al. (2018) considered dust trapping in radial pressure bumps as a mechanism to produce dust rings in five different systems observed by the ALMA large program DSHARP campaign. The constrained α∕St differs in each ring but it is larger than ~10⁻². The relatively weak level of turbulence suggested from the observations is also preferable for the formation of ring structures by secular gravitational instability (Takahashi & Inutsuka 2014, 2016; Tominaga et al. 2019).

The possible origin of the turbulence and the α values vary according to the previous studies. Considering the distribution of the turbulence strength in a disk, we may constrain the types of planets that are finally formed. Growth of the planets in the strong turbulence region is expected to be slower than in the weakturbulence region (Figs. 10 and 11). However, since the accretion probability changes significantly depending on the Stokes number, both α and St are needed to discuss planet formation. We discuss the size distribution of the solid materials in a disk in the following section.

4.3.2 Size distribution of pebbles

The planet-induced gas flow significantly reduces the accretion probability of the smaller pebbles (St ≲ 10⁻²) when pebble accretion occurs in 2D or the planetary mass is small, m = 0.03 (Figs. 9–11). The size distribution of the solid materials in a disk may also regulate the types of the planets that are finally formed.

Okuzumi & Tazaki (2019) has shown that the dust particles covered with the H₂ O-ice mantle in the HL Tau disk can grow up to ~50 mm (St ~ 5 × 10⁻²) between the H₂O and CO₂ snow lines. Inside the H₂O snow line, where the dust grains exist as nonsticky silicate grains, the maximum dust size is ~1 mm (St ~ 10⁻³). Outside the CO₂ snow line, since the CO₂ ice is as nonsticky as silicate grains (Musiolik et al. 2016a,b), the upper limit of the size of the dust is ~0.1–1 mm (St ~ 10⁻³–10⁻²).

Given such a size distribution, the growth of the planets inside the H₂ O snow line may be expected to be slower than that between the H₂O and CO₂ snow lines if we consider the reduction of accretion probability of small pebbles due to the planet-induced flow.

4.3.3 Planet formation via pebble accretion in the planet-induced gas flow

Now we introduce the formation scenario of planetary systems based on our results and above discussions (Fig. 13). We adopted three hydrodynamic instabilities as the possible origins of the turbulence: ZVI, COV, and VSI. We divided the disk into three sections according to the previous studies and assumed turbulence strength in each section as: α ~10⁻⁵ (≲ 1 au), α ~10⁻³ (~ 1–10 au), and α ~ 10⁻⁴ (≳ 10 au) (Malygin et al. 2017; Lyra & Umurhan 2019). The outermost region (> 100 au) may be MRI-active, but here we set α ~ 10⁻⁴ from the standpoint of observational constraints. Given the size distribution of the solid materials in a disk (Okuzumi & Tazaki 2019), we assumed that the pebbles have St ~ 10⁻³ (≲1 au), St ~ 3 × 10⁻² (~1–10 au), and St ~ 3 × 10⁻³ (≳10 au), which determine the pebble accretion regime. The pebble accretion regime changes according to the pebble size and the distance from the central star (Eqs. (11), (12), (15)). Even in the region around ~ 1 au, pebble accretion occurs in the Epstein regime if St ≤ 10⁻² (Eq. (15)).

Here we consider the growth of the protoplanet with m = 0.03, which is the lower limit of planetary mass used in our study. Since we only study the shear regime of pebble accretion, we do not consider theformation of proto-cores or proto-planets via pebble accretion. We propose a possible scenario for the formation of the planetary system as follows:

• 1

  The rocky terrestrial planets or super-Earths are formed in the inner region of the disk, ≲1 au, where α and pebble size are small. In the region around ~1 au, since the accretion occurs in the Epstein regime (Eq. (15)), the accretion probability is P_acc ~ 3 × 10⁻⁴ when we assume that α = 10⁻⁵ and St = 10⁻³ (Fig. 11a). The achieved accretion probability is too small to grow the planets to much larger sizes within the typical lifetime of the gas disk (~3–10 Myr). In this case we would expect that the planets that are finally formed are the rocky terrestrial planets. Within the growth track, they may experience inward migration or giant impacts. In the region very close to the star (≲ 0.1 au), pebble accretion occurs in the Stokes regime. In this case we would have P_acc ~ 10⁻² when we assume α = 10⁻⁵ and St = 10⁻³ (Fig. 10a). The accretion regime is switched during inward migration, which may lead to the formation of more massive planets like super-Earths.

• 2

  The gas giants are formed in the middle region of the disk, ~ 1–10 au, where α and pebble size becomes larger than those in the inner region. In this region, since the accretion occurs in the Epstein regime (Eq. (15)), the accretion probability is P_acc ~ 3 × 10⁻² when we assume α = 10⁻³ and St = 3 × 10⁻² (Fig. 11a). As the planets grow, the accretion probability becomes larger. We expect that the planets might exceed the critical core mass (>10 M_⊕) within the typical lifetime of the disk.

• 3

  The ice giants are formed in the outer region. Since the accretion occurs in the Epstein regime (Eq. (15)), the accretion probability is P_acc ~ 9 × 10⁻³ when we assume α = 10⁻⁴ and St = 3 × 10⁻³ (Fig. 11a). With the accretion of a moderate amount of icy pebbles, the planets would eventually become ice giants without evolving into gas giants.

Our scenario may be a natural explanation for the occurrence rate of the exoplanets: an abundance of super-Earths at < 1 au (Fressin et al. 2013; Weiss & Marcy 2014) and a possible peak in the occurrence of gas giants at ~ 2–3 au (Johnson et al. 2010; Fernandes et al. 2019), as well as the architecture of the Solar System.

The dichotomy between the inner super-Earths and outer gas giants may be the result of the reduction of the pebble isolation mass (Eq. (A.10)) in the inner region of the inviscid disk (Fung & Lee 2018). However, the small pebbles (St ≲ 10⁻³) may not be trapped at the local pressure maxima, and they continue to contribute to the growth of the planet (Bitsch et al. 2018). Since the accretion probability of small pebbles decreases significantly in the inner region of the disk, our scenario has the potential to explain the dichotomy even if the influx of small pebbles does not stop at the pressure maxima.

The spacial variety of turbulence strength and pebble size may induce the formation of planetary seeds. Grains would accumulate where the turbulence strength or pebble size changes and the pile-up may trigger planetesimal formation via direct growth and streaming instability (Kretke & Lin 2007; Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007; Morbidelli et al. 2015; Ida et al. 2016).

Fig. 13 Formation scenario of a planetary system. Top: schematic picture of the protoplanetary disk. The brown and white filled circles denote the planets and drifting pebbles. Bottom: assumed Stokes numbers, the origins of the turbulence (ZVI, COV, and VSI), and the resulting accretion probability in the Stokes (solid line) and in the Epstein (dashed line) regimes.

5 Conclusions

We investigated the influence of the 3D planet-induced gas flow on pebble accretion. We considered nonisothermal, inviscid gas flow and performed a series of 3D hydrodynamical simulations on a spherical polar grid that has a planet placed at its center. We then numerically integrated the equation of the motion of pebbles in 3D using hydrodynamical simulation data, which included the Coriolis, tidal, two-body interaction and gas drag forces. Three-types of orbital calculation of pebbles are conducted in this study: in the unperturbed shear flow (Shear case), in theplanet-induced gas flow in the Stokes regime (PI-Stokes case), and in the planet-induced gas flow in the Epstein regime (PI-Epstein case). The subject of the range of dimensionless planetary mass in this study is m = 0.03–0.3, which corresponds to a size ranging from three Mars masses to a super-Earth-sized planet, M_pl = 0.36–3.6 M_⊕, orbiting a solar-mass star at 1 au. Since the planetary masses are larger than the transition mass, we only considered the shear regime of pebble accretion. We summarize our main findings as follows.

• 1

  The trajectories of pebbles in the planet-induced gas flow differ significantly from those in the unperturbed shear flow, in particular for the pebbles with St ≤ 10⁻¹ when m = 0.1 (Fig. 3). Most pebbles coming from within the vicinity of the planetary orbit move away from the planet along the horseshoe flows, which is consistent with a previous study (Popovas et al. 2018). The outflow of the gas deflects the pebble trajectories and inhibits small pebbles from accreting. The 3D trajectories of pebbles also show a similar trend as seen in 2D because the horseshoe streamlines extend in the vertical direction, keeping their configuration, and the vertical scale of the outflow is ~ R_Bondi (Fig. 5). The region that is perturbed by the gravity of the planet can be scaled by the Bondi radius. Therefore, even if the pebbles with a certain Stokes number are hardly influenced by the planet-induced gas flow when the planetary mass is small, thesame-sized pebbles are highly influenced when the planetary mass is large (Fig. 4).

• 2

  The width of the accretion window (Eq. (16)) and the accretion cross section (Eq. (17)) in the planet-induced gas flow becomes smaller than those in the unperturbed shear flow (Figs. 6 and 8). The reduction of these quantities becomes more prominent when the Stokes number is small, or when in the Epstein regime rather than the Stokes regime.

• 3

  The accretion probability of pebbles (Eq. (18)) in the PI-Stokes case matches that in the Shear case when St ≳ 3 × 10⁻³–10⁻², or decreases when St falls below the preceding value in 2D. On the contrary, the accretion probability matches or is slightly larger than that in the Shear case when m ≥ 0.1, or smaller when m = 0.03 for St ≲ 3 × 10⁻² in 3D. In the planet-induced gas flow, the accretion cross section shifts to the right as a whole (Fig. 7), which leads to the increase of the relative velocity between the pebbles and the planet. The reduction of the accretion cross section and the increase of relative velocity cancel each other out. Thus, the suppression of pebble accretion in the PI-Stokes case can be seen only when the Stokes number is small, St ~ 10⁻³, in 2D or the planetary mass is small, m = 0.03, in 3D. In the Epstein regime, however, since the accretion cross section becomes more significant than those in the PI-Stokes case, the increase of the relative velocity does not fully offset its reduction. Therefore, the accretion probability in the PI-Epstein case becomes smaller than that in the Shear case, both in 2D and in 3D and regardless of assumed St and m.

Assuming the global structure of the disk in terms of the distribution of the turbulence parameter and the size distribution of the solid materials in a disk based on previous studies (Malygin et al. 2017; Lyra & Umurhan 2019; Okuzumi & Tazaki 2019), we proposed a formation scenario of planetary systems (Fig. 13). The 3D planet-induced gas flow does affect the pebble accretion and may be helpful to explain the distribution of exoplanets (the dominance of super-Earths at < 1 au (Fressin et al. 2013; Weiss & Marcy 2014) and a possible peak in the occurrence of gas giants at ~ 2–3 au (Johnson et al. 2010; Fernandes et al. 2019), as well as the architecture of the Solar System.

Acknowledgements

We would like to thank an anonymous referee for constructive comments. We thank Athena++ developers: James M. Stone, Kengo Tomida, and Christopher White. This study has greatly benefited from fruitful discussions with Satoshi Okuzumi. Numerical computations were in part carried out on Cray XC30 at the Earth-Life Science Institute and on Cray XC50 at the Center for Computational Astrophysics at the National Astronomical Observatory of Japan. This work was supported by JSPS KAKENHI Grant number 17H01175, 17H06457, 18K13602, 19H01960, 19H05072.

Appendix A Analytical estimation of pebble accretion

A.1 Width of accretion window and accretion cross section

The width of the accretion window in the unperturbed shear flow is expressed by $$b_x=b_{x,0}\exp \left[-{\left(\frac{\text{St}}{2}\right)}^{0.65}\right],$$(A.1)

where the exponential factor is the cutoff parameter given by Ormel & Kobayashi (2012) and b_x,0 is the solution to (Ormel & Klahr 2010) $$b_{x,0}^3+\frac{2}{3}{v}_{\text{hw}}{b}_{x,0}^2-\frac{8}{3}m\text{St}=0,$$(A.2)

which can be divided into two formulas as follows (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Guillot et al. 2014; Ida et al. 2016; Sato et al. 2016): $b_{x,0}\backsimeq \{\begin{array}{c}2\sqrt{\frac{\text{St}}{{\tau }_{\text{B}}}}{R}_{\text{Bondi}},\left(M_{\text{pl}}<M_{\text{t}}:\text{headwind regime}\right)\\ 2{\text{St}}^{1/3}{R}_{\text{Hill}},\left(M_{\text{pl}}>M_{\text{t}}:\text{shear regime}\right)\end{array}$

where τ_B is the gas-free core-crossing time associated with the Bondi radius, ${\tau }_{\text{B}}=\frac{R_{\text{Bondi}}}{v_{\text{hw}}}.$(A.5)

Again, we only considered the shear regime in this study, since the masses of the embedded planets (Eq. (2)) are much larger than the transition mass (Eq. (5)).

We derived the height of the accretion window with a similar procedure of the derivation of b_x,0 in Ormel & Klahr (2010). The gravitational force acting on the pebble at z = b_z,o is $F_{\text{grav}}=m/b_{z,0}^2$. In the shear regime, the approach velocity of pebbles can be written by v_p,∞ = −3∕2xe_y, which gives the interaction timescale as t_int = b_z,0∕v_p,∞. Gas drag can be neglected if t_int < St, but otherwise cannot. When t_int > St the pebble velocity equilibrates to v_p = F_dragSt. According to Ormel & Klahr (2010), we equate v_p with v_p,∞∕4, and then obtain $$b_{z,0}=\sqrt{\frac{8m\text{St}}{3x}}.$$(A.6)

Applying the same cutoff parameter in Eq. (A.1), the height of the accretion window is described by $$b_z=b_{z,0}\exp \left[-{\left(\frac{\text{St}}{2}\right)}^{0.65}\right].$$(A.7)

From Eqs. (A.1) and (A.7), we calculate the accretion cross section by $$A_{\text{acc}}=2{\displaystyle {\int }_0^{b_z}{\displaystyle {\int }_0^{b_x}\text{d}}}x\text{d}z=16{\text{St}}^{2/3}{R}_{\text{Hill}}^2\exp \left[-\frac{3}{2}{\left(\frac{\text{St}}{2}\right)}^{0.65}\right].$$(A.8)

Equations (A.1) and (A.8) are plotted in Figs. 6 and 8 with dashed-dotted lines. Both equations show excellent agreement with the numerically calculated width of the accretion window and accretion cross section in the Shear case.

A.2 Pebble isolation mass

When a growing planet opens up a gap in a disk and generates a pressure maxima outside of its orbit, pebble accretion halts (Lambrechts et al. 2014). This happens when the core mass reaches pebble isolation mass: $$M_{\text{iso}}^{L14}\approx 20{\left(\frac{a}{5\text{au}}\right)}^{3/4}{M}_{\oplus }.$$(A.9)

A subsequent study derived a detailed pebble isolation mass: $$M_{\text{iso}}=25f_{\text{fit}}{M}_{\oplus }+\frac{{\Pi }_{\text{crit}}}{\lambda }{M}_{\oplus },$$(A.10)

where λ ≈ 0.00476∕f_fit, Π_crit = α∕2St, and $$f_{\text{fit}}={[\frac{H/a}{0.05}]}^3[0.34{(\frac{3}{\log (\alpha )})}^4+0.66][1-\frac{\frac{\partial \ln P}{\partial \ln a}+2.5}{6}],$$(A.11)

(Bitsch et al. 2018).

A.3 Accretion rate

The accretion rate of pebbles, Ṁ_p, is given analytically by (Ormel 2017): ${\stackrel{.}{M}}_{\text{p}}=\{\begin{array}{c}{\stackrel{.}{M}}_{2\text{D}},\left(2\text{D},\text{analytical}\right)\\ {\stackrel{.}{M}}_{2\text{D}}\frac{2b_x}{b_x+H_{\text{p}}\sqrt{8/\pi }},\left(3\text{D},\text{analytical}\right)\end{array}$

where Ṁ_2D is calculatedby $${\dot{M}}_{\text{2D}}=2{\displaystyle {\int }_0^{b_x}{\Sigma }_{\text{p}}}{v}_{\text{p},\infty }\text{d}x.$$(A.14)

We use a slightly different formula of 3D accretion rate described in Ormel et al. (2017), where the factor of two is not included. Equations (A.12) and (A.13) are plotted in Fig. 9 with dashed-dotted lines. Both equations show excellent agreement with the numerically calculated accretion rate in the Shear case.

Appendix B Interpolation of simulation data

In a series of numerical results obtained from our grid simulations using Athena++ code, all physical quantities were given as discrete data on the center of each grid. In this study, we interpolated the gas velocity and the density at the position using the bilinear interpolation method. We note that when we interpolate the gas velocity and the density, since the hydro simulations are performed in a spherical polar coordinates, we transformed the position of the pebble in Cartesian coordinates to that in spherical coordinates. Given the position of the pebble at an arbitrary time, $r=(x,y,z)=\left(r(t),\theta (t),\varphi (t)\right)$, a certain physical quantity at the position of the pebble, q(r, t), is described as the following interpolation formula: $$q(r,t)=(1-\zeta )q_{\text{itp},1}+\zeta q_{\text{itp},2},$$(B.1)

where q_itp,1 and q_itp,2 are the interpolated physical quantities on the top and bottom surface of the cell including the pebble, $$\begin{array}{ll}{q}_{\text{itp},1}=\hfill & \left(1-\eta \right)\left(1-\xi \right)q_{i,j,k}+\eta \left(1-\xi \right)q_{i+1,j,k}\hfill \\ \hfill & +\eta \xi q_{i+1,j,k+1}+\left(1-\eta \right)\xi q_{i,j.k+1},\hfill \\ q_{\text{itp},2}=\hfill & \left(1-\eta \right)\left(1-\xi \right)q_{i,j+1,k}+\eta \left(1-\xi \right)q_{i+1,j+1,k}\hfill \\ \hfill & +\eta \xi q_{i+1,j+1,k+1}+\left(1-\eta \right)\xi q_{i,j+1,k+1},\hfill \end{array}$$

where q_i,j,k is a certain physical quantity at the center of the grid, the subscripts denote the grid number, and ζ, η, ξ are given by $$\zeta =\frac{\theta (t)-{\theta }_j}{{\theta }_{j+1}-{\theta }_j},\eta =\frac{r(t)-r_i}{r_{i+1}-r_i},\xi =\frac{\varphi (t)-{\varphi }_k}{{\varphi }_{k+1}-{\varphi }_k}.$$(B.4)

After the interpolation, we integrate Eq. (9) in Cartesian coordinates. We then have the gas velocity after conversion from spherical polar coordinates to Cartesian coordinates using Athena ++.

Appendix C Independence of accretion rate from orbital radius

From Eqs. (3), (4), (19), and (20), the surface density of pebbles can be expressed by $${\Sigma }_{\text{p}}=\frac{{\dot{M}}_{\text{disk}}}{2\pi \times \frac{\text{d}\ln P}{\text{d}\ln a}}\frac{1+{\text{St}}^2}{\text{St}}\frac{1}{H^2\Omega }.$$(C.1)

When P ∝ a^−q, where q is an arbitrary index, which holds for typical disk models including the MMSN model, the pressure gradient, d ln P∕dlna, on the right-hand side of Eq. (C.1) becomes constant. Therefore, $${\Sigma }_{\text{p}}\propto H^{-2}{\Omega }^{-1}.$$(C.2)

From Eqs. (21), (22), and (C.1), the pebble density can be described by $${\rho }_{\text{p}}\propto \frac{{\Sigma }_{\text{p}}}{H}\propto H^{-3}{\Omega }^{-1}.$$(C.3)

In this study, the length scale is normalized by the disk scale height, and therefore the integrals in Eqs. (23), (24) are proportional to H and H², respectively, when the dimensionless planetary mass, m, is fixed. Thus, from Eqs. (23), (24), (C.3), and (C.4), the accretion rates in both 2D and 3D are proportional to $$\begin{array}{l}{\dot{M}}_{\text{p},\text{2D}}\propto {\Sigma }_{\text{p}}{v}_{\text{p},\infty }H\propto H^{-2}{\Omega }^{-1}\times H^2\Omega =\text{const}.,\\ {\dot{M}}_{\text{p},\text{3D}}\propto {\rho }_{\text{p}}{v}_{\text{p},\infty }{H}^2\propto H^{-3}{\Omega }^{-1}\times H^3\Omega =\text{const}.\end{array}$$

The accretion rates both in 2D and 3D are independent of the orbital radius, a, for a fixed dimensionless planetary mass, m.

References

All Tables

All Figures

Fig. 1 Schematic picture of orbital calculation of pebbles. A planet is located at the origin of the co-rotational frame. The dashed line represents the outer boundary of the hydrodynamical simulations. The starting point of orbital calculation is beyond the outer boundary of hydrodynamical simulations. Its y-component is fixed at 40RHill. The initial velocity of the pebble is the same with the Keplerian shear velocity. The gas velocity is assumed to be the speed of the Keplerian shear both inside and outside of rout in the shear flow case in the Stokes regime (Shear case), but is switched to the gas velocity obtained from the hydrodynamical simulations within rout in the planet-induced flow case in the Epstein regime (PI-Epstein case), and in the Stokes regime (PI-Stokes case).

In the text

	Fig. 2 Streamlines of 3D planet-induced gas flow around the planet with m = 0.1 at t = 150. The red, green, and blue solid lines are the recycling streamlines, the horseshoe streamlines, and the Keplerian shear streamlines, respectively. The sphere is the Bondi sphere of the planet.
In the text

Fig. 3 Trajectories of pebbles in Shear case (top) and PI-Stokes case (bottom) with different Stokes numbers at midplane around an embedded planet with m = 0.1. We set zs = 0 for all cases. The red and blue solid lines correspond to the trajectories of pebbles which accreted and did not accrete onto the planet, respectively. The dashed and dotted circles show the Hill and the Bondi radius of the planet, respectively. The black dots at the center of each panel denote the position of the planet. The interval of pebbles at their initial locations is 0.05 [H].

In the text

Fig. 4 Trajectories of pebbles in the PI-Stokes case with St = 10−2 at the midplane around the embedded planets with m = 0.03 (left panel) and 0.3 (right panel). We set zs = 0. The red and blue solid lines correspond to the trajectories of pebbles which accreted and did not accrete onto the planet, respectively. The dashed and dotted circles show the Hill and the Bondi radius of the planet, respectively. The sizes of the Hill radius are 0.22 [H] for m = 0.03 and 0.46 [H] for m = 0.3. The black dots at the center of each panel denote the position of the planet. The interval of pebbles at their initial locations is 0.05 [H].

In the text

Fig. 5 Trajectories of pebbles with St = 10−2 around an embedded planet with m = 0.1 in the PI-Stokes case. The height of the initial position of the pebbles is zs = 0.7RHill. The red and blue solid lines correspond to the trajectories of pebbles which accreted and did not accrete onto the planet, respectively. The black dot denotes the position of the planet. The sphere around the planet represents the Hill radius of the planet. We only plot the trajectories within the region where r < 10RHill. The interval of pebbles at their initial locations is 0.05 [H].

In the text

	Fig. 6 Width of accretion window in the midplane, wacc(0), as a function of the Stokes number in the PI-Stokes case (solid lines), PI-Epstein case (dashed-lines), and the Shear case (dotted lines). The dashed-dotted lines correspond to the analytical estimation for the Shear case expressed by Eq. (A.4). Colors indicate the mass of the planet: m = 0.03 (red), m = 0.1 (yellow), and m = 0.3 (blue).
In the text

Fig. 7 Accretion cross section with the different Stokes numbers for the planet with m = 0.1 in the Shear case (top) and the PI-Stokes case (bottom). We assumed α = 10−3 and the dust-to-gas ratio is equal to 10−2. Color represents the density of the pebbles expressed by Eq. (21) normalized by the gas density. The two panels in panels a and e show the enlarged outlines of accretion cross sections. We note that the color contour is saturated for the ρp ≲ 10−3.

In the text

Fig. 8 Accretion cross section as a function of the Stokes number as a function of the Stokes number in the PI-Stokes case (solid lines), PI-Epstein case (dashed-lines), and the Shear case (dotted lines). The dashed-dotted lines correspond to the analytical estimation for the Shear case expressed by Eq. (A.8). Colors indicate the mass of the planet: m = 0.03 (red), m = 0.1 (yellow), and m = 0.3 (blue).

In the text

Fig. 9 Accretion rate (left vertical axis) and probability (right vertical axis) as a function of the Stokes number in the PI-Stokes case (solid lines), PI-Epstein case (dashed-lines), and the Shear case (dotted lines). The dashed-dotted lines correspond to the analytical estimation for the Shear case expressed by Eqs. (A.12) and (A.13). Top: 2D case. Bottom: 3D case (α = 10−3). Colors indicate the mass of the planet: m = 0.03 (red), m = 0.1 (yellow), and m = 0.3 (blue).

In the text

	Fig. 10 Dependence of accretion rate (left vertical axis) and accretion probability (right vertical axis) on the turbulent parameter, α. Panels a–c: results of the planet with m = 0.03, 0.1, and 0.3 in the PI-Stokes case (solid lines) and Shear case (dotted lines).
In the text

	Fig. 11 Same as Fig. 10, but the dashed lines show the accretion probability in the PI-Epstein case.
In the text

	Fig. 12 Accretion probability as a function of the planetary mass and the Stokes number for the turbulence parameters α = 10−6, 10−4, and 10−2 (bottom to top) in the Shear case (left column), the PI-Stokes case (middle column), and the PI-Epstein case (right column). The contours represent the accretion probabilities.
In the text

	Fig. 13 Formation scenario of a planetary system. Top: schematic picture of the protoplanetary disk. The brown and white filled circles denote the planets and drifting pebbles. Bottom: assumed Stokes numbers, the origins of the turbulence (ZVI, COV, and VSI), and the resulting accretion probability in the Stokes (solid line) and in the Epstein (dashed line) regimes.
In the text