Turbulence plays a key role in the transport of pebblesized particles. It also affects the ability of pebbles to be accreted by protoplanets because it stirs pebbles out of the disk midplane. In addition turbulence suppresses pebble accretion once the relative velocities become too high for the settling mechanism to be viable. Following Paper I we aim to quantify these effects by calculating the pebble accretion efficiency  using threebody simulations. To model the effect of turbulence on the pebbles we derive a stochastic equation of motion SEOM applicable to stratified disk configurations. In the strong coupling limit ignoring particle inertia the limiting form of this equation agrees with previous works. We conduct a parameter study and calculate  in 3D varying pebble and gas turbulence properties and accounting for the planet inclination. We find that strong turbulence suppresses pebble accretion through turbulent diffusion agreeing closely with previous works. Another reduction of  occurs when the turbulent rms motions are high and the settling mechanism fails. In terms of efficiency the outer disk regions are more affected by turbulence than the inner regions. At the location of the HSUB2SUBO iceline planets around lowmass stars achieve much higher efficiencies. Including the results from Paper I we present a framework to obtain  under general circumstances.
['planets and satellites: formation', 'protoplanetary disks', 'methods: numerical', 'Astrophysics - Earth and Planetary Astrophysics']
Catching drifting pebbles. II. A stochastic equation of motion for pebbles
{'II. A stochastic equation of motion for pebbles': 'Chris W. Ormel, Beibei Liu  \nAnton Pannekoek Institute (API), University of Amsterdam, Science Park 904,1090GE Amsterdam, The Netherlands e-mail: [c.w.ormel,b.liu]@uva.nl  \nApril 11, 2018', 'ABSTRACT': 'Turbulence plays a key role in the transport of pebble-sized particles. It also a GLYPH<11> ects the ability of pebbles to be accreted by protoplanets because it stirs pebbles out of the disk midplane. In addition, turbulence can suppress pebble accretion once the relative velocities become too high for the settling mechanism to be viable. Following Paper I, we aim to quantify these e GLYPH<11> ects by calculating the pebble accretion e GLYPH<14> ciency " using three-body simulations. To model the e GLYPH<11> ect of turbulence on the pebbles, we derive a stochastic equation of motion (SEOM) applicable to stratified disk configurations. In the strong coupling limit (ignoring particle inertia) the limiting form of this equation agrees with previous works. We conduct a parameter study and calculate " in 3D, varying pebble and gas turbulence properties and accounting for the planet inclination. We find that strong turbulence suppresses pebble accretion through turbulent di GLYPH<11> usion, agreeing closely with previous works. Another reduction of " occurs when the turbulent rms motions are high and the settling mechanism fails. In terms of e GLYPH<14> ciency, the outer disk regions are more a GLYPH<11> ected by turbulence than the inner regions. At the location of the H2O iceline, planets around low-mass stars achieve much higher e GLYPH<14> ciencies. Including the results from Paper I, we present a framework to obtain " under general circumstances.  \nKey words. planets and satellites: formation - protoplanetary disks - methods: numerical', '1. Introduction': 'It is widely believed that turbulence plays an important role in the evolution of protoplanetary disks. For a long time the magneto-rotational instability (MRI) (Balbus & Hawley 1991) has been regarded as the leading candidate in driving the disk\'s angular momentum transport. More recently, disk wind models have regained traction (Bai et al. 2016; Suzuki et al. 2016; Gressel 2017), where the turbulence in the midplane regions is limited to hydrodynamic instabilities such as the vertical shear instability (Nelson et al. 2013; Stoll & Kley 2014). Turbulence, in addition, is important in shaping the outcome of the early coagulation process. Already at low mach numbers, turbulence dominates the relative velocity between particles (Völk et al. 1980; Ormel & Cuzzi 2007; Pan & Padoan 2010). It is also the only explanation for why we infer vertical structure (e.g., flared versus settled geometry), since turbulence allows small particles to be lifted from the disk midplane regions. Indeed, with the current state-of-the-art models small (micron-size) particles are produced in the midplane through collisions between pebbles and boulder sized particles before they di GLYPH<11> use upwards (Birnstiel et al. 2010, 2011; Krijt & Ciesla 2016).  \nBecause of its subsonic nature, obtaining observational evidence of turbulence is hard. Disks like TW Hya, HD 163296, and DM Tau have been modeled by several groups (Hughes et al. 2011; Guilloteau et al. 2012; Flaherty et al. 2015, 2017) with turbulent Mach numbers inferred from the rather quiescent GLYPH<24> 0.01 to the more vigorous GLYPH<24> 0 : 1. However, it should be emphasized that constraining the turbulent rms velocity ( GLYPH<27> ) by these single-line profiles requires that the temperature profile be known to great precision. Generally, uncertainties a GLYPH<11> ecting GLYPH<27> are limited by the  \nabsolute flux calibration and spectral resolving power (Teague et al. 2016). More indirect methods of obtaining GLYPH<27> employ the appearance of the dust disk in ALMA imagery. Applied to HL tau this indicates that the pebbles are settled into the midplane, resulting in a vertical turbulent di GLYPH<11> usivity parameter GLYPH<11> z GLYPH<24> 10 GLYPH<0> 4 (Pinte et al. 2016). 1 Flock et al. (2017) conclude this is consistent with a magnetized disk models that feature a \'dead\' midplane.  \nBoth classical planetesimal-driven models for planet formation (Safronov 1969; Pollack et al. 1996) as well as the more recent pebble accretion model (Ormel & Klahr 2010; Lambrechts & Johansen 2012) are greatly a GLYPH<11> ected by turbulence. Both models operate best under low-turbulence conditions. The runaway growth phase for the classical, planetesimal-driven accretion paradigm can only operate once the planetesimals start out with close to zero velocity dispersions, but stochastic forcing by turbulence-triggered density fluctuations (Ida et al. 2008; Nelson & Gressel 2010; Gressel et al. 2011, 2012; Okuzumi & Ormel 2013) excites planetesimals to random velocities higher than their escape velocity. This implies that planetesimals have to be born large or that turbulence has to be weak (Ormel & Okuzumi 2013; Kobayashi et al. 2016). Similarly, the e GLYPH<14> cacy of pebble accretion to grow planets also depends on the turbulence. As pebbles will be stirred away from the midplane, it reduces the number of pebbles left to be accreted (Ormel & Klahr 2010; Guillot et al. 2014; Morbidelli et al. 2015). A second, less known, e GLYPH<11> ect is that turbulent forcing may provide particles with  \nan additional relative motion, which could also suppress accretion.  \nIn this work, we consider simultaneously the e GLYPH<11> ects of turbulent di GLYPH<11> usivity and turbulent velocity. We do this by deriving a stochastic equation of motion (SEOM) for pebble-sized particles. Simply put, the SEOM is an extension of the Newtonian equation of motion, but with an additional stochastic component. For planets, stochastic forces have been invoked as a means to cross mean motion resonances (Rein & Papaloizou 2009; Paardekooper et al. 2013). As detailed by Rein & Papaloizou (2009) the stochastic motion is characterized by two parameters: the di GLYPH<11> usivity DP and the correlation time t corr. The latter is crudely the time over which the stochastic force changes its direction. For pebbles, we adopt a similar model, where now the stochastic motions are driven by aerodynamical coupling to the turbulent gas. However, in the few studies that have considered stochastic e GLYPH<11> ects for pebble-sized particles, it is often assumed that turbulence does not feature a correlation time, i.e., t corr is assumed less than any other timescale in the problem (Ciesla 2010; Zsom et al. 2011; Krijt & Ciesla 2016). This implies \'white noise\' behavior, i.e., that the particle is displaced in a random direction at every time. This approximation is known as the strong coupling limit (SCA).  \nA key goal of this paper is to test how the SCA fares in the light of the more accurate SEOM. We find that the SCA is generally applicable, as long as both the particle stopping time t stop and the turbulent correlation time are su GLYPH<14> ciently small. In addition, we will study the e GLYPH<11> ect of a vertically varying turbulent gas di GLYPH<11> usivity, Dzz ( z ), to investigate when it is viable to stir a fraction of pebble size particles to the disk surface.  \nOur main thrust will be to apply the SEOM and the SCA methods to calculate pebble accretion e GLYPH<14> ciencies in threedimensional (3D) settings. In Liu & Ormel (2018, henceforth Paper I) we have defined " as the probability that a pebble, drifting towards the star, will be accreted by a single planet(esimal) 2 . A very small value of " implies that a large number of pebbles are needed to grow the planet, while " close to unity implies that pebble accretion is a very e GLYPH<14> cient accretion process. In Paper I we used planar 3-body calculation (star, planet, pebble) to calculate " in two dimensions. We then investigated how this " 2D changed as function of planet properties (mass and eccentricity), disk properties (position, radial drift velocity), and pebble properties (stopping time). In this work, we extend these calculation to the vertical dimension by including the planet\'s inclination and disk turbulence. With the " 3D of this paper and the " 2D of Paper I, we then obtain a general recipe for the pebble accretion e GLYPH<14> ciency ( " ) of a single planet.  \nThe plan of the paper is the following. In Sect. 2 we derive the SEOM. This section, as well as Appendix A, are rather technical and may be skipped by readers more interesting in the physical applications. In Sect. 3 we apply our newly developed SEOMto find vertical density distributions and show that our results are consistent with previous numerical and analytical studies. We apply the SEOM and SCA to pebble accretion in Sect. 4. We find " for a variety of settings (planet mass and inclination, particle and disk properties) and present a framework to obtain " under general circumstances (including the results found in Paper I). A comparison with previous studies is presented in Sect. 5. We summarize our findings in Sect. 6.', '2.1. Advection-diffusion equation': 'In this work we model turbulence motion of particles and gas by an advection-di GLYPH<11> usion equation  \n@GLYPH<26> P @ t + r GLYPH<1> GLYPH<26> P v GLYPH<0> r GLYPH<1> GLYPH<26> gas D P r GLYPH<26> P GLYPH<26> gas ! = 0 (1)  \nwhere GLYPH<26> P is the density of particles or gas species, v the systematic (drift) velocity, GLYPH<26> gas the gas density, and D P the particle diffusivity tensor whose elements are denoted Dij . Importantly, the di GLYPH<11> usion term acts on the gradient of the concentration ( GLYPH<26> P =GLYPH<26> gas): it tends to erase concentration gradients and vanishes when the concentration is uniform.  \nIn this work, we will restrict di GLYPH<11> usion to operate only in the vertical ( z ) direction, considering only D P ; zz. Furthermore, we adopt the vertically isothermal solution for the gas density  \nGLYPH<26> gas = GLYPH<6> gas H gas p 2 GLYPH<25> exp 2 6 6 6 6 6 4 GLYPH<0> 1 2 z H gas ! 2 3 7 7 7 7 7 5 (2)  \nwhere GLYPH<6> gas is the gas surface density and H gas the pressure scaleheight. Under these conditions Eq. (1) can be manipulated  \n@GLYPH<26> P @ t + @ @ z 0 B B B B @ v z GLYPH<0> D P ; zz z H 2 gas 1 C C C C A GLYPH<26> P = @ @ z D P ; zz @GLYPH<26> P @ z (3)  \n(Ciesla 2010) 3 . For small particles (including pebbles) the vertical velocity v z equals the settling velocity, v z = GLYPH<0> z GLYPH<10> 2 t stop with t stop the stopping time and GLYPH<10> the Keplerian orbital frequency. The gas density no longer appears in Eq. (3), but the di GLYPH<11> usivity appears at two places. On the RHS the di GLYPH<11> usive term is responsible for spreading the particle concentration, resulting in a broader distribution of GLYPH<26> P ( z ). However, the additional advection term GLYPH<0> D P ; zzz = H 2 gas - a consequence of imposing Eq. (2) - counteracts this, enforcing the particle layer to remain stratified with a finite dispersion at all times (Ciesla 2010).  \nFor a distribution of particles GLYPH<26> P ( z ) dz = GLYPH<6> gives the fraction of the particles within the interval [ z ; z + dz ]. For a single particle, P ( z ) = GLYPH<26> ( z ) = GLYPH<6> similarly denotes the probability of finding the particle within [ z ; z + dz ] where P ( z ) is the probability density. We will use this identification below to obtain the correct, statistical properties of our single-particle (Lagrangian) stochastic model.', '2.2. Stochastic equation of motion (SEOM)': 'The stochastic equation of motion is given by the following set of stochastic di GLYPH<11> erential equations (SDEs):  \nd x = v d t (4a)  \nd v = 0 B B B B B @ F g + GLYPH<0> v + v gas + p Dzz = t corr GLYPH<16> t e z + v hs t stop 1 C C C C C A d t (4b)  \nd GLYPH<16> t = GLYPH<0> GLYPH<16> t t corr d t + r 2 t corr d Wt (4c)  \nwhere x is position and v the velocity of a particle. The particle is subject to gravitational forces F g and gas drag forces. The latter have been expressed in terms of the stopping time, GLYPH<1> v = t stop, where GLYPH<1> v is the relative gas-particle velocity. In Eq. (4c) Wt denotes a Wiener process (Brownian motion) and the corresponding di GLYPH<11> erential is d Wt GLYPH<24> p d t N (0 ; 1) where N (0 ; 1) is the normal distribution with zero mean and unity variance.  \nApart from v three velocity terms appear in Eq. (4b):  \n- 1. A laminar gas velocity v gas. In our case, the gas velocity operates in the azimuthal direction  \nv gas = (1 GLYPH<0> GLYPH<17> ) v K e GLYPH<30> (5)  \nwhere v K = p G ( M ? + Mp ) = r and GLYPH<17> represents the disk radial pressure gradient  \nGLYPH<17> = GLYPH<0> 1 2 @ log P @ log r ! midplane H gas r ! 2 (6)  \nis assumed constant (Nakagawa et al. 1986).  \n- 3. A correction term v hs\n- 2. A turbulent velocity v turb. In Eq. (4b) this has been written in terms of an rms value ( GLYPH<27> z = p Dzz = t corr) and a nondimensional stochastic variable GLYPH<16> t . Here t corr and Dzz are, respectively, the correlation time and di GLYPH<11> usivity of the turbulent gas. In Eq. (4b) the turbulent forcing acts only in the vertical dimension.  \nv hs = GLYPH<0> Dzzz H 2 gas e z + 1 2 D 0 zz : (7)  \nwhere D 0 zz = @ D zz =@ z . This is needed to enforce that Eq. (4) satisfies the hydrostatic balance condition, which assumption has entered the advection-di GLYPH<11> usion Eq. (3). It also accounts for spatial gradients in Dzz . We derive it below.  \nFinally, Eq. (4c) is a stochastic di GLYPH<11> erential equation (SDE) about a quantity GLYPH<16> t . This can be thought of as the normalized strength of the turbulent velocity. Specifically, Eq. (4c) describes an Ornstein-Uhlenbeck process (Uhlenbeck & Ornstein 1930) with zero mean ( h GLYPH<16> t i = 0), unity variance ( h GLYPH<16> 2 t i = 1) and correlation time t corr. On long timescales GLYPH<16> t will be normally distributed; events separated by GLYPH<1> t GLYPH<29> t corr are uncorrelated. However, times separated by GLYPH<1> t GLYPH<28> t corr will feature a similar value of GLYPH<16> t and hence a similar turbulent gas velocity.', '2.3. Turbulent correlation time': 'In Eq. (4) we are at liberty to choose t corr, which can be identified with the correlation time (or lifetime) of the turbulent eddies. A smaller t corr (while keeping Dzz fixed) implies a more vigorous turbulent forcing (larger GLYPH<27> z ), while a long t corr implies that the turbulence is characterized by weaker but larger and longer-lived eddies. It is customary to adopt the Shakura & Sunyaev (1973) GLYPH<11> -parameterization for the turbulent viscosity  \nGLYPH<23> T = GLYPH<11> H 2 gas GLYPH<10> : (8)  \nWe will adopt a similar parameterization for the gas di GLYPH<11> usivity, i.e., Dzz = GLYPH<11> zH 2 gas GLYPH<10> where GLYPH<11> z reflects the di GLYPH<11> usivity of the gas, not angular momentum transport. Using Dzz = t 2 corr GLYPH<27> z the turbulent rms velocity becomes  \nGLYPH<27> z = r Dzz t corr = GLYPH<11> 1 = 2 z H gas GLYPH<10> p t corr GLYPH<10> : (9)  \nFollowing Dubrulle et al. (1995), Cuzzi et al. (2001) and Johansen et al. (2006) we usually adopt t corr = GLYPH<10> GLYPH<0> 1 and hence GLYPH<27> z = GLYPH<11> 1 = 2 z H gas GLYPH<10> . In Sect. 3.2 we also consider models where t corr is longer.  \nThe following qualifications will be adopted towards the turbulence strength:  \n- -laminar for GLYPH<11> z = 0;\n- -moderately turbulent for 10 GLYPH<0> 4 < GLYPH<11> z < 10 GLYPH<0> 2 ;\n- -weakly turbulent for GLYPH<11> z < 10 GLYPH<0> 4 ;\n- -strongly turbulent for GLYPH<11> z > 10 GLYPH<0> 2 .', '2.4. Strong coupling approximation (SCA)': 'In the strong coupling approximation t stop is assumed small. It can be shown that Eq. (4) then simplifies with the equation of motion being described by a single SDE  \nd x = " F t stop + v gas + v hs + D 0 zz e z 2(1 + t stop = t corr) # d t + p 2 D zz e z d Wt (10)  \n(see Appendix A for the derivation). In Eq. (10) we have allowed D zz to depend on position, which would give rise to an additional advective term (the fourth term in the square brackets). Equation (10) is analogous to Smoluchowski (1916) equation for the stochastic motion of a massless particle subject to a fluctuating force.  \nWe are now in a position to obtain the hydrostatic correction term v hs. SDEs of the form  \nd x = A ( x )d t + B ( x )d Wt (11)  \ncan equivalently be cast in terms of an equation for the evolution of the probability density P ( x ; t ) - i.e., a Fokker-Planck equation  \n@ P ( x ; t ) @ t + @ @ x A ( x ) P = 1 2 @ 2 @ x 2 B 2 ( x ) P : (12)  \n(van Kampen 1992) 4 . Applied to Eq. (10), the Fokker-Planck equation for the probability density P ( z ; t ) reads  \ndP ( z ; t ) dt + @ @ z Fzt stop + 1 2 D 0 zz + v hs ! P = @ 2 @ z 2 D zz P (13)  \nwhere we consider the limit t stop GLYPH<28> t corr. The RHS of Eq. (13) can be expanded as @ @ z D zz @ @ z P + @ @ z D 0 zz P . Identifying the probability density P ( z ) with the density GLYPH<26> P of Eq. (3), Fzt stop with v z , and using that D P ; zz = D zz for strongly coupled particles (Völk et al. 1980), we obtain the hydrostatic correction term, Eq. (7). With this correction term, the SCA for the particle position in 1D reads  \nd z = 2 6 6 6 6 4 GLYPH<0> z GLYPH<10> 2 t stop GLYPH<0> D zz z H 2 gas + D 0 zz 3 7 7 7 7 5 d t + p 2 D zz d Wt : (14)  \nas was already derived by Ciesla (2010) and also used in Zsom et al. (2011) and Charnoz et al. (2011). Comparing the second and third terms on the RHS, we obtain that the turbulent gradient e GLYPH<11> ect becomes important when D zz changes on scales less than GLYPH<1> z GLYPH<24> D zz = D 0 zz = H 2 gas = z .  \n0  \n6  \nFig. 1. top : Normalized distributions of the vertical position z for the integration of the tracer case, using the strong coupling approximation (SCA) with t stop = 0. The vertical height z is recorded after every 1 GLYPH<10> GLYPH<0> 1 for t = 10 5 GLYPH<10> GLYPH<0> 1 . Bars give the simulated distribution while the analytic normal - distribution is shown by the black dashed line. bottom : temporal evolution of the vertical position for GLYPH<11> z = 10 GLYPH<0> 2 (black) and GLYPH<11> z = 10 GLYPH<0> 2 (blue).  \n<!-- image -->  \nWe reflect on our findings. Equation (10) with v hs equal to Eq. (7) describes the stochastic motion of a particle experiencing drag and turbulent forces, with the turbulence characterized by a correlation time t corr and a (possibly spatially dependent) gas di GLYPH<11> usivity D zz. Under the assumption of small t stop and small t corr we obtain Eq. (14), consistent with Ciesla (2010). However, these equations provide no model for the particle velocity; and they will fail when the SCA-conditions no longer materialize (long t stop or long t corr) - i.e., when the particle\'s inertia matters. In these cases Eq. (4) provides a more general description of stochastic motion of pebble-sized particles.', '3. Vertical diffusion': 'We test our algorithms - the stochastic equation of motion (SEOM; Eq. (4)) and the strong coupling approximation (SCA; Eq. (14)) - for tracer particles ( t stop = 0) in Sect. 3.1 and massive particles (Sect. 3.2) for a variety of stopping times and GLYPH<11> z .', '3.1. Tracer particle': 'Tracer particles should have a vertical distribution identical to the gas, Eq. (2). Because t stop = 0 for tracer particles, the SEOM, as described in Eq. (4), contains a singularity and is not applicable. Therefore, we adopt the SCA method. In Eq. (10) we take  \nt stop = 0, GLYPH<17> = 0, Dz = GLYPH<11> zH 2 gas GLYPH<10> . The choices for H gas and GLYPH<10> are arbitrary.  \nIn Fig. 1a we show the distribution of the vertical position of a single particle for GLYPH<11> z = 10 GLYPH<0> 4 , 10 GLYPH<0> 3 and 10 GLYPH<0> 2 . These have been obtained by storing the vertical positions after every 1 GLYPH<10> GLYPH<0> 1 for a total time of t max = 10 5 GLYPH<10> GLYPH<0> 1 . In addition we plot the expected distribution according to Eq. (2). The distributions are normalized such that they integrate to unity. Clearly, the distributions among the GLYPH<11> z di GLYPH<11> er, with GLYPH<11> z = 10 GLYPH<0> 2 best matching the expected normal distribution GLYPH<11> z = 10 GLYPH<0> 4 the worst. The origin of these di GLYPH<11> erences is the di GLYPH<11> erent number of independent samples that are obtained among the GLYPH<11> z . Since our sampling time is only GLYPH<1> t = 1 GLYPH<10> GLYPH<0> 1 , much smaller than the di GLYPH<11> usion time t di GLYPH<11> = 1 =GLYPH<11> z GLYPH<10> , sequential samples (in time) will be strongly correlated; only t max = t di GLYPH<11> = 10 5 GLYPH<11> z samples will be independent. Hence, the higher GLYPH<11> , the better the correspondence to a Gaussian distribution.', '3.2. Massive particle': "Next we consider the vertical distribution of massive particles ( t stop > 0), obtained by the SEOM and the SCA methods. This means that in Eq. (4) we put F g = GLYPH<0> GLYPH<10> 2 z e z . We further take t corr = 1 GLYPH<10> GLYPH<0> 1 unless mentioned otherwise. Integrations ran for 10 6 GLYPH<10> GLYPH<0> 1 and the sampling period was GLYPH<1> t = 10 GLYPH<10> GLYPH<0> 1 . The results are shown in Fig. 2. We fix GLYPH<11> z = 10 GLYPH<0> 2 but vary GLYPH<28> s = t stop GLYPH<10> = 10 GLYPH<0> 2 (left panel), 10 GLYPH<0> 1 (center), and 10 0 (right).  \nIn all panels we compare the results for the SEOM (gray curves) with the SCA of Eq. (10) (blue). The thin black curve corresponds to a normal distribution with pebble aspect ratio  \nhP = r GLYPH<11> z GLYPH<11> z + GLYPH<28> s GLYPH<24> GLYPH<0> 1 = 2 h gas (15)  \n(Youdin & Lithwick 2007) 5 where h gas = H gas = r and  \nGLYPH<24> = 1 + GLYPH<28> s ( GLYPH<10> t corr) 2 GLYPH<28> s + GLYPH<10> t corr : (16)  \nIn the limit of GLYPH<10> t corr GLYPH<28> 1 or GLYPH<28> s GLYPH<28> ( GLYPH<10> t corr) GLYPH<0> 1 , GLYPH<24> GLYPH<25> 1 and the pebble aspect ratio reduces to hP = p GLYPH<11> z = ( GLYPH<11> z + GLYPH<28> s ) h gas (Dubrulle et al. 1995).  \nIn the above, we assumed that t corr GLYPH<25> GLYPH<10> GLYPH<0> 1 , which is applicable in the ideal limit of MRI-turbulence (Sano et al. 2004; Johansen et al. 2006; Carballido et al. 2011). However, for nonideal e GLYPH<11> ects as ambipolar di GLYPH<11> usion, the correlation time is expected to be longer (Bai & Stone 2011; Zhu et al. 2015). From  \nFor small stopping times ( GLYPH<28> s = 10 GLYPH<0> 2 ; left panels) the SCA and the SEOM overlap. From a numerical perspective the SCA is preferable as it is computationally much less intensive than the SEOM-method. For GLYPH<28> s = 10 GLYPH<0> 1 the distributions slightly di GLYPH<11> er, as can best be seen from the lower panels. For GLYPH<28> s = 1 particles the di GLYPH<11> erences between the two methods amount to several tens of percents at z = 0, while towards the tails of the distribution the relative di GLYPH<11> erence is larger even. The SEOM, however, is in perfect agreement with the Youdin & Lithwick (2007) theory on di GLYPH<11> usive transport. The reason is that, like Youdin & Lithwick (2007), the SEOM accounts for the vertical oscillation (epicyclic motion) of particles, whereas the SCA does not. The SCA does not account for the pebble's inertia and also does not involve a correlation time. From Eq. (15) we deduce that the SCA becomes invalid for GLYPH<28> s > ( t corr GLYPH<10> ) GLYPH<0> 1 .  \nFig. 2. Vertical distribution for particles of di GLYPH<11> erent stopping times: GLYPH<28> s = t stop GLYPH<10> = 10 GLYPH<0> 2 ( left ), 10 GLYPH<0> 1 ( center ) and 1 ( right ). Histograms plot the numerically obtained distributions with the stochastic equation of motion (SEOM; gray and red) and strong coupling approximation (SCA; blue) methods. Long t corr runs are shown with red histograms (the t corr = 10 2 GLYPH<10> GLYPH<0> 1 , GLYPH<28> s = 0 : 1 run is displayed in panel c). Thin curves gives the normal distribution with the scaleheight of Eq. (15) (Youdin & Lithwick 2007). Note the di GLYPH<11> erent scaling among the panels.  \n<!-- image -->  \nFig. 3. Density distribution P ( z ) obtained from a vertically varying diffusivity. The di GLYPH<11> usivity profile in terms of GLYPH<11> z is given by the red curve. The particle stopping time is GLYPH<28> s = 10 GLYPH<0> 2 and is taken independent of z . The distributions obtained from integrating the stochastic equation of motions (Eq. (4); black histogram) and the strong coupling approximation (Eq. (10); gray histograms) are consistent. The dark blue curve gives the normal distribution, for the particle scaleheight evaluated in the midplane (i.e., Eq. (15) with GLYPH<11> = 10 GLYPH<0> 3 ). The light blue curves gives the gas distribution, scaled by a factor 0.1.  \n<!-- image -->  \nEq. (15) it is clear that pebbles will be more strongly stratified for large values of the turbulent correlation time t corr. In Fig. 2b a case with a 10 times longer correlation times (but still the same GLYPH<11> z = 10 GLYPH<0> 2 vertical di GLYPH<11> usivity, implying larger, longer-lived but less vigorous eddies) is presented. Its vertical distribution is much narrower than the canonical t corr = 1 GLYPH<10> GLYPH<0> 1 turbulence. A case with t corr = 10 2 GLYPH<10> GLYPH<0> 1 is also shown in Fig. 2c. Clearly, a degeneracy between turbulent correlation time t corr, di GLYPH<11> usivity ( GLYPH<11> z ), and stopping time ( t stop) is present.", '3.3. Vertical gradient in the diffusivity': 'As an application of a more convoluted model, we consider a vertically dependent di GLYPH<11> usion. In terms of GLYPH<11> z we adopt  \nGLYPH<11> z ( z ) = GLYPH<11> mid + GLYPH<11> surface( z = 2 H gas) 4 1 + ( z = 2 H gas) 4 (17)  \nwhere GLYPH<11> mid is the di GLYPH<11> usivity in the midplane and GLYPH<11> surface is the di GLYPH<11> usivity in the upper regions. Such layered accretion (Gammie 1996) when the turbulence is confined to the upper regions, although our parameterization in Eq. (17) is completely arbitrary. We choose GLYPH<11> surface = 0 : 1 and GLYPH<11> mid = 10 GLYPH<0> 3 . This profile is plotted in Fig. 3 by the red curve. We further choose GLYPH<28> s = 10 GLYPH<0> 2 , such that GLYPH<28> s =GLYPH<11> z > 1 in the midplane (indicating settling) and GLYPH<28> s =GLYPH<11> z < 1 in the upper regions (indicating coupling to the gas).  \nIn Fig. 3 the gray and black histograms show the normalized distribution of z obtained with the SCA and SEOM methods, respectively. The methods give consistent result. Clearly, the high j z j regions are sparsely sampled as the probability to find a particle at these heights is low. However, the fact that pebbles can be stirred to these heights at all may be surprising given the low GLYPH<11> mid. This is illustrated with the dark blue curve in Fig. 3, which plots P ( z ; h mid) where h mid GLYPH<25> p GLYPH<11> mid =GLYPH<28> sh gas GLYPH<25> 0 : 3 h gas. In fact, pebbles at j z j & 2 H gas follow the gas distribution (light blue curve). Altogether P ( z ) can be approximated as the sum of two distributions. The majority of the pebbles ( GLYPH<25> 90%) follow the midplane distributions (given by h mid), but about 10% of the pebbles follow the distribution given by the gas scaleheight.', '3.4. Local replenishment of small grains?': 'The ability of turbulence to stir GLYPH<24> mm-sized pebbles from the midplane to many gas scaleheights may o GLYPH<11> er an explanation for the persistent presence of small, (sub)-micron size particles in the disk surface, as deduced from near-IR observations (e.g., Juhász et al. 2010). From a theoretical perspective, the presence of small particles is problematic as they should quickly coagulate among themselves and then settle to the disk midplane  \n(Nakagawa et al. 1986; Tanaka et al. 2005; Dullemond & Dominik 2005). This implies that the grains are replenished. The most common theory postulates that the replenishment occurs in the disk midplane regions. Here grains are produced by highvelocity collisions among pebble-sized particles, which are subsequently transported (by di GLYPH<11> usion) to the disk surface (Birnstiel et al. 2010). However, this is a rather indirect route to replenish small grains. First, it is doubtful if pebbles in the midplane will fragment; the gas may not be su GLYPH<14> ciently turbulent. Second, it takes grains a time GLYPH<24> 1 =GLYPH<11> z GLYPH<10> to di GLYPH<11> use, which is rather long, again when GLYPH<11> z is small; these small grains may simply collide before reaching the surface (Krijt & Ciesla 2016).  \nAlternatively, in a disk with a suitable di GLYPH<11> usivity profile ( GLYPH<11> z ( z )), it is possible to di GLYPH<11> use a small number of pebbles to the disk surface. There, due to the much stronger turbulent velocity field as compared to the midplane, collisions will undoubtedly be catastrophic. To cement these ideas, a coupled transportcollision / fragmentation model need to be considered.', '4. 3D pebble accretion': "Following Paper I we calculate the accretion e GLYPH<14> ciency ( GLYPH<15> ) by conducting a series of N-body integrations to follow the trajectory of pebbles as they drift from orbits exterior to the planet to orbits interior to it. The pebble accretion e GLYPH<14> ciency GLYPH<15> is then found simply by counting the fraction of particles that settle to the planet. While Paper I investigated the role of the planet's eccentricity, we fix ep = 0 here and instead investigate the role of turbulence and the planet inclination (Sect. 4.5). To this e GLYPH<11> ect we let the pebble experience a stochastic motion in the vertical direction, as outlined in Eq. (4). The initial vertical position z ini is given by the steady-state distribution characterized by the pebble scaleheight hP . The initial radial position is set by  \nr ini = ap + 3 b shear + GLYPH<1> r syn( GLYPH<17>; GLYPH<28> s ; N ) + 0 : 1 r Hill (18)  \nwhere b shear = GLYPH<28> 1 = 3 s r Hill is the impact parameter for pebble accretion in the Hill regime and GLYPH<1> r syn( GLYPH<28> s ; GLYPH<17>; N ) the distance a pebble drifts after N synodical orbits. 6 We choose N = 3. When the pebble radius has drifted to a distance rp GLYPH<0> r Hill we stop the calculation. It is then counted as a miss.  \nThe fact that particles are only kicked in the z -direction, allows us to restrict the computations to a narrow ring. In contrast, when we would have considered the general case (turbulence operating in all dimensions), a much larger computational domain would be required because of the possibility of multiple encounters (similar to the eccentric case in Paper I). This complication is the key reason why we consider only turbulence in the vertical dimension.  \nSimilar to Paper I, we express lengths in terms of the disk radius rp and times in units of GLYPH<10> GLYPH<0> 1 . The key parameters are:  \n- 6 We calculate GLYPH<1> r syn from the equation  \nZ rp +GLYPH<1> r syn rp GLYPH<10> p GLYPH<0> vGLYPH<30> ( r 0 ) = r 0 v r d r 0 = 2 GLYPH<25> N (19)  \nwhere v r and vGLYPH<30> are the radial and azimuthal drift velocities (e.g., Weidenschilling 1977). Approximating the integral to second order in GLYPH<1> r syn, we obtain  \nGLYPH<1> r syn = GLYPH<0> GLYPH<17> + p GLYPH<17> 2 + (2 GLYPH<25> N )4 GLYPH<17>GLYPH<28> A GLYPH<28>;GLYPH<17> A GLYPH<28>;GLYPH<17> rp (20)  \nwhere A GLYPH<28>GLYPH<17> = 3 2 (1 + GLYPH<28> 2 s ) GLYPH<0> GLYPH<17> .  \nArticle number, page 6 of 15  \n- -qp = Mp = M ? the planet-to-stellar mass ratio;\n- -h gas = H gas = r the disk aspect ratio at the location of the planet;\n- -GLYPH<17> , a measure of the radial drift velocity of the pebbles (Eq. (6));\n- -GLYPH<28> s = t stop GLYPH<10> , the dimensionless stopping time;\n- -GLYPH<11> z , a proxy for the gas vertical di GLYPH<11> usivity Dzz ;\n- -t corr, the turbulent correlation time, which together with GLYPH<11> z determines the magnitude of the turbulent rms velocities, Eq. (9).  \nAnalogous to Paper I, we consider two integration methods:  \n- -The SEOM, which integrates the particle velocity (Eq. (4));\n- -The hybrid method, which uses the SCA, but switches to the SEOM when the particle is in the vicinity of the planet. Here we take the criterion to switch to the SEOM to be p ( GLYPH<1> x ) 2 + ( GLYPH<1> y ) 2 < 2 r Hill, where GLYPH<1> x and GLYPH<1> y are the distances between the planet and pebble in the x and y Cartesian coordinates. Note the absence of the vertical position in this criterion. The reason is that for some parameter combinations (small planet mass and high GLYPH<11> z ) the vertical step size can easily become larger than the Hill radius in the SCA method.  \nAs explained in Paper I, the SCA assumes that gas and particles are well coupled. It does not capture e GLYPH<11> ects that take place on timescales GLYPH<1> t less than t stop + t corr. On small GLYPH<1> t the particle moves ballistically, while the SCA model keep exhibiting random walk behavior at all (time)scales. Within the Hill sphere, where the numerical timestep will become small, the particle trajectories are therefore incorrect. The strong fluctuations in velocity space also complicates the numerical integration.  \nDi GLYPH<11> erent from Paper I, we do not account for the ballistic regime. Ballistic encounters are encounters in which the pebble is not captured by the settling mechanism, but where accretion occurs by virtue of the pebble hitting the surface of the target. Computationally, we can easily distinguish between ballistic and settling encounters by assigning an arbitrary small physical size to the planet (while keeping its mass). All accretion then occurs through the settling mechanism.", '4.1. Standard model and analytical fits': 'For the standard model we take qp = 3 GLYPH<2> 10 GLYPH<0> 7 (a 0 : 1 M GLYPH<8> mass planet for a solar-mass star), GLYPH<17> = 10 GLYPH<0> 3 , h gas = 0 : 03, t corr = GLYPH<10> GLYPH<0> 1 and vary GLYPH<11> z and GLYPH<28> s . The choices for GLYPH<17> and h gas approximately correspond to a disk location of 1 au; both values will generally be higher in the outer disk. In Fig. 4 symbols give the mean value of " obtained from our numerical integrations, " = N set = N tot, where out of N tot integrations N set pebbles settled to the planet. Error bars correspond to the Poisson error, p N set = N tot. The pebble accretion e GLYPH<14> ciency can be converted into the pebble growth mass, M P ; grw, defined as  \nM P ; grw = qpM ? " (21)  \n(values labeled on the right y -axis). This is the amount of pebbles needed to e -fold the planet\'s mass. Finally, solid curves gives our fit to the data, which will be discussed in the subsequent sections and summarized in Sect. 4.7. The fitting expression is appropriate only for GLYPH<28> s . 1.  \nClearly, where " is higher, fewer integrations are needed to obtain a minimum signal-to-noise. In our integrations N tot is not  \n⊕  \nFig. 4. Pebble accretion e GLYPH<14> ciency vs dimensionless stopping time for several values of the vertical turbulence strength, parameterized by GLYPH<11> z (colors). Open symbols give GLYPH<15> obtained from directly integrating the stochastic equation of motion, while crosses give the results from the hybrid algorithm. Crosses are o GLYPH<11> set by 20% to the right for clarity. Error bars correspond to the Poisson counting error on the number of hits. For GLYPH<11> z = 0 the triangles give " resulting from the local calculations and the black line gives our fit to the 2D-limit ( GLYPH<15> 2d), which we obtained in Paper I. The dotted lines gives GLYPH<15> 3D accounting just for the density correction Eq. (22). The colored solid lines give " accounting for all 3D / 2D-e GLYPH<11> ects (Eq. (40)) as explained in the main text. The right y -axis converts " into the pebble growth mass - the total amount of pebbles needed to e -fold the mass of the planet - assuming a solarmass star.  \n<!-- image -->  \nfixed, but di GLYPH<11> erent for each run in order to obtain a certain signalto-noise ratio, which is determined by the number of settling encounters. Hence, most of the computational e GLYPH<11> ort is spent in the runs where " is small. Results of the SEOM method are shown by the open circles, while the crosses (slightly o GLYPH<11> set) denote the results from the hybrid method. The two methods give consistent results (see also Paper I). For small GLYPH<28> s the SEOM method becomes computationally ine GLYPH<14> cient as it takes the pebble a long time to drift to the interior disk and the integration becomes very sti GLYPH<11> . The hybrid method removes this latter problem and is the method of choice for small GLYPH<28> s . Remarkably, the hybrid method gives acceptable results up to GLYPH<28> s = 0 : 5.  \nResults for the non-turbulent (2D) limit ( GLYPH<11> z = 0; black and gray symbols) were already discussed in Paper I. The e GLYPH<14> ciency decreases towards increasing particle stopping times GLYPH<28> s = 1 particles are accreted at the lowest e GLYPH<14> ciency - because the faster drift by the higher GLYPH<28> s particles outweighs the larger linear cross section. Consequently, the pebble growth mass M P ; grw is large; many pebbles are needed in order to grow the planet. The steepening of the slope that can be noticed at small GLYPH<28> s is the result of a transition from the shear regime (velocities dominated by the Keplerian shear) at high GLYPH<28> s to the headwind regime (velocities dominated by the gas sub-Keplerian motion, GLYPH<17>v K ) at small GLYPH<28> s .  \nWhen GLYPH<11> z > 0 (colored symbols) pebbles are stirred to higher regions, reducing the local density of pebbles in the midplane. This reduces " with respect to the 2D case. As can be seen in Fig. 4 pebbles are most a GLYPH<11> ected when they are small ( GLYPH<28> s GLYPH<28> 1) and when the turbulence is strong (high GLYPH<11> z ), which is of course natural. Hence, in the 3D case there is a preferred pebble aerodynamic size where " ( GLYPH<28> s ) peaks, which occurs approximately at the point when the pebble accretion impact parameter equals the pebble scaleheight, i.e., at the transition of the 2D and 3D regimes. For heavier particles " decreases because of more rapid radial drift, whereas for small particles " decreases because of a  \nreduced local density. However, when GLYPH<11> z & GLYPH<28> s the pebble scaleheight has reached that of the gas, hP GLYPH<25> h gas, and no further reduction is possible. As a result, the curves eventually converge when GLYPH<28> s becomes very small, as is seen in Fig. 4 in the bottomleft corner.  \nIn Ormel (2017), as well as Paper I, we derived that the pebble accretion e GLYPH<14> ciency in the 3D limit reads  \n" 3D = A 3 qp GLYPH<17> hP f 2 set (22)  \nwhere A 3 is a numerical constant, hP the pebble aspect ratio, and f set - the settling fraction - a modulation factor which becomes less than unity when the settling criteria (slow encounters) is no longer fulfilled. The characteristic velocity beyond which settling encounters disappear is  \nv GLYPH<3> = qp GLYPH<28> s ! 1 = 3 v K (23)  \n(Ormel & Klahr 2010; Paper I). Qualitatively, when the approach velocity 7 GLYPH<1> v GLYPH<28> v GLYPH<3> settling is fully operational ( f set = 1), while for GLYPH<1> v GLYPH<29> v GLYPH<3> settling (and therefore pebble accretion) are no longer viable ( f set = 0). Quantitatively, the settling modulation function is fitted empirically by an exponential function (Ormel &Klahr 2010; Visser & Ormel 2016). In Paper I we adopted  \nf set ; I = exp 2 6 6 6 6 6 4 GLYPH<0> a set GLYPH<1> v v GLYPH<3> ! 2 3 7 7 7 7 7 5 : (24)  \nwith a set = 0 : 5. The subscript \'I\' indicates that this expression is valid for a laminar disk (Paper I). We will refine it below (Sect. 4.6) accounting for a turbulence velocity field.  \n⊕  \nFig. 5. E GLYPH<14> ciency of pebble accretion of a 0 : 01 M GLYPH<8> planet ( qp = 3 GLYPH<2> 10 GLYPH<0> 8 ) at 1 au (left) and at 30 au (right). At 1 au h gas = 0 : 03 and GLYPH<17> = 10 GLYPH<0> 3 , while at 30 au we take h gas = 0 : 07 and GLYPH<17> = 5 GLYPH<2> 10 GLYPH<0> 3 . Right axis gives the planet growth mass. At 30 au turbulence more significantly a GLYPH<11> ects the pebble accretion e GLYPH<14> ciency.  \n<!-- image -->  \nThe reduction of " by f set enters quadratically in 3D because the cross section is two-dimensional. By virtue of the rather high planet mass, f set nevertheless evaluates to unity for most runs in Fig. 4. In Fig. 4 the dotted lines give Eq. (22), where we took f set = 1 and evaluated h peb according to Eq. (15). With A 3 = 0 : 39 this matches the numerical results well for small GLYPH<28> s . At high GLYPH<28> s , Eq. (22) clearly overestimates the pebble accretion e GLYPH<14> ciency; the settling e GLYPH<14> ciency is then given by its planar limit.  \nOnly for GLYPH<11> z = 0 : 1 tend the e GLYPH<14> ciencies to lie below the expression given by Eq. (22). The reason is that now the pebble velocity becomes dominated by turbulent motions since GLYPH<11> 1 = 2 z h gas & GLYPH<17> , which suppresses accretion through settling as encounters become too fast for settling, f set < 1, because of a high turbulent velocity. We present a model to include for turbulence e GLYPH<11> ect in f set in Sect. 4.6.', '4.2. Pebble accretion in the outer disk': 'For planets on circular orbits, the pebble accretion e GLYPH<14> ciency is fully determined by the five dimensionless parameters qp , h gas, GLYPH<17> , GLYPH<28> s , and GLYPH<11> z . In the outer disk, the aspect ratio h gas and (as a consequence) GLYPH<17> are usually higher, resulting in much lower e GLYPH<14> -ciencies. This is illustrated in Fig. 5, which shows " ( GLYPH<28> s ; GLYPH<11> z ) for the inner disk (left) and outer disk (right) for a qp = 3 GLYPH<2> 10 GLYPH<0> 8 planet (0 : 01 M GLYPH<8> for a solar-type star). For the outer disk run we have increased h gas and GLYPH<17> by factors of GLYPH<25> 2 and GLYPH<25> 5, respectively. For a standard passively irradiated disk model this would correspond to an increase by a factor 30 in rp ; e.g., we contrast the situation at 1 au (left) with 30 au (right).  \nE GLYPH<14> ciencies in the outer disk are always lower - in particular, the 3D e GLYPH<14> ciencies. In addition, the transition to the 2D limit occurs at a longer stopping time. Both e GLYPH<11> ects are caused by the larger gas scaleheight. A further consequence of a larger h gas is that turbulence velocities become higher compared to the critical threshold v GLYPH<3> (both velocities are lower in the outer disk, but whereas GLYPH<27> z / cs / h gas v K , v GLYPH<3> / v K ). Consequently, for the 30 au run, settling already fails for GLYPH<11> = 10 GLYPH<0> 2 , which manifests itself by the flattening and decrease of the curves.  \nIt has been suggested that pebble accretion is an e GLYPH<11> ective mechanism to grow planets in the outer disk (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch et al. 2015; Jo-  \nhansen & Lambrechts 2017). Compared to planetesimal accretion the key advantage is that capture radii are large due to the large Hill radii, whereas planetesimals su GLYPH<11> er from scattering a negative feedback to the growth of planets (Kobayashi et al. 2010). Nevertheless, as Fig. 5 illustrates, the e GLYPH<14> ciency of pebble accretion in the outer disk is lower than in the inner disk. For example, growing planets in strongly turbulent ( GLYPH<11> z > 10 GLYPH<0> 2 ) disks may require hundreds, if not thousands, of Earth masses in pebbles - numbers that seem rather large in the light of recent ALMA observations (Ansdell et al. 2017; Miotello et al. 2017). Still, invoking (a combination of) low turbulence (as suggested by HL tau; Pinte et al. 2016), pressure bumps (Pinilla et al. 2012), or massive disks (Bitsch et al. 2018a) there is enough leeway to grow planets through pebble accretion in the outer disks. But from an e GLYPH<14> ciency perspective, it is more conducive to grow them in the inner disk.', '4.3. Efficiency at H 2 O iceline for different stellar mass': "There have been a number of studies that argue for the H2O iceline as a preferential site for the formation of the (first) generation of planetesimals and planetary embryos (Cuzzi & Zahnle 2004; Ros & Johansen 2013; Ida & Guillot 2016; Banzatti et al. 2015; Schoonenberg & Ormel 2017; Dra¸˙zkowska & Alibert 2017). In many of these works, the abundance of ices increases just outside the iceline by condensation of H2O vapor that diffused back over the iceline. In addition the surface density increases through a 'tra GLYPH<14> c jam' e GLYPH<11> ect when evaporating ice boulders liberate much smaller grains. For these reasons it is worthwhile to consider the e GLYPH<14> ciency of pebble accretion at the snowline. However, the snowline locations vary with stellar mass. In disks of lower mass stars the snowline will be much closer in, where the disk aspect ratio is likely to be smaller.  \nIn Figure 6 the pebble accretion e GLYPH<14> ciency is plotted at the H2O iceline, contrasting a solar type star (bottom) with a late Mstar (top). The planet mass is 0 : 1 M GLYPH<8> in both cases. Two e GLYPH<11> ects conspire to render e GLYPH<14> ciencies much higher for icelines around low-mass stars. First, the aspect ratio and GLYPH<17> are lower because the iceline lies further in. Second, a planet(esimal) of the same mass around a lower mass star will have a higher mass ratio ( qp ). The pebble accretion cross section is then larger because of the  \nFig. 6. Pebble accretion e GLYPH<14> ciency for a 0.1 M GLYPH<8> planet at the location of the H2O iceline. Top: Iceline around a 0.1 M ? mass star ( h gas = 0 : 03, GLYPH<17> = 10 GLYPH<0> 3 ; qp = 3 GLYPH<2> 10 GLYPH<0> 6 ). Bottom: iceline of a solar mass star ( h gas = 0 : 05, GLYPH<17> = 3 GLYPH<2> 10 GLYPH<0> 3 , and qp = 3 GLYPH<2> 10 GLYPH<0> 7 ).  \n<!-- image -->  \n⊕  \nreduced Keplerian shear and headwind velocities. Figure 6b illustrates the e GLYPH<11> ect for a 0 : 1 M GLYPH<12> M-star for (the same) 0 : 1 M GLYPH<8> mass planet with h gas = 0 : 03 and GLYPH<17> = 10 GLYPH<0> 3 ; i.e., qp is higher by a factor of 10, h gas lower by a factor of 1.6, and GLYPH<17> lower by a factor 3 compared to the solar-type star. Clearly, M-star pebble accretion e GLYPH<14> ciencies are both higher and are less sensitive to variations in GLYPH<28> s and GLYPH<11> z . Therefore, (late type) M-stars can e GLYPH<14> -ciently convert their pebble-sized building blocks into planetary systems. These finding confirm our earlier analytical estimates on a pebble formation origin of the TRAPPIST-1 system (Ormel et al. 2017).", '4.4. Dichotomy of the solar system': 'It was argued by Morbidelli et al. (2015) that pebble accretion is more e GLYPH<14> cient beyond the H2O iceline because icy pebbles, being less prone to collisional fragmentation, are larger (higher GLYPH<28> s ). Specifically, Morbidelli et al. (2015) adopted GLYPH<28> s = 10 GLYPH<0> 1 : 5 for icy pebbles outside the iceline and GLYPH<28> s = 10 GLYPH<0> 2 : 5 for silicate pebbles (just) interior to the iceline and adopted GLYPH<11> z = 10 GLYPH<0> 3 . The curve corresponding to these parameters is the purple curve in the bottom panel of Fig. 6. The positive slope indeed indicates that larger (icy) pebbles have a higher " than the smaller (silicate) pebbles, a situation that generally applies to the 3D limit (Eq. (22) where the pebble aspect ratio hP decreases with higher GLYPH<28> s . Morbidelli et al. (2015) concluded that embryos just outside  \nthe H2O iceline will outcompete inner embryos 8 , arguing that the great dichotomy of the solar system - small Mars next to big Jupiter - is well explained under a pebble accretion scenario.  \nThere is a corollary to this hypothesis, not (explicitly) addressed by Morbidelli et al. (2015). As Fig. 6 shows the e GLYPH<14> ciencies of the GLYPH<28> s = 10 GLYPH<0> 1 : 5 pebbles ( " GLYPH<25> 0 : 4%) are very small: > 99% of the pebbles drift past the planet to enrich the inner solar system. This means that the formation of Jupiter\'s core is associated with hundreds of Earth masses in pebbles drifting into the terrestrial planet region. Evidently, in the scenario outlined by Morbidelli et al. (2015) most pebbles could not have been accreted by the bodies in the inner solar system but must have ended up in the young Sun. This puts a constraint on the structure of the inner disk, i.e., it had to be transparent to pebble drift. No dense planetesimals belts or long-lived pressure maxima (which could trigger formation of super-Earths; Chatterjee & Tan 2014) could have existed in the inner solar system during Jupiter\'s formation.', '4.5. Inclined planets': 'When planet(esimal)s move on orbits inclined with respect to the disk, the pebble accretion e GLYPH<14> ciency can be significantly reduced (Johansen et al. 2015; Levison et al. 2015). Planets of high inclination ip , such that ip > hP , only interact with pebbles over a fraction GLYPH<24> ip = hP of their orbits. Inclined planets, therefore, have a similar e GLYPH<11> ect on the settling e GLYPH<14> ciencies as turbulence: the higher the inclination, the less it interacts with the pebbles. In addition, inclined planets encounter pebbles at an additional velocity ( GLYPH<24> ip v K ), which suppresses accretion once it exceeds v GLYPH<3> .  \nSince the planet\'s inclination has a similar e GLYPH<11> ect on " as turbulent stirring of pebbles - both reduce the amount of interaction between the two-components - an e GLYPH<11> ective scaleheight may be defined as h e GLYPH<11> \' ip + hP . A more precise estimate is obtained by averaging the pebble density over the phase of the planet\'s orbit  \nThese e GLYPH<11> ects are illustrated in Fig. 7, where " is plotted as function of inclination ip (curves) for the laminar disk ( GLYPH<11> z = 0; all pebbles in the midplane), GLYPH<11> z = 10 GLYPH<0> 4 , and GLYPH<11> z = 10 GLYPH<0> 2 . In general, the higher the inclination, the lower the accretion efficiencies. For GLYPH<11> z = 0 inclination e GLYPH<11> ects already become visible for ip GLYPH<24> 10 GLYPH<0> 3 , whereas for GLYPH<11> z = 10 GLYPH<0> 2 the curves only diverge for ip > 10 GLYPH<0> 2 (inclinations are given in radians). Accretion of large GLYPH<28> s -particles in particular are suppressed because of the increase in the approach velocity and the decrease in f set. The planet moves too fast through the pebble plane.  \n1 p 2 GLYPH<25> h e GLYPH<11> = Z 2 GLYPH<25> 0 N ( ip sin t ; hP ) d t 2 GLYPH<25> (25)  \nwhere N ( z ; hP ) is the normal distribution with standard deviation hP and t the phase (mean anomaly). In the limit of ip GLYPH<28> hP , N (0 ; hP ) = 1 = p 2 GLYPH<25> hP and h e GLYPH<11> is equal to the pebble scaleheight. More generally, the formal solution to Eq. (25) reads  \nh e GLYPH<11> = hp exp h i 2 p = 4 h 2 P i I 0( i 2 p = 4 h 2 P ) (26)  \nwhere I 0( x ) is the modified Bessel function of the first kind. For practical purposes Eq. (26) may be approximated as  \nh e GLYPH<11> GLYPH<25> s h 2 P + GLYPH<25> i 2 p 2 GLYPH<16> 1 GLYPH<0> exp h GLYPH<0> ip = 2 hP iGLYPH<17> (27)  \nFig. 7. Planet inclination dependence on pebble accretion e GLYPH<14> ciency. We plot " for our standard parameters ( qp = 3 GLYPH<2> 10 GLYPH<0> 7 , GLYPH<17> = 10 GLYPH<0> 3 and h gas = 0 : 03) as function of GLYPH<28> s (x-axis), GLYPH<11> z (panels) and planet inclination ip (colors). Curves give our fit (Sect. 4.7). The plotted planet inclination values are ip = 0 (black, top), 10 GLYPH<0> 3 , 2 GLYPH<2> 10 GLYPH<0> 3 , 5 GLYPH<2> 10 GLYPH<0> 3 , 10 GLYPH<0> 2 , 2 GLYPH<2> 10 GLYPH<0> 2 , 5 GLYPH<2> 10 GLYPH<0> 2 and 0 : 1 radians. Many of the low ip points and curves overlap.  \n<!-- image -->  \nIn the limit ip GLYPH<29> hP , h e GLYPH<11> \' p GLYPH<25>= 2 ip GLYPH<25> 1 : 25 ip , implying that our first guess ( h e GLYPH<11> = ip ) is o GLYPH<11> by 25%. The reason is that the inclined planet spends most of its time near its end points, whereas it moves quickly through the midplane regions.', '4.6. Role of turbulence - refinement of f set': "We refine the expression for f set (Eq. (24)) accounting for a, possibly anisotropic, turbulent velocity field. Let GLYPH<1> v be the nonturbulent component of the approach velocity with its component, GLYPH<1> v i , pointing in the i th direction. Concerning the planar direction, we already obtained GLYPH<1> v in Paper I, which we now relabel GLYPH<1> vy :  \nGLYPH<1> vy = max( v cir ; v ecc) (28)  \nwhere  \nv cir v K = GLYPH<17> 1 + a cir qp GLYPH<28> s =GLYPH<17> 3 + a sh( qp GLYPH<28> s ) 1 = 3 (29)  \nis the approach velocity in the circular limit. It combines the headwind (approach velocity dominated by GLYPH<17>v K ) and the shear (approach velocity determined by the Keplerian shear velocity) regimes. In addition  \nv ecc = aeep v K (30)  \nis the eccentric velocity. In the above formulae a cir, ae and a sh are all fit constants (see Table 1).  \nSimilarly, for the vertical approach velocity we have  \nGLYPH<1> v z = aiip v K (31)  \nform the planet's inclination. The order of unity prefactor ai is again obtained numerically.  \nWeassume a tri-axial Gaussian velocity distribution of width GLYPH<27> P , centered on GLYPH<1> v . Explicitly, for component i the approach velocity v i is normally distributed  \nP ( v i ) = 1 GLYPH<27> P ; i p 2 GLYPH<25> exp 2 6 6 6 6 6 4 GLYPH<0> ( v i GLYPH<0> GLYPH<1> v i ) 2 2 GLYPH<27> 2 P ; i 3 7 7 7 7 7 5 (32)  \nArticle number, page 10 of 15  \nwhere GLYPH<27> P ; i is the turbulent rms velocity in direction i . Following Youdin & Lithwick (2007) we take  \nGLYPH<27> P ; i = r GLYPH<11> i GLYPH<10> t corr + GLYPH<28> s GLYPH<24> GLYPH<0> 1 = 2 h v K (33)  \nfor the pebble rms velocity where GLYPH<24> is defined in Eq. (16). From Eq. (24) we can write for the accretion probability  \nf acc = exp 2 6 6 6 6 6 4 GLYPH<0> a set v 2 x + v 2 y + v 2 z v 2 GLYPH<3> 3 7 7 7 7 7 5 : (34)  \nThe new, distribution-averaged f set is then obtain by integration over the velocity distribution:  \nf set = Z f acc P ( v x ) P ( vy ) P ( v z ) d v xd vy d v z = Y i exp 2 6 6 6 6 6 4 GLYPH<0> a set GLYPH<1> v 2 i v 2 GLYPH<3> + a turb GLYPH<27> 2 P ; i 3 7 7 7 7 7 5 v GLYPH<3> q ( v 2 GLYPH<3> + a turb GLYPH<27> 2 P ; i ) : (35)  \nFormally, the integration gives a turb = 2 a set, but we relax this constant in order to obtain the best fit to the simulated data.  \nEquation (35) features the following limits:  \n- 1. GLYPH<27> P ; x ; GLYPH<27> P ;y; GLYPH<27> P ; z GLYPH<28> v GLYPH<3> . Turbulence is unimportant and the laminar form of f set is retrieved, Eq. (24).\n- 2. Isotropic turbulence, GLYPH<27> P ; x = GLYPH<27> P ;y = GLYPH<27> P ; z = GLYPH<27> P . This simplifies Eq. (35) to  \nf set = exp 2 6 6 6 6 4 GLYPH<0> a set GLYPH<1> v 2 v 2 GLYPH<3> + a turb GLYPH<27> 2 P 3 7 7 7 7 5 v 3 GLYPH<3> ( v 2 GLYPH<3> + a turb GLYPH<27> 2 P ) 3 = 2 (36)  \n3. GLYPH<27> P ; x = GLYPH<27> P ;y = 0 with turbulence only operating in the vertical direction, which is (by construction) the case considered in this paper and may be applicable to the vertical shear instability (Stoll et al. 2017). Hence, we have used here  \nf set = exp 2 6 6 6 6 6 4 GLYPH<0> a set 0 B B B B B @ GLYPH<1> v 2 y v 2 GLYPH<3> + GLYPH<1> v 2 z v 2 GLYPH<3> + a turb GLYPH<27> 2 P ; z 1 C C C C C A 3 7 7 7 7 7 5 v GLYPH<3> q v 2 GLYPH<3> + a turb GLYPH<27> 2 P ; z : (37)  \nTable 1. Breakdown of the expressions and parameters involved in the pebble accretion e GLYPH<14> ciency.  \n- 4. GLYPH<27> P ; x ' GLYPH<27> P ;y ' GLYPH<27> P ; z ' GLYPH<27> P GLYPH<29> v GLYPH<3> +GLYPH<1> v , the turbulence-dominant limit. In this case  \nf set = 0 B B B B B @ v GLYPH<3> a 1 = 2 turb GLYPH<27> P 1 C C C C C A 3 : (38)  \nTurbulence reduces the accretion, but not exponentially. In a turbulence-dominated velocity field there is always a fraction of particles with velocities low enough to accrete by settling.", '4.7. Summary of the pebble accretion efficiency fit': 'We provide an executive summary on how to generally obtain " . Quantities involving the recipe and corresponding numerical constants are given in Table 1.  \nThe settling e GLYPH<14> ciencies in the 2D and 3D limits read  \n" 2D ; set = A 2 GLYPH<17> s qp GLYPH<28> s GLYPH<1> v v K f set (39a)  \n" 3D ; set = A 3 qp GLYPH<17> h e GLYPH<11> f 2 set (39b)  \nApart from the f set term, the 3D expression is independent of the pebble-particle relative velocity; a characteristic feature of pebble accretion in the 3D limit (Ormel 2017; Paper I). The e GLYPH<11> ective scaleheight, h e GLYPH<11> , Eq. (27), accounts for the reduced interaction between pebble and planetesimal by either turbulent stirring of pebbles or planetesimal inclination. The relative motion or pebble approach velocity (barring turbulence) GLYPH<1> v is given in Eq. (28)  \nfor the planar motion and Eq. (31) for the vertical motion due to planet inclination.  \nThe other velocity-dependent term is the settling fraction f set, which becomes important when either laminar or turbulent velocities exceed the critical settling velocity v GLYPH<3> = ( qp =GLYPH<28> s ) 1 = 3 v K . In the laminar case ( GLYPH<27> = 0) f set is given by Eq. (24), while in the turbulent case ( GLYPH<27> , 0) it is given by Eq. (35). When v GLYPH<3> significantly exceeds GLYPH<27> + GLYPH<1> v the settling fraction evaluates to unity.  \nFinally, we find that the 3D (this Paper) and 2D (Paper I) can be combined as  \n" set = GLYPH<16> " GLYPH<0> 2 2d ; set + " GLYPH<0> 2 3D ; set GLYPH<17> GLYPH<0> 1 = 2 : (40)  \nwhich ensures a smooth transition. All curves shown in Figs. 4-9 follow this recipe.  \nWhen f set GLYPH<28> 1 (e.g., GLYPH<1> v GLYPH<29> v GLYPH<3> or GLYPH<27> GLYPH<29> v GLYPH<3> ) the ballistic regime (see Paper I) takes over. In Paper I we found that the combined e GLYPH<14> ciency may be given  \n" = f set " set + (1 GLYPH<0> f set) " bal (41)  \nHowever, in the 3D limit, growth through ballistic interactions (gravitational focusing or geometric sweepup) is quite slow as pebble accretion cross section are generally much larger than their ballistic counterparts.', '4.8. Neglected effects': 'We list several neglected e GLYPH<11> ects, which may change " :  \n- 1. Aerodynamical deflection. For planetesimals, we have not accounted for the change in v gas in the vicinity of the planetesimals and correspondingly ignored any aerodynamic deflection (Visser & Ormel 2016). This reduction, however, only becomes important when particles are very small, t stop < R =v hw, and only when the encounters operate in the ballistic regime. Under these conditions, turbulence may in fact overcome the aerodynamic barrier (Homann et al. 2016).\n- 2. Pre-planetary atmospheres. We have not accounted for the e GLYPH<11> ects of the early primordial atmosphere that forms around massive bodies in gaseous disks. As a rule of thumb, this affects the density and flow pattern out to a Bondi radius, Rb = GMp = c 2 s . Planets as massive as the thermal mass ( h 3 M ? ) will have such extensive envelopes that our constant-density assumption will break down. A complex flow structure may prevent small particles from accreting to the planet (Ormel 2013), perhaps after their evaporation (Brouwers et al. 2018; Chambers 2017).\n- 3. Pressure maxima and resonances. Such massive planets also a GLYPH<11> ect the pressure profile of the disk, resulting in a pressure maximumat a distance GLYPH<24> H gas from the planet. For this pebble isolation mass, accretion will terminate (Lambrechts et al. 2014; Bitsch et al. 2018b). Similarly, GLYPH<28> s > 1 pebbles can be stopped at resonant location (Weidenschilling & Davis 1985; Picogna et al. 2018).\n- 4. Gas radial flow. In viscous disks, gas flows radial at a velocity GLYPH<24>GLYPH<0> GLYPH<23>= r and adds to the drift velocity of pebbles. This e GLYPH<11> ect hence starts to dominate radial drift motions for GLYPH<11>GLYPH<23> > GLYPH<28> s . It is straightforward to adjust expressions for " (see, e.g., Ida et al. 2016). Similarly, the drift velocity depends on vertical position, v r ( z ), because GLYPH<28> s increases towards the disk surface. We can correct " accordingly, i.e., by using a vertically averaged v r (Takeuchi & Lin 2002; Kanagawa et al. 2017).  \nTable 2. Collected e GLYPH<14> ciency expressions used in recent studies.  \nNotes. In comparing the expressions we adopt the shear-dominated velocity regime in the 2D limit (massive planets). ( a ) Using their Equations (1), (29) and (31).  \n( c ) Did not obtain the correct leading term in the case when the approach velocity was given by the Keplerian shear.  \n- ( b ) In the 2D limit the expression followed Lambrechts & Johansen (2014). Bitsch et al. (2015) also adopts these expressions.  \n( d ) Their Equation (59).  \n- 5. By design, the simulations of this work only considered turbulence in the vertical direction - a simplification that allowed us to carry out a numerical parameter study in a controlled way. By modeling vertical turbulence we fully account for the di GLYPH<11> usion e GLYPH<11> ect on " - i.e., the reduction of the midplane pebble density. This will be una GLYPH<11> ected by addition of planar turbulence. However, the turbulent velocity e GLYPH<11> ect - i.e., the reduction of " because of too fast encounters - is determined by all velocity components. In particular, situation where GLYPH<27> x ; GLYPH<27>y GLYPH<29> GLYPH<27> z are not covered by our integration. Nevertheless, in Sect. 4.6 we formulated a general recipe for f set in a anisotropic velocity field.', '5. Discussion': 'Table 2 compiles expressions for " that have been used in, or derived from, recent studies. For simplicity, we ignore the f set term and list the numerical prefactor belonging to the leading term for the 3D and the 2D limits. In the 2D limit the shear-dominated velocity regime has been adopted (valid for large planets), while in the 3D limit the expression is independent of GLYPH<1> v . We emphasize that this work gives the correct numerical prefactor, as it is calibrated against numerical simulations. On the other hand, the existing literature expressions often employed scaling arguments, where typically the 3D rate is estimated to be a fraction b set = HP of the 2D rate where b set is the pebble accretion impact parameter. It is therefore quite remarkable that the existing literature prefactors lie so close to our calculated values. Nevertheless, in the 2D-limit the literature expression turn out to be too high, up to 50%. This may still be significant, since a factor of two difference means that planet formation by pebble accretion takes twice as long and requires twice the number of pebbles.  \nLike our study Chambers (2014) considered the cases where planet eccentricity (inclinations) dominate and accounts for the velocity e GLYPH<11> ect of turbulence. However, he only considered the turbulent rms velocity for the relative motion between planet and pebble, which suppresses pebble accretion exponentially ( f set GLYPH<29> 1) once GLYPH<27> GLYPH<29> v GLYPH<3> . This is illustrated in Fig. 8 by the gray curve, where the settling fraction is plotted as function of planet mass. Parameters are chosen such that turbulent rms velocities GLYPH<24> GLYPH<11> 1 = 2 z h gas v K are similar to the laminar headwind velocity, GLYPH<17>v K . When f set is calculated by adding in quadrature the laminar and turbulent rms velocities (as in Eq. (24)) it ensures exponential behavior at low qp . However, accounting for a velocity distribution changes the functional behavior of f set and, for the adopted parameters, we obtain the counter intuitive result that turbulence  \nFig. 8. Settling fraction f set as function of planet mass for GLYPH<11> z = 10 GLYPH<0> 2 , h gas = 0 : 05, GLYPH<17> = 5 GLYPH<2> 10 GLYPH<0> 3 , GLYPH<28> s = 0 : 1 and M ? = 1 M GLYPH<12> . Turbulence is assumed to be isotropic and the correlation time t corr = GLYPH<10> GLYPH<0> 1 . Curves give f set when it is calculated based on the laminar-only velocity (dashed), turbulent rms velocity (gray), and for a distribution of velocities (solid).  \n<!-- image -->  \nmass [  \nM  \n]  \nincreases f set (at low qp ; solid curve). Accounting for the distribution, there are always a few particles slow enough for the settling mechanism to operate. This finding is important especially in the outer disk, where even at f set GLYPH<24> 10 GLYPH<0> 3 settling interactions dominate over ballistic interactions.  \nRecently, Xu et al. (2017) measured pebble accretion rates from laminar and MRI-turbulent hydrodynamical simulations. They expressed their result in terms of a dimensionless quantity k abs, which is the ratio of the mass accretion rate ( ˙ M ) normalized to the Hill accretion rate 3 R 2 H GLYPH<10>GLYPH<6> P . In these units, our expressions for the 2D and 3D e GLYPH<14> ciencies (Eq. (39)) read: 9  \nk abs ; 2D = 2 : 8 GLYPH<28> 1 = 2 s v t 0 B B B B B @ qp GLYPH<17> 3 ! 1 = 3 + a cir GLYPH<28> s qp GLYPH<17> 3 ! 4 = 3 1 C C C C C A GLYPH<0> 1 + a sh GLYPH<28> 1 = 3 s (42a)  \nk abs ; 3D = 3 : 4 q 1 = 3 p GLYPH<28> s hP = 3 : 4 0 B B B B @ qp h 3 gas 1 C C C C A 1 = 3 GLYPH<28> s hP = h gas (42b)  \nwhere we did insert GLYPH<1> v in Eq. (42a), but have for clarity omitted the f set modulation factor.  \nAs a note in passing, Eq. (42a) depends, apart from GLYPH<28> s , only on the quantity qp =GLYPH<17> 3 . The reason is that the pebble equation of motion, expressed in Hill units, only contains this parameter 10 . Likewise, the 3D-rates also depend on a single, but di GLYPH<11> erent, parameter (assuming h e GLYPH<11> can be quantified in terms of the stopping time). In the general case, then, two parameters (or three if GLYPH<28> s is included) are necessary to calculate the Hill accretion rate. In their local shearing box simulations, Xu et al. (2017) choose to fix the thermal mass qp = h 3 gas and the quantity GLYPH<17>= h gas. For the Hill accretion rate ( k abs), the problem is then fully specified. However, to calculate the e GLYPH<14> ciencies ( " ) the degeneracy that existed among qp , GLYPH<17> , and h gas is broken qp , GLYPH<17> , and h gas each need to be specified. This is because the accretion rate is a local quantity, whereas in order to calculate " the radial pebble flux must be known, i.e., the disk circumference should be specified.  \nFig. 9. Hill-normalized accretion rates for a planet of thermal mass qp = h 3 gas = 3 GLYPH<2> 10 GLYPH<0> 3 and GLYPH<17>= h gas = 0 : 1. Data points give k abs obtained from hydrodynamical laminar (2D), ideal MRI (ID), and resistive (ambipolar di GLYPH<11> usion) MRI (AD) simulations reported by Xu et al. (2017). Solid lines give our fits, where we adopt GLYPH<27> from the same simulations. Data kindly provided by Ziyan Xu.  \n<!-- image -->  \nIn Fig. 9 symbols correspond to the runs conducted by Xu et al. (2017, their Figure 4a) for a thermal mass of qp = h 3 gas = 3 GLYPH<2> 10 GLYPH<0> 3 and a headwind parameter of GLYPH<17> = 0 : 1 h gas. Blue symbols correspond to the hydrodynamic runs (non-turbulent), green symbols to the ambipolar di GLYPH<11> usion runs, and red symbols to the ideal MRI runs. As expected, the accretion rate decreases in the turbulent case and more so in the ideal MRI run than in the AD-restive run. Solid curves give the corresponding k abs obtained from Eq. (42), f set, and the 2D and 3D averaging formula (Eq. (40)). 11 Fits are plotted until GLYPH<28> s = 1 beyond which they loose their validity.  \nThe results of our analytical model agree well with Xu et al. (2017). For the turbulent runs, we find that the reduction of the accretion rate is mostly due to turbulent di GLYPH<11> usion (turbulence lofting pebbles away from the midplane) rather than the turbulent velocity e GLYPH<11> ect. This explains why the ideal MRI run features lower accretion rates than the ambipolar di GLYPH<11> usion runs. Towards GLYPH<28> s = 0 : 1-1 we see that the ambipolar di GLYPH<11> usion run and the hydrodynamic run converge, as the accretion becomes 2D and turbulent di GLYPH<11> usion does not enter k abs. However, we also find that the ideal MRI run does not entirely converge on the hydrodynamic run. This can be attributed to the turbulent velocity e GLYPH<11> ect; f set . 1 because GLYPH<27> & v GLYPH<3> for the GLYPH<28> s \' 0 : 1-1 particles in the MRI run. In other words, settling is marginally failing. We anticipate a much stronger reduction for either smaller planets or more vigorous turbulence.', '6. Summary': 'In this work we have developed a general framework for the stochastic equation of motion (SEOM) of pebble-sized particles. The SEOM is described by the particle\'s stopping time, by the gas di GLYPH<11> usivity D gas and correlation time t corr, and by gravitational forces. Using the SEOM we have investigated the vertical transport of particles in disks.  \n- 1. From the SEOM we obtain the strong coupling approximation (SCA) by taking the limit of t stop ! 0, as was used in previous studies (Ciesla 2010; Zsom et al. 2011). For small t stop and t corr the SEOM becomes consistent with the SCA.\n- 2. Exploring the e GLYPH<11> ect of a vertically dependent di GLYPH<11> usivity, GLYPH<11> z ( z ), such that the midplane di GLYPH<11> usivity GLYPH<11> mid GLYPH<28> GLYPH<28> s and the surface di GLYPH<11> usivity GLYPH<11> surface GLYPH<29> GLYPH<28> s , we find a two component distribution for the pebbles. Although most of the pebbles are concentrated in the midplane (following the scaleheight given by GLYPH<11> mid), a small fraction of pebbles are distributed according to the gas scaleheight. Then, high-velocity collisions among these pebbles could explain the prolonged presence of small, micron-sized grains in the disk surface.  \nAs its main application, we have integrated pebble trajectories to find the pebble accretion e GLYPH<14> ciency " .  \n- 3. Compared to the laminar disk, turbulence reduces " in two ways: it diminishes the local density of pebbles at the midplane through di GLYPH<11> usion and it increases the rms velocity of particles. The latter becomes important for low planet masses, where fulfilling the settling condition becomes more di GLYPH<14> cult.\n- 5. Together with the 2D expressions already derived in Paper I, we have formulated a general prescription for " as a function of pebble properties (aerodynamical size), planet properties (mass, inclination, eccentricity), and disk properties (pressure profile, gas density, turbulence properties). This prescription is summarized in Sect. 4.7.\n- 4. Because of the disk geometry, pebble accretion is more efficient in the inner disk. From an e GLYPH<14> ciency perspective, accretion around low-mass stars is also favored because pebble capture radii are larger around low-mass stars.  \nFinally, we remark that " has been defined with respect to a single planet. A small " therefore does not necessarily imply that pebble accretion is globally ine GLYPH<14> cient. Indeed, in the case of a planetesimal belt (or multiple planets) the total filtering efficiency may well reach unity, whereas the individual " GLYPH<28> 1 (Guillot et al. 2014). This raises the question of how pebble accretion proceeds when multiple seeds are present.  \nWhile it is likely to result in a more e GLYPH<14> cient pebble sweepup, a multi-seed scenario could also suppress planet growth because of the mutual dynamical excitation among the protoplanets (Levison et al. 2015). In a following work, we will investigate the e GLYPH<14> cacy of pebble accretion when pebbles interact with a narrow planetesimal belt (Liu et al., in prep). In such cases, the combined planetesimal and pebble coagulation can best be studied with N-body techniques, where the pebble accretion rate on the N bodies is obtained using the prescription for " summarized in Sect. 4.7.  \nAcknowledgements. The authors are supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). We thank Xuening Bai and Ziyan Xu for sharing and discussing their simulations, and Ramon Brasser and Sebastiaan Krijt for proofreading the manuscript. The comments from the referee improved the clarity of the manuscript.', 'References': 'Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240  \nBai, X.-N. & Stone, J. M. 2011, ApJ, 736, 144  \nBai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152  \nBalbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214  \nBanzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15  \nBirnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79  \nBirnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11  \nBitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 Bitsch, B., Lambrechts, M., & Johansen, A. 2018a, A&A, 609, C2 Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018b, ArXiv e-prints:1801.02341 Brouwers, M. G., Vazan, A., & Ormel, C. W. 2018, A&A, 611, A65 Carballido, A., Bai, X.-N., & Cuzzi, J. N. 2011, MNRAS, 415, 93 Chambers, J. 2017, ApJ, 849, 30 Chambers, J. E. 2014, Icarus, 233, 83 Charnoz, S., Fouchet, L., Aleon, J., & Moreira, M. 2011, ApJ, 737, 33 Chatterjee, S. & Tan, J. C. 2014, ApJ, 780, 53 Ciesla, F. J. 2010, ApJ, 723, 514 Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496  \n- Dra¸˙zkowska, J. & Alibert, Y. 2017, A&A, 608, A92  \nCuzzi, J. N. & Zahnle, K. J. 2004, ApJ, 614, 490  \n- Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237\n- Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99 Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2017, ApJ, 835, 230\n- Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971 Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150  \nGammie, C. F. 1996, ApJ, 457, 355  \nGressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291  \n- Gressel, O. 2017, in Journal of Physics Conference Series, Vol. 837, Journal of Physics Conference Series, 012008  \nGressel, O., Nelson, R. P., & Turner, N. J. 2012, MNRAS, 422, 1140  \n- Guilloteau, S., Dutrey, A., Wakelam, V., et al. 2012, A&A, 548, A70\n- Guillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72\n- Homann, H., Guillot, T., Bec, J., et al. 2016, A&A, 589, A129  \nHottovy, S., McDaniel, A., Volpe, G., & Wehr, J. 2015, Communications in  \nMathematical Physics, 336, 1259  \nHughes, A. M., Wilner, D. J., Andrews, S. M., Qi, C., & Hogerheijde, M. R.  \n2011, ApJ, 727, 85  \n- Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292\n- Ida, S. & Guillot, T. 2016, A&A, 596, L3\n- Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72  \nJohansen, A. & Lambrechts, M. 2017, Annual Review of Earth and Planetary  \nJohansen, A., Klahr, H., & Mee, A. J. 2006, MNRAS, 370, L71', 'Sciences, 45, 359': 'Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Ad-  \nvances, 1, 15109 Juhász, A., Bouwman, J., Henning, T., et al. 2010, ApJ, 721, 431 Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142 Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icarus, 209, 836 Kobayashi, H., Tanaka, H., & Okuzumi, S. 2016, ApJ, 817, 105 Krijt, S. & Ciesla, F. J. 2016, ApJ, 822, 111 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 Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322 Liu, B. & Ormel, C. W. 2018, ArXiv e-prints:1803.06149 (Paper I) Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113 Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258,', '418': 'Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610  \nNakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 Nelson, R. P. & Gressel, O. 2010, MNRAS, 409, 639  \nOkuzumi, S. & Ormel, C. W. 2013, ApJ, 771, 43  \nOrmel, C. W. 2017, in Astrophysics and Space Science Library, Vol. 445, Astro-  \nOrmel, C. W. 2013, MNRAS, 428, 3526  \nphysics and Space Science Library, ed. M. Pessah & O. Gressel, 197  \nOrmel, C. W. & Klahr, H. H. 2010, A&A, 520, A43  \nOrmel, C. W. & Cuzzi, J. N. 2007, A&A, 466, 413  \nOrmel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1  \nPaardekooper, S.-J., Rein, H., & Kley, W. 2013, MNRAS, 434, 3018  \nOrmel, C. W. & Okuzumi, S. 2013, ApJ, 771, 44  \nPan, L. & Padoan, P. 2010, Journal of Fluid Mechanics, 661, 73  \nPinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114 Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62  \nPicogna, G., Stoll, M. H. R., & Kley, W. 2018, ArXiv e-prints:1803.08730  \nRein, H. & Papaloizou, J. C. B. 2009, A&A, 497, 595  \nSafronov, V. S. 1969, Evolution of the Protoplanetary Cloud and Formation of Earth and the Planets, ed. V. S. Safronov (Moscow: Nauka. Transl. 1972 NASA Tech. F-677)  \nRos, K. & Johansen, A. 2013, A&A, 552, A137  \nSano, T., Inutsuka, S.-i., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321  \nShakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337  \nSchoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21  \nSmoluchowski, M. V. 1916, Zeitschrift für Physik, 17, 557  \n- Stoll, M. H. R., Kley, W., & Picogna, G. 2017, A&A, 599, L6 Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A,  \nStoll, M. H. R. & Kley, W. 2014, A&A, 572, A77', '596, A74': 'Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49 Uhlenbeck, G. E. & Ornstein, L. S. 1930, Physical Review, 36, 823 van Kampen, N. G. 1992, Stochastic Processes in Physics and Chemistry Visser, R. G. & Ormel, C. W. 2016, A&A, 586, A66 Völk, H. J., Jones, F. C., Morfill, G. E., & Roeser, S. 1980, A&A, 85, 316 Weidenschilling, S. J. 1977, MNRAS, 180, 57 Weidenschilling, S. J. & Davis, D. R. 1985, Icarus, 62, 16 Xu, Z., Bai, X.-N., & Murray-Clay, R. A. 2017, ApJ, 847, 52 Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588 Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81 Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73  \nTakeuchi, T. & Lin, D. N. C. 2002, ApJ, 581, 1344', 'Appendix A: Proof of Equation (10)': 'In deriving Eq. (10), we follow the proof outlined by Hottovy et al. (2015). The first step is to write Eq. (4) as:  \nd x = v d t (A.1a)  \nt stopd v = F d t GLYPH<0> γ v d t + GLYPH<27> d Wt (A.1b)  \nIn this equation, the stochastic variable GLYPH<16> t has been combined with the velocity into one vector v = ( v; GLYPH<16> t ). Similarly, x contains an additional parameter, say GLYPH<29> , but this is entirely dummy. Comparing with Eq. (4) we therefore have:  \n˜ F = v gas + v hs + F g t stop 0 ! (A.2a)  \nγ = 1 GLYPH<0> p Dz = t corr 0 t stop = t corr ! (A.2b)  \nGLYPH<27> = 0 B B B B B @ 0 q 2 t 2 stop = t corr 1 C C C C C A (A.2c)  \nIt must be emphasized that we only treat here the z coordinate. When the full 3D equation of motion is considered, v will be vector of length six and γ a matrix of 36 elements, containing all entries of the di GLYPH<11> usion tensor Dij .  \nIn the limit of t stop ! 0, Hottovy et al. (2015) proves that x can be described by the SDE  \nx = h γ GLYPH<0> 1 ˜ F + S i d t + γ GLYPH<0> 1 GLYPH<27> d W t (A.3)  \nwhere γ GLYPH<0> 1 is the inverse of γ and S is the noise-induced drift term - a vector whose i th component is defined  \nS i = X j ; l @ @ xl h ( γ GLYPH<0> 1 ) i j i J jl (A.4)  \nwith J the solution of the Lyapunov equation  \nγ J + J γ y = GLYPH<27>GLYPH<27> y : (A.5)  \nSince we consider here only the z coordinate, Eq. (A.3) reduces too  \nS 1 = X j @ @ z ( GLYPH<13> GLYPH<0> 1 )1 j J j1 : (A.6)  \nThe equations can now be solved. Inverting γ we find  \nγ GLYPH<0> 1 = 1 p Dzztc = t stop 0 t corr = t stop ! (A.7)  \nFurthermore, solving the Lyapunov equation, we find  \nJ = 0 B B B B B B B B B B B B @ t stop t stop + t corr D ( z ) t stop t stop + t corr p D ( z ) t corr t stop t stop + t corr p D ( z ) t corr t stop 1 C C C C C C C C C C C C A (A.8)  \nwith which the noise-induced drift term becomes  \nS 1( z ) = 1 2 D 0 ( z ) 1 + t stop = t corr (A.9)  \nand we retrieve Eq. (10).', 'Appendix A.1: Stratonovich interpretation': "A feature peculiar to stochastic integrals is that equations of the form  \nZ t +GLYPH<1> t t B ( x )d Wt (A.10)  \nare ill-defined when B is a function of position. In contrast to ODEs, it does matter whether the integrand is evaluated at t (i.e., B ( x ) = B ( x [ t ]) - Ito's choice), at t + GLYPH<1> t , or whether we let  \nZ GLYPH<1> t t B ( x )d Wt GLYPH<25> 1 2 B x [ t ] + x [ t + GLYPH<1> t ] 2 ! ( Wt +GLYPH<1> t GLYPH<0> Wt ) (A.11)  \n(Stratonovich' choice). In contrast to ODEs, these definitions will produce di GLYPH<11> erent results when B is not constant (van Kampen 1992). Specifically, with Stratonovich', rather than Ito's interpretation for the stochastic integral a term  \nGLYPH<0> 1 2 @ @ x B @ B @ x : (A.12)  \nshould be added to the RHS of the Fokker-Planck Equation (12).  \nIn deriving Eq. (10), as well as the conversion from Eq. (10) to Eq. (13) we have followed Ito's interpretation. On the other hand, under Stratonovich interpretation, the RHS of Eq. (10) and the RHS of Eq. (13) gain an additional term GLYPH<0> 1 2 D 0 :  \nd x t = F t stop + v gas + v hs GLYPH<0> t stop D 0 gas 2( t stop + t corr) + p 2 D d Wt : (A.13)  \nHence, the SDE for the strong coupling approximation depends on the interpretation rule (a fact not highlighted by Ciesla 2010 or Zsom et al. 2011). Which form should we choose?  \nMore practically, these di GLYPH<11> erences are expressed in our choice of the numerical integration scheme. When Ito's interpretation is adopted, the stochastic equation must be integrated with a corresponding numerical scheme, of which the Euler method is the simplest example. Similarly, the Stratonovich equation should be integrated with an appropriate numerical scheme, where the midpoint scheme (Heun's method) is the simplest example. In our code, where we use an Runge-Kutta method, which is a generalization of a midpoint scheme, we therefore adopt Eq. (A.13), instead of Eq. (11).  \nThe underlying reason between the Ito and Stratonovich interpretations reflects the nature of the stochastic forcing (van Kampen 1992). Ito's interpretation would hold, for example, when the stochastic forcing amounts to a series of infinitely short 'pulses' with every pulse completely independent. This interpretation is often used in finance However, it may be argued that for physical problems, where the correlation time is never really infinitely small, Stratonovich amounts to the more correct model (van Kampen 1992)."}