How dust fragmentation may be beneficial to planetary growth by pebble accretion

Pebble accretion is an emerging paradigm for the fast growth of planetary cores. Pebble flux and pebble sizes are the key parameters used in the pebble accretion models. We aim to derive the pebble sizes and fluxes from state-of-the-art dust coagulation models, understand their dependence on disk parameters and the fragmentation threshold velocity, and the impact of those on the planetary growth by pebble accretion. We use a one-dimensional dust evolution model including dust growth and fragmentation to calculate realistic pebble sizes and mass flux. We use this information to integrate the growth of planetary embryos placed at various locations in the protoplanetary disk. Pebble flux strongly depends on disk properties, such as its size and turbulence level, as well as on the dust aggregates fragmentation threshold. We find that dust fragmentation may be beneficial to planetary growth in multiple ways. First of all, it prevents the solids from growing to very large sizes, for which the efficiency of pebble accretion drops. What is more, small pebbles are depleted at a slower rate, providing a long-lasting pebble flux. As the full coagulation models are computationally expensive, we provide a simple method of estimating pebble sizes and flux in any protoplanetary disk model without substructure and with any fragmentation threshold velocity.


Introduction
In the classical paradigm of planet formation a significant fraction of solids is very quickly converted into kilometer-sized planetesimals. Some of those planetesimals continue to grow rapidly via the runaway growth (e.g., Wetherill & Stewart 1989;Ida & Makino 1993;Kokubo & Ida 1996). However, this fast stage of growth is soon over as the embryo starts to stir its feeding zone and the accretion transitions to oligarchic growth stage, slowing down with embryo mass (e.g., Kokubo & Ida 1998;. By the end of the oligarchic growth, the embryo absorbs all of the material in its immediate vicinity and reaches the isolation mass, which is as low as Mars mass in the terrestrial planets region, but increases with the distance to the star (Kokubo & Ida 2002). The isolation mass may be high enough to reproduce the cores of giant planets but the core growth timescale in the planetesimal-driven scenario becomes prohibitively long to allow for accretion of gaseous atmosphere outside of Jupiter location (Thommes et al. 2002(Thommes et al. , 2003Levison et al. 2010;Johansen & Bitsch 2019).
These drawbacks of the classical planetesimal-driven model motivated development of the alternative scenario called pebble accretion. In this paradigm, the embryo is accreting centimetersized pebbles rather than planetesimals. Gravity and gas drag act together to enhance the cross section of a planetary embryo for its encounters with pebbles and thus speed up the growth (Ormel & Klahr 2010;Lambrechts & Johansen 2012). Because pebbles are rapidly drifting through the protoplanetary disk, the size of the feeding zone increases, allowing the embryos to grow to larger sizes before the pebble flux is halted by planet-disk interactions .
In the pebble accretion scenario, the outcome of planet formation is determined by the sizes and flux of the pebbles drifting through the protoplanetary disk. Kretke & Levison (2014) used large pebbles with a size distribution between 0.5 m and 5 m (corresponding to the Stokes number of 0.25 and 25) and a constant pebble flux of 0.13 M ⊕ yr −1 . They found that reproducing the giant planets in the Solar system is challenging as pebble accretion tends to convert too many planetesimals into large embryos. Levison et al. (2015a) showed that this difficulty may be mitigated if pebbles form gradually so that the growing planetesimals have time to interact and stir each other. Recently, Lambrechts et al. (2019) showed that a high pebble flux leads to formation of giant planets, while a lower flux leads to formation of super-Earths. However, many authors studying the pebble-driven planet growth for simplicity define a fixed size or Stokes number of pebbles and the value of the pebble flux (or the pebble to gas flux ratio) as arbitrary parameters (typically fixed in orbital distance and/or time), independent of the protoplanetary disk model Liu et al. 2019;Lambrechts et al. 2019;Ogihara & Hori 2020;Brasser & Mojzsis 2020;Wimarsson et al. 2020).  presented an algorithm for calculating the pebble flux and pebble Stokes number, which is often used by the community (Levison et al. 2015a,b;Ida et al. 2016;Bitsch & Johansen 2016;Matsumura et al. 2017;Brügger et al. 2018;Bitsch et al. 2019a;Izidoro et al. 2019). In this model, which is based on a simplified scenario of dust evolution, it is assumed that all dust grains grow until the growth timescale becomes longer than the drift timescale. Since the growth timescale strongly increases with distance, the "pebble formation front" moves outward, leading to subsequent portions of solids decoupling from the gas and maintaining inward pebble flux until the front reaches the outer edge of the disk. This picture is roughly consistent with the outcome of dust coagulation models as long as dust growth can proceed without being hindered by fragmentation. However, the increase of collision speeds with aggregate sizes often prevents dust growth to reach those sizes.
Laboratory experiments show that collisions of dust aggregates become destructive already at speeds much lower than predicted for the protoplanetary disk environment. The exact fragmentation threshold speed remains uncertain and it was shown to strongly depend on the porosity and composition of the aggregates (Wada et al. 2011;Meru et al. 2013;Blum 2018). Until recently, it was believed that the icy aggregates outside of the water ice line are significantly more sticky than the dry aggregates inside of the ice line (Aumatell & Wurm 2014;Gundlach & Blum 2015). However, newer laboratory data do not support this view, demonstrating that the icy aggregates break easily (Gundlach et al. 2018;Steinpilz et al. 2019;Musiolik & Wurm 2019).
Nevertheless, most of the works concerned with planet growth via pebble accretion published to date neglect the effect that fragmentation of pebbles has on their sizes and flux (the exceptions are Chambers 2016, Guilera et al. 2020, and Venturini et al. 2020a. In this paper, for the first time, we study the growth of planetary embryo by pebble accretion in connection with a self-consistent dust evolution model considering full size distribution obtained in detailed dust coagulation simulation. This paper is organized as follows. We present the numerical model used in this work in Sect. 2. We present the resulting planet growth tracks and their dependence on the parameters that are the most important to defining dust evolution outcome in Sect. 3. We present a simple way of predicting a realistic pebble flux and sizes which is valid in a fragmentation dominated as well as in drift dominated protoplanetary disk in Sect. 4. In Sect. 5 we discuss our work and we summarize it in Sect. 6.

Protoplanetary disk model
We model protoplanetary disk surrounding a solar mass star with effective temperature of T eff = 5772 K. In the initial condition, we set the gas surface density as a function of the radial distance r to the often used, self-similar solution to the viscous evolution equations (Lynden-Bell & Pringle 1974) where M disk is the initial disk mass, and r c is the characteristic radius, and the dust surface density to where Z is the global solids-to-gas ratio. In all the models presented in this paper, we consider disks with exactly the same initial mass budget. M disk is set to 0.2 M and Z has the standard value of 0.01, which corresponds to the initial mass in solids of about 650 M ⊕ . We consider two values of the characteristic radius, r c = 30 au in the compact disk model and r c = 300 au in the large disk model. The initial gas surface density profiles are presented in the upper panel of Fig. 1 of the disk models are gravitationally stable, with the Toomre Q > 2, as shown in the bottom panel of Fig. 1.
The gas disk evolves viscously, which means that the characteristic radius increases, while the surface density drops with time. The speed of viscous evolution is determined by the α parameter. We set the default value to α = 10 −4 . In Section 3.3, we show the effects of varying the α parameter. In this work, we do not consider disk dispersal by photoevaporation.
The pebble accretion rate is sensitive to the vertical scale height of gas H gas and thus to the disk temperature structure. In this paper, we take into account both stellar irradiation and viscous heating using the simple prescription proposed by Ida et al. (2016, their equation 7 and 8). Most of the disk is heated by stellar irradiation, leading to a shallow temperature profile T irr ∝ r −3/7 . In the inner part of the disk, viscous heating may change the temperature profile to a steeper function of the radial distance, T vis ∝ r −9/10 . In the viscous heating regime, the temperature depends on the gas mass accretion rate, and thus this effect is particularly important in the compact disk model, as visible in the middle panel of Fig. 1, and in models with α = 10 −3 , where the gas mass flux, defined asṀ gas = 3παΣ gas H 2 gas Ω K , with Ω K being the Keplerian frequency, is higher.

Dust evolution
We investigate the dust evolution on a logarithmic radial grid with 120 grid points spaced between 1 au and 1000 au. We compute dust evolution in the azimuthally and vertically averaged framework using the DustPy code, a Python-based version dust coagulation model similar to the code presented by Birnstiel et al. (2010). At every radial distance to the star, dust mass grid with 9 mass bins per mass decade is built and the Smoluchowski equation is used to solve for dust growth and fragmentation. We assume that the dust grains initially have radius of one micron and internal density of 1.25 g cm −3 . The collisions are driven by the Brownian motion, turbulence prescribed with the α param- eter, radial and azimuthal drift, and dust settling. We consider various values of the fragmentation threshold speed v f , between 1 m/s and ∞, corresponding to no fragmentation at all. In line with the newest laboratory results, we do not employ a transition at the water ice line, but rather we use a single value of v f independent on the radial distance.
Radial drift of particles in every mass bin i is computed based on their local Stokes number (Weidenschilling 1977 where v gas is the local gas velocity (driven by gas viscous evolution) and v η is the maximum velocity of the radial drift driven by the local pressure gradient ∂ r P gas . We calculate the v η as v η = ∂ r P gas where ρ gas = Σ gas /( √ 2πH gas ) is the midplane density of gas. We take into account both the Epstein and the Stokes drag regimes when calculating the Stokes number St i . In this work, we do not include the backreaction from dust to gas.
Further details about the dust evolution code can be found in Birnstiel et al. (2010) and in forthcoming Stammler & Birnstiel (in prep.).

Planet growth by pebble accretion
We consider growth of planetary embryos with starting mass of 0.001 M ⊕ (approximately 0.1 Moon mass, or 7 masses of the largest asteroid, Ceres). We model the growth of embryos at various radial locations and do not consider their migration or planet-disk interaction. We do not subtract the dust accreted by an embryo from the dust evolution code, which means that we assume that each model corresponds to a single planet. We consider various possibilities for the point at which we insert the embryos into the simulation t 0 : either we let the pebble accretion start to grow from the beginning of the simulation, or after 10 5 years and after 10 6 years. We assume that the embryos are on circular orbits with zero inclination.
To calculate the pebble accretion efficiency, we use the method provided by  and  who used 3-body calculations to obtain a general recipe for the pebble accretion efficiency of a single planet 1 . In Fig. 2, we show the example of pebble accretion efficiencies obtained using this recipe for a range of planetary embryo masses and sizes at 5 au and at 30 au in the initial condition of the compact disk model with the turbulence parameter α = 10 −4 . As can be seen in this figure, the recipe takes into account both the 2D regime, where the embryo is large enough to accrete from the complete layer of pebbles, and the 3D regime, where the small embryo only has access to some fraction of the pebble layer. The transition from 3D to 2D regime is marked with the dashed-dotted line. The grey dashed line shows the transition between the headwind and shear regimes. The headwind regime, also called Bondi regime, is valid for small embryos or small pebbles, when the pebbles are accreted only from the embryo proximity and the relative velocity between the pebbles and the embryo is determined by the gas flow. Larger embryos enter the shear regime, also called Hill regime, when the pebbles are accreted from much larger area and their approach velocity is determined by the Keplerian shear. The transition between the headwind and shear regimes does not impact the pebble accretion efficiency in the 3D regime, while in the 2D regime it slightly changes the dependence of ε PA on the Stokes number. This effect becomes more important in the outer regions of the disk, when the transition between headwind and shear regime shifts to higher embryo masses. As explained in Ormel (2017), in the 3D regime the pebble accretion efficiency increases with the Stokes number because dust settling is more efficient and the embryo accretes from a denser midplane layer, while in the 2D regime this effect is canceled by the increasing drift speed of pebbles and the decreasing embryo-pebble interaction time. The dotted line in Fig. 2 shows the pebble accretion onset, that is the minimum embryo mass for which the aerodynamic deflection becomes important (Visser & Ormel 2016). For reference, the typical planetesimal size formed in the streaming instability models is about 100 km (Schäfer et al. 2017;Li et al. 2019a;Klahr & Schreiber 2020). Interestingly, these smallest embryos could only grow by accretion of small pebbles and with an extremely low efficiency. The pebble accretion onset shifts toward even larger embryo masses in the outer parts of the disk. Figure 2 shows that the pebble accretion efficiency is mostly dependent on the embryo mass and not on the pebble size, as long as the Stokes number stays between 10 −3 and 1, in the typical "pebble-size" regime. A relatively large embryo, with significant fraction of an Earth mass, is able to grow quickly by accreting a considerable fraction of the pebble flux. It is worth noting that multiple small embryos could grow at the same time, collectively capturing similar fraction of pebbles as a single large embryo (Bitsch et al. 2019b). The growth rate of the planetary embryo depends both on the pebble accretion efficiency and the pebble flux. In this paper, we compute pebble fluxes consistently with the pebble size resulting from our full dust coagulation model.
Our pebble accretion algorithm works as follows. At the planetary embryo location a e , we calculate the Stokes number of the dust particles from every mass bin, St i . Then we calculate the pebble accretion efficiency ε PA (m e , a e , St i ) for the current embryo mass m e and multiply it by the flux of solids corresponding to each bin of the mass gridṁ p (a e , St i ) obtained in the dust evolution code. In the next step, we sum the contributions from each dust mass bin to obtain the embryo growth ratė The growth stops when the embryo reaches the pebble isolation mass. Following , we calculate the pebble isolation mass as where H gas /r is disk aspect ratio. It is worth noting that Ataiee et al. (2018) and Bitsch et al. (2018b) reported that the pebble isolation mass may be increased in turbulent disks while Zormpas et al. (2020) found that the pebble isolation mass is decreased when a realistic equation of state and radiative cooling are taken into account.

Results
We performed seven different models, varying the fragmentation threshold velocity, disk size, and turbulence level. In each of the models, we considered growth of planetary embryos by pebble accretion at 1 au, 5 au, 10 au, 20 au, 30 au, 40 au, and 50 au. However, in none of the models did an embryo located outside of 20 au reach the pebble isolation mass.

Fragmentation threshold
We first focus on the impact of the fragmentation threshold value on the planetary growth. Figure 3 presents results of four runs, all performed with the compact disk model with the turbulence strength α = 10 −4 and different fragmentation threshold speeds v f = 1 m/s, v f = 2 m/s, v f = 10 m/s, and v f = ∞, corresponding to no fragmentation. The left panel of Fig. 3 presents the growth of planetary embryos at 1 au, 5 au, 10 au, and 20 au. At each location, we tested different embryo introduction times t 0 : at the beginning of the simulation (t 0 = 0), after 10 5 years, and after 10 6 years of dust evolution.
The runs with v f = 10 m/s and no fragmentation produced very similar results except for the inner region of the disk (r ≤ in the runs with fragmentation than in the run without fragmentation. In the run without fragmentation, the embryo located at 1 au does not reach its pebble isolation mass irrespective on the introduction time, while in the run with v f = 10 m/s, the innermost embryo reaches the isolation mass after only tens of thousands of years. The reason for the slow growth of the close-in embryos in the model without fragmentation is that pebbles break through the drift barrier and continue growing quickly over the Stokes number of unity, in the Stokes drag regime. This process has been previously described in the literature, see e.g. Brauer et al. The breakthrough only happens in the inner part of the disk, as the outer part is still dominated by the radial drift barrier. Figure 4 presents the average Stokes numbers and grain sizes as a function of radial distance at t = 2 · 10 5 years of evolution in the four models with different fragmentation threshold. It is worth noting that the flux-averaged dust size (and Stokes number) is typically a bit larger than the mass-averaged value, except for the grains that break through the St = 1 barrier. This is because the radial drift velocity increases up to St = 1 and then decreases again (see Eq. 3). The grain size of ∼100 m reached in the inner region in the run without fragmentation is an effect of the upper limit of the size grid used in the dust coagulation calculation rather than a physical value. The large solids do not drift efficiently, which manifests as the low pebble flux at 1 au in the right panel of Fig. 3. Since the grain growth is limited to the upper boundary of the size grid, those results should be treated as approximate, as the pebble flux would drop even more if the growth could proceed without a limit.
Except for reducing the pebble flux, breaking through the St > 1 barrier has another negative impact on planet growth rate. The efficiency of pebble accretion significantly drops for such large dust aggregates (see Fig. 2), as technically they do not count as pebbles anymore. Thus, even though the total mass of solids passing around the close-in embryo over the disk lifetime is essentially the same, in the run without fragmentation the planet does not reach the pebble isolation mass. This result shows that some level of dust fragmentation may be beneficial for the planetary growth via pebble accretion. Thus, we performed additional runs, lowering the fragmentation speed even more, to v f = 2 m/s and v f = 1 m/s, the latter value being roughly consistent with the conclusions from laboratory work (Güttler et al. 2010). As visible in the right panel of Fig. 3, decreasing the fragmentation speed leads to pebble flux reduction. In the runs with the low fragmentation speed, grains stay at smaller sizes (see Fig. 4) and drift slower, leading to lower but long-lasting mass flux, as depletion of small aggregates takes longer. This effect is particularly important if the planetary embryo takes a long time to form. In the left panel of Fig. 3, one can see that if embryos are introduced at 1 Myr of dust evolution, they grow the most in the run with the lowest fragmentation threshold, v f = 1 m/s. Interestingly, for embryos introduced at t 0 = 10 5 years, the v f = 2 m/s is the most favorable value in terms of planetary growth. This is because, while the flux in the runs with the higher fragmentation threshold is already decreasing at that point, the flux in the v f = 2 m/s run stays high for another couple of hundred thousands years (see the right panel of Fig. 3). At the same time, the flux is factor of two higher in the v f = 2 m/s run than in the v f = 1 m/s case.

Disk size
In the results presented in Sect. 3.1, the pebble flux lasts for at most 3 Myrs. However, millimeter and centimeter-sized pebbles are observed in protoplanetary disks older than several Myrs (Wilner et al. 2005;Ricci et al. 2010).  Appendix A) noticed this problem and speculated that limited dust growth in the outer part of the disk may be the solution to retaining pebbles for long timescales. They have also shown that large initial disk size facilitates long-lasting pebble flux, a solution often used in pebble accretion models Sato et al. 2016;Bitsch et al. 2018a). While large disks are bright and thus relatively easy to image at high angular resolution (Andrews et al. 2018), there is evidence that most of the disks may, in fact, be small (Ansdell et al. 2018;Long et al. 2019;Maury et al. 2019;Tobin et al. 2020;Trapman et al. 2020). In this section, we study the impact of disk size on the pebble flux and resulting planetary growth. Figure 5 presents results of two models where the same initial mass was distributed either in a compact disk with the critical radius r c = 30 au or in a much larger disk with r c = 300 au. The other simulation parameters are the same: the fragmentation speed is v f = 10 m/s and the turbulence strength is α = 10 −4 .
The compact disk produces a very high pebble flux, well over 10 −3 M ⊕ /yr, which, however, only lasts for about 3 · 10 5 years and declines over the next couple of hundred thousands years (see the right panel of Fig. 5). With such a high pebble flux, the embryos introduced at the beginning of the simulation take less than one hundred thousand years to reach their isolation mass. The embryo placed at 1 au reaches its isolation mass in little over ten thousand years, which is the fastest growth we observe in any of our runs.
Pebble flux obtained in the large disk model is almost one order of magnitude lower than in the compact disk model, but lasts proportionally longer. Thus, even embryos introduced late during the simulation, at t 0 = 10 6 years, are able to obtain their isolation mass, at least inside of 10 au. Thus, extending the initial size of the disk gives similar results to lowering the fragmentation speed, as we discussed in the previous section. There is however, one important difference, which manifests itself in the planet growth results at and outside of 10 au. In the compact disk, a significant fraction of the mass is initially located inside of 10 au. The mass reservoir available for wide-orbit planets is significantly higher in the large disk than in the compact disk. Thus, for example the embryo introduced at t 0 = 10 5 years at 10 au, reaches its isolation mass in the large disk model but not in the compact disk model. With the large mass reservoir available even at wide orbital distances, it may be surprising that no embryo reaches its isolation mass outside of 10 au in the large disk model. This is because the pebble accretion efficiency for a constant embryo size significantly decreases with radial distance (see Fig. 2). In this paper, we consider the initial embryo mass of 0.001 M ⊕ irrespective of its location (and starting time). More massive embryos could probably still reach pebble isolation mass in the outer parts of the disk.

Turbulence
The turbulence level in protoplanetary disk is a key unknown. Observational constraints on the disk dispersal timescales suggest turbulence levels of α ≈ 10 −2 if the dispersal would be driven solely by accretion driven by turbulent transport of angular momentum (Hernández et al. 2007). However, more recent disk models suggest that the turbulence level may vary within the disk and that large regions of the disk, particularly near to the midplane, may have low levels of turbulence (Lesur et al. 2014;Bai 2016). This would be consistent with the disk images obtained in millimeter wavelengths, where the pebbles appear to be well-settled (Pinte et al. 2016;Villenave et al. 2020), and also consistent with the low upper limits set by measurements of turbulent line broadening (Flaherty et al. 2015(Flaherty et al. , 2020Teague et al. 2016Teague et al. , 2018. Given the uncertainties and difficulties in constraining the turbulence level observationally, in this section, we consider the turbulence strength α as a free parameter. Figure 6 shows results of three runs performed on the backdrop of the compact disk model with different levels of turbulence. We keep the fragmentation speed at v f = 1 m/s in all the runs. As visible in the left panel of Fig. 6, low turbulence promotes the planetary growth irrespective on the distance to the star and the starting time of the embryo. This has multiple reasons. First of all, lower turbulence level leads to better settling of the pebbles and higher pebble accretion efficiency (Ormel 2017). Lower turbulence speed promotes growth to larger pebble sizes, which increases the pebble flux, although this effect is stopped when the fragmentation starts to be dominated by the differential radial drift rather than turbulence (this happens for α 10 −4 , see the right panels of Fig. 6).
It is worth noting that in case of the highest turbulence parameter value that we considered, α = 10 −3 , no embryo reaches the pebble isolation mass. If the turbulence is strong, with α ≥ 10 −3 , the fragmentation threshold velocity needs to be much higher than v f = 1 m/s, which we considered in models presented in Fig. 6, to allow planetary embryos to benefit from pebble accretion.

Simple prediction of the pebble flux
As we showed above, planetary growth by pebble accretion is very sensitive to disk parameters and details of dust evolution. The pebble flux depends on the initial mass budget, disk size, turbulence level, and fragmentation threshold speed. The pebble accretion efficiency is sensitive to pebble size and settling. The DustPy simulations described above are relatively expensive, typically taking several days of computations. In this paper, we propose a simple and efficient method that predicts the size and flux of pebbles for arbitrary protoplanetary disk without substructure, the pebble predictor 2 .
The general mechanics of pebble predictor is closely related to two-pop-py presented by Birnstiel et al. (2012), albeit it a full time integration of the disk is not performed. The pebble predictor uses the initial state of the gas and dust disk surface density Σ gas,0 and Σ dust,0 , temperature T , the turbulence strength α, and the fragmentation threshold velocity v f to predict the flux-2 pebble predictor is publicly available at https://github.com/astrojoanna/pebble-predictor. Version 1.0 of the script used in this paper is permanently available at https://doi.org/10.5281/zenodo.4383153. averaged Stokes number and total flux of pebbles at every location within the disk and at any time point. It is assumed that all the dust grains are initially in the form of micron-sized monomers and they grow in collisions driven by turbulence, which leads to the growth timescale of where Ω K is the Keplerian frequency and r is the radial distance to the central star. The growth proceeds as an exponential function of time t where St 0 is the Stokes number corresponding to the micronsized monomers, until one of the growth barriers is encountered: fragmentation or radial drift. Fragmentation may be driven by turbulence or the differential radial drift. The maximum Stokes number of dust aggregates possible to achieve with respect to the turbulent fragmentation is where f f is a parameter, v f is the fragmentation threshold velocity, α is turbulence strength parameter, and c s is the sound speed. If the radial drift speeds dominate over the turbulent speeds, the maximum Stokes number is where v η is the maximum drift speed (see Eq. 4). We adopt f f = 0.37. We assume that the growth needs to be f d/g = 30 times faster than drift in order for the grains to be not impacted by the radial drift (see Okuzumi et al. 2012, their Sect. 4). We calculate the size of grains that are impacted by radial drift by checking the condition τ drift = f d/g · τ growth , which leads to where η = v η /v K is the pressure gradient parameter. The resulting, representative Stokes number St in every radial bin is then chosen as a minimum of St ini , St f , St df , and St drift , however, it is not allowed to fall below its initial value St 0 .
We use the first part of Eq.
(3) to calculate the radial drift speed v r corresponding to the representative Stokes number and we can calculate the pebble flux aṡ where Σ dust would be the surface density of dust at a given time. This equation is however only correct if the grains grow faster than they drift. In such a case, there are always "enough" large grains, as the drift is depleting them at a slower rate than the growth could replenish them. In the opposite case, the pebble flux needs to be limited by taking into account that the supply of large grains can only be replenished at the growth timescale τ growth . Thus we limit the inward drift velocity v r to v r = min 2v η St which gives us a good estimate of the flux both in the fragmentation dominated and drift dominated regions of the disk. As mentioned above, pebble predictor does not perform a full time integration of the disk. We neglect the evolution of gas disk and we approximate the dust surface density evolution by keeping track of the mass budget. Using the initial Σ dust,0 profile supplied as an input, we calculate the initial mass M out,0 outside of every radial grid point. Using a time grid supplied as an input, pebble predictor calculates how much mass is left outside of every radial grid point at every time point i using the values obtained in the previous time step: The dust surface density used in Eqs (11-12) is then approximated as Σ dust,i = Σ dust,0 · M out,i /M out,0 . Since the time steps supplied in the input time grid are not necessarily fulfilling the CFL condition, this routine cannot be treated as a rigorous integration. Its goal is correcting the pebble flux for the remaining mass budget rather than recovering the actual surface density evolution. Figure 7 shows the comparison of the pebble predictor and DustPy results for various setups varying the disk size, fragmentation speed, and turbulence level. The representative Stokes number from pebble predictor is benchmarked against the flux averaged Stokes number from the full coagulation calculation. On the top right panel of Fig. 7, we additionally show the flux and Stokes number prediction from the analytical model proposed by . In this model, it is assumed that the pebble flux is set by the initial dust surface density at the pebble formation front, which spreads radially with time. The general outcome of this model fits our results well, although the growth timescale is slightly underestimated and the pebble flux starts later than in the numerical results. The most significant difference reveals at later times, as the pebble flux in  model abruptly stops when the pebble formation front reaches the outer edge of the disk. In reality, however, disk viscous spreading and dust diffusion modify the shape of the outer disk edge before the grains there reach pebble sizes. What is more, not all the pebbles at a given location grow and start drifting at once, which extending the duration of the pebble flux as well. Finding analytical prediction that would encompass all the relevant processes is not feasible. Thus, the pebble predictor is a semianalytic model relying on some empirical scaling.
The main limitations of pebble predictor, causing some discrepancies between its predictions and the full simulation results can be summarized as follows.
pebble predictor assumes that the growth is driven by turbulence (only one of the turbulent regimes calculated by Ormel & Cuzzi 2007, to be precise). This assumption breaks in the case of low turbulence, when the collisions are mostly driven by differential drift. For very small particles, the turbulent speeds fall into a different regime (see also the discussion in Powell et al. 2019, their appendix B), which leads to the need to modifying the growth timescale prescription with the empirical dependencies on turbulent strength α and radial distance in Eq. (7). -Only the Epstein drag regime is considered, which leads to error in estimating both the size and pebble flux where the growth enters the Stokes regime, so in the inner parts of disks in runs where dust fragmentation is not considered (see Sec. 3.1). -Gas disk evolution is not included, leading to inaccuracies particularly at the outer regions of the disk, where disk spreading slows down the flow of pebbles. -Similarly, dust advection with the gas flow is not considered in pebble predictor, which may lead to issues with estimating both the flux and pebble size correctly in runs when a significant fraction of solids is coupled to the gas, for example when the turbulence level is high and the fragmentation threshold is low. -Diffusion is not considered, which is particularly important at the outer edge of the disk, where the dust-to-gas ratio gradient is strong (Birnstiel & Andrews 2014). This leads to the characteristic peak in the dust flux visible in the runs involv-ing compact disk and high fragmentation threshold, a feature which is not reproduced by the pebble predictor. -Because the dust advection with gas flow and diffusion are not included, pebble predictor should not be applied to disks with substructure (see Sect. 5.5). A full model, such as DustPy or a more advanced simplified method such as two-pop-py developed by Birnstiel et al. (2012) may be better choice in such case.
Considering all the limitations, the pebble predictor does surprisingly well predicting the dependencies of pebble flux and size on the fragmentation speed, disk size, and turbulence level, as demonstrated in Fig. 7.

Limitations of the (full) model
In this paper, we coupled dust evolution calculated with the DustPy code and embryo growth by pebble accretion. Our models have certain limitations, some of them may impact the planetary growth results discussed in this paper. These are listed below.
-We do not consider that dust grains are composed of various materials and that evaporation of volatile components could cause decrease of the pebble flux in the inner regions of the disk.
-We assume that all the dust aggregates are compact, with constant internal density, and neglect the effects of porosity. -The structure of gas disk is independent of the evolution of solids. This is valid as long as the solids-to-gas ratio stays low, which is true in almost all of the models except for the one without fragmentation, when the grains break through the St = 1 barrier and pileup in the inner part of the disk. -Our models correspond to one-planet-per-disk approach as we do not decrease the pebble reservoir by the amount accreted by each planet. Similarly, we do not include perturbations to disk structure caused by the growing planets. Such multiplanet effects could be important as many embryos would grow from the same pebble reservoir and the outermost embryo would block the pebble flux when it reaches its pebble isolation mass. -In this paper, we assume that the disk is fully turbulent, with a constant α value. Modern disk models predict that large parts of the disk may be free from turbulence and the accretion is primarily driven by magnetized disk winds (Lesur et al. 2014;Bai 2017). In such disks, the viscous heating would not be relevant, leading to lower temperatures in the inner part of the disk than we assume in our models (Mori et al. 2019). Lower temperatures typically promote pebble accretion as the grains grow to larger sizes and are more settled.
On the other hand, the pebble isolation mass is lower for colder disks. -Similarly, we do not include the gas disk dispersal by effects other than viscous accretion. Magnetic and photoevaporative winds could completely change the disk structure, particularly in the later phases of its evolution, which would impact the pebble fluxes and probably abruptly stop planetary growth by pebble accretion.

Planet migration
Perhaps the most significant limitation of our models is that we do not include planetary migration. Interaction between a sufficiently massive planetary core and the gas disk leads to radial migration (Ward 1997). A typical timescale of the inward migration for Earth-mass planet at 1 au is 10 5 years (Ogihara et al. 2015), comparable to the timescale of planetary growth by pebble accretion. Many authors have studied the interplay between the planetary migration and planet growth in pebble accretion scenario. Bitsch et al. (2015) showed that planets mostly migrate in the type I regime while accreting pebbles, but the migration direction may be both inward and outward, depending on their radial location. Planets that reach masses allowing for gas accretion typically outgrow the outward migration region and significantly change their radial location, by tens of au. Nevertheless, planet growth by pebble accretion is generally fast enough for planets to reach the pebble isolation mass before falling onto the central star due to migration (Bitsch & Johansen 2016;Ndugu et al. 2018;.
A connection between planet migration and pebble accretion including the self-consistent dust evolution model will be a subject to our future work.

Validity of the single dust size approach in pebble accretion models
In the models presented in this paper, we used dust evolution model that includes distribution of dust sizes. For every size bin, we use the solid flux calculated by the DustPy code, and multiply it by the pebble efficiency factor ε PA , corresponding to this pebble size, to get the pebble accretion rate. We sum up contribution from each size bin to calculate planet accretion rate, see Eq. (5) 3 . However, since this type of model is computationally expensive, many authors employ simplified models of dust evolution. For example, the two population model for dust evolution proposed by Birnstiel et al. (2012) is recently commonly used to compute dust evolution in protoplanetary disk as it allows us recovering the dust density evolution at a relatively low computational cost, without solving the full coagulation equation (see. e.g., Drążkowska et al. 2016;Cridland et al. 2017;Tamfal et al. 2018;Charnoz et al. 2019;Gárate et al. 2020).
In this paper, we propose even simpler way of predicting pebble size and flux based on one representative grain size at every distance, the pebble predictor (see Sect. 4). It is justified to ask whether these simple models can be reliably used to calculate planetary growth via pebble accretion even though they do not calculate the full size distribution. We tested this by using the results of pebble predictor to calculate embryo growth and compared the results to the results obtained in the DustPy models in Fig. 8. Although there are some differences, the results in terms of the time of reaching the pebble isolation mass or the final embryo mass usually agree within one order of magnitude. The biggest differences are for planets introduced late, at t 0 = 10 6 years. This is because the pebble predictor does not include the gas disk evolution.

Early planet formation
Studies of mass reservoirs of planet forming disks infer that planet(-esimal) formation should start early (Greaves & Rice 2011;Najita & Kenyon 2014;Manara et al. 2018;Tychoniec et al. 2020). There is growing evidence that planet formation process is indeed well underway in young disks (Harsono et al. 2018;Segura-Cox et al. 2020). Pebble accretion may certainly be an avenue to fast planet formation, as young, massive disks should facilitate very high pebble fluxes.
Our results show that the pebble flux is the highest just after the dust aggregates at a given orbital distance reach their maximum size. Thus, if a large enough embryo forms early enough, it may benefits from this high pebble flux and reach its pebble isolation mass already during early stages of disk evolution (Tanaka & Tsukamoto 2019). The pathway to such an early formation of massive planetesimals is not yet understood. The fastest route to planetesimal formation, the streaming instability, needs an enhanced solids-to-gas ratio to operate (Johansen et al. 2009;Bai & Stone 2010). Enhancing their density requires significant redistribution of solids, which may take a long time (Drążkowska et al. 2016). The cold finger effect at water snow line could lead to formation of some planetesimals already during the disk buildup stage (Drążkowska & Dullemond 2018). However, the planetesimals formed via the streaming instability are most likely not large enough to accrete pebbles and intermediate stage of planetesimal accretion is necessary (Liu et al. 2019).

Substructures
In this paper, we focused on the smooth disk model. However, there is observational evidence that substructures are ubiquitous in protoplanetary disks (ALMA Partnership et al. 2015;Long et al. 2018;Andrews et al. 2018). Substructures are thought to modify the evolution of dust, allowing for solids retention over long timescales, consistent with observations (Pinilla et al. 2012b;Li et al. 2019b). They may be a preferential location of planetesimal and planet formation (Haghighipour & Boss 2003;Brauer et al. 2008b). In particular, planetesimal formation may be triggered at the outer edge of a gap carved by a preexisting planet Eriksson et al. 2020;Shibaike & Alibert 2020).
In terms of pebble flux, substructures cause pileup of pebbles, likely speeding up local growth of planet located inside of the pressure bump (Guilera et al. 2020;Morbidelli 2020). Existence of a long lasting substructure may be a way to enable growth of giant planets at large orbital distance. At the same time, however, overall pebble flux through the disk inward from a pressure bump would decrease as some, potentially significant, fraction of solids is held in the substructure.
Comprehensive knowledge of the substructure and detailed dust models are necessary to calculate the resulting pebble flux. One dimensional models of dust evolution in a pressure bump showed that dust growth can proceed to larger sizes inside of a pressure bump (Brauer et al. 2008b;Pinilla et al. 2012a). Dust trapping is hindered by diffusion and advection of small grains with the gas flow. There is a critical size of dust, which depends on the details of the bump structure, below which the grains cannot be trapped (Pinilla et al. 2016;Weber et al. 2018). This effects can, in principle, be modeled with DustPy, which includes both dust diffusion and advection with the gas flow, however the results must be treated with caution, as recent two dimensional models show that the planet induced spiral wakes may modify the growth pattern . The simple flux prediction we proposed in Sect. 4 does not include diffusion and grains advection with the gas flow and thus is not suitable for disks with substructures. We leave development of tools that would be useful in such disks for future work.

Summary
In this paper, we studied planetary growth by pebble accretion in a protoplanetary disk without substructure. We investigated how the most uncertain parameters, the fragmentation velocity, disk size, and turbulence strength, impact the growth rate of embryos placed at different radial locations and times during the dust evolution. We found that planetary growth strongly depends on dust evolution and the pebble sizes and fluxes should be calculated self-consistently with the disk structure. An important byprod-uct of this paper is the pebble predictor, publicly available script that predicts the Stokes number and flux of pebbles in any smooth protoplanetary disk model for a given fragmentation threshold velocity at a very low computational cost (see Sect. 4).
One of our conclusions is that the planet growth by pebble accretion may be extremely fast if a large embryo is formed in the disk while the pebble flux is still high. Compact disk may produce very high pebble flux, leading to the close-in embryos reaching their pebble isolation mass within 10 4 -10 5 years, depending on the fragmentation speed, when the turbulence is not too strong (α ≤ 10 −4 ). High turbulence levels universally reduce the planetary growth rate, while the effects of fragmentation threshold and disk size are not as obvious.
Dust fragmentation keeps dust aggregates at pebble sizes, for which pebble accretion is the most efficient. Without fragmentation, dust growth would proceed over the pebble size limit, making the solids useless for growing planets by pebble accretion. Strong fragmentation keeps the pebbles small and decreases their drift speed, lowering pebble flux and slowing down the growth of embryos. At the same time, the low fragmentation threshold promotes keeping significant fraction of solids in the protoplanetary disk over long time and enables growth of embryos that formed late.
Compact disks are favorable for very rapid growth of closein planets. In contrast, large disks are needed for formation of giant planets on wide orbits. Even with the large disk model, growing cores of giant planets outside of 10 au is challenging as the efficiency of pebble accretion drops with radial distance. One possible solution to this problem may be existence of dust traps that could facilitate planet growth by pebble accretion.