# Challenges in Forming the Solar System's Giant Planet Cores via Pebble Accretion

## Abstract

Though mass rocky/icy cores are commonly held as a prerequisite for the formation of gas giants, theoretical models still struggle to explain how these embryos can form within the lifetimes of gaseous circumstellar disks. In recent years, aerodynamic-aided accretion of “pebbles,” objects ranging from centimeters to meters in size, has been suggested as a potential solution to this long-standing problem. While pebble accretion has been demonstrated to be extremely effective in local simulations that look at the detailed behavior of these pebbles in the vicinity of a single planetary embryo, to date there have been no global simulations demonstrating the effectiveness of pebble accretion in a more complicated, multi-planet environment. Therefore, we have incorporated the aerodynamic-aided accretion physics into LIPAD, a Lagrangian code that can follow the collisional / accretional / dynamical evolution of a protoplanetary system, to investigate how pebble accretion manifests itself in the larger planet formation picture. We find that under generic circumstances, pebble accretion naturally leads to an “oligarchic” type of growth in which a large number of planetesimals grow to similar-sized planets. In particular, our simulations tend to form hundreds of Mars- and Earth-mass objects between 4 and 10 AU. While the merging of some oligarchs may allow them to grow massive enough to form giant planet cores, leftover oligarchs lead to planetary systems that cannot be consistent with our own solar system. We investigate various ideas presented in the literature (including evaporation fronts and planet traps) and find that none easily overcome this tendency toward oligarchic growth.

###### Subject headings:

Planets and Satellites: Formation## 1. Introduction

The core accretion model (Pollack et al., 1996) has been the most successful model for giant planet formation for almost 20 yr now, and yet it still faces significant challenges. According to this model, giant planets emerge in two steps: first a massive rocky/icy core forms and only then can a hydrogen/helium envelope be accreted from the surrounding gaseous protoplanetary disk. This solid core must form quickly, within the few-million-year lifetime of the gas disk (Haisch et al., 2001; Hernández et al., 2008; Williams, 2012). The exact required mass of this core is still under active debate, but values around 10 Earth masses () are generally invoked (although it is possible that masses as low as 5 may be sufficient; Hori & Ikoma (2011)).

Unfortunately, it is quite challenging to build 10 cores quickly using standard models (Ngo et al., 2012). The crux of the issue is that, as a proto-core grows, it excites the eccentricities and inclinations of surrounding planetesimals. This dynamical excitement increases the relative velocity of close encounters, decreasing the effectiveness of gravitational focusing and thus the likelihood of actual collisions. The result is core formation timescales that are too long. In order to get around this problem, some authors have suggested that, if planetesimals are small, then the eccentricities will be effectively damped as a result of either mutual collisions or gas drag (e.g. Goldreich et al., 2004). This damping yields dynamically cold planetesimals that are strongly affected by gravitational focusing and hence should have large capture cross sections (e.g. Rafikov, 2004, ; Figure 1). Unfortunately, if planetesimals lose sufficient energy to stay dynamically cold, then this is the same condition for the planets to open gaps in the planetesimal disk (Levison et al., 2010). Accretion will stop. Therefore we are left with a conundrum; neither large nor small planetesimals seem to be able to effectively form giant planet cores.

However, if one looks at even smaller objects, pebbles or boulders with stopping times () comparable to the orbital period, a possible solution emerges. These objects (which we call pebbles as a shorthand) are very effectively damped by gas drag so that they will be dynamically quite cold. This drag also causes these pebbles to drift extremely rapidly in a normal pressure-supported gaseous disk (Adachi et al., 1976; Weidenschilling, 1977a). As the pebble’s drift timescale is much shorter than nearby a planet’s synodic period, the planet cannot effectively open gaps in the pebble disk. Now, if an object is too small (its stopping time is short compared to the time for it to encounter the planetesimal) then it will be so well coupled to the gas that it follows the gas and flows around the planetesimal and will not be efficiently accreted. However, if the stopping time is comparable to the encounter time then the particle can be deflected out of the gas stream line and thereby lose enough energy during the encounter to be gravitationally captured by the planetesimal. This means that as these pebbles drift by a planetesimal they will spiral inward and can be accreted very efficiently, with accretion cross sections as large as the entire Hill sphere (Johansen & Lacerda, 2010; Ormel & Klahr, 2010, hereafter OK10).

Furthermore, there may be good reason to believe that the early protoplanetary disk was in fact filled with these types of pebbles. While binary collisions appear effective to grow dust grains from micron to pebble sizes, once the solids reach pebble size their growth is stifled for three reasons. First, as particles grow larger their sticking efficiency is reduced, possibly leading to a regime in which collisions are more likely to lead to “bouncing” rather than accretion (Zsom et al., 2010). Additionally, as particles grow larger they can attain larger relative velocities, increasing the likelihood that collisions will be destructive (Blum & Wurm, 2000, 2008). Furthermore, if the do grow then, as discussed above, they will drift inward very rapidly.

However, given enough of these pebble-sized objects, one can expect interactions between the solids and the gas (such as the streaming instability Youdin & Goodman, 2005; Youdin & Johansen, 2007) to concentrate pebbles into clumps that can become gravitationally unstable, directly forming large planetesimals (Johansen et al., 2007). If one looks at the present-day asteroid belt and Kuiper Belt there is evidence that the planetesimals were in fact originally somewhere around 100 km in size. (Morbidelli et al., 2009; Nesvorný et al., 2010). Therefore, it seems plausible that the early solar system consisted only of large ( km-sized) planetesimals and a sea of small pebbles.

Indeed, if a substantial fraction of the disk’s solids are in pebbles, then as the pebbles rapidly drift past planetesimals with very large capture cross sections, the accretion rate can be extremely high. Lambrechts & Johansen (2012, hereafter LJ12) coined the term “pebble accretion” to describe this interesting process in which pebbles are swept up by existing planetesimals and postulated that this very rapid accretion could lead to the formation of the cores of giant planets.

In this first of a series of papers we investigate the feasibility of forming the giant planets in our solar system from an initial population of planetesimals and pebbles. This work builds on the seminal works of OK10 and LJ12 which presented, in detail, how pebbles are accreted by a single planetesimal or protoplanet. These researchers have found that pebble accretion can be very efficient, making it a promising mechanism to form the cores of giant planets. Indeed, LJ2012 postulated that if there were only a couple of planetesimals large enough to efficiently accrete pebbles, then those planetesimals could directly form the giant planet cores. However, these models have not yet demonstrated whether pebble accretion can be the dominant mechanism of planet growth in a more comprehensive planet formation calculation with many interacting planetesimals.

In this first paper we are interested in investigating a relatively simple scenario, inspired by the aforementioned theoretical work and designed to focus on the process of planetesimal formation by the streaming instability followed by pebble accretion. In this scenario we assume that dust grows effectively and efficiently up to pebble size, leading to a large pebble-to-gas ratio in which a substantial fraction of the solids in the disk are in the form of pebbles. These pebbles then settle toward the midplane and, owing to the oww streaming instability, form clumps that become gravitationally unstable (Youdin & Goodman, 2005; Johansen et al., 2007). Then, because this planetesimal formation is not efficient, there will be a large population of remnant pebbles (Johansen et al., 2009). These pebbles will drift inward and be accreted by the planetesimals, leading to the formation of planets (LJ2012). We furthermore assume that the large pebble-to-gas ratio required by the streaming instability is difficult to achieve and therefore only occurs once in the disk lifetime. As such, we will model a single epoch of pebble formation followed rapidly by a single epoch of planetesimal formation. After the remnant pebbles have either drifted out of the disk or been accreted, we assume that not enough mass is ever again in pebbles in order to impact the planet formation process, through either planetesimal formation or pebble accretion. While this model is simplistic and perhaps naive, as it does not address the challenge of how this large mass of pebbles would have formed in a timescale short compared to the drift timescale, it yields initial conditions that will allow us to directly investigate the problem of pebble accretion as a natural outcome of the leftovers of planetesimal formation.

We note that this work differs from recent work by Chambers (2014) in that we are interested in investigating the initial proposal by LJ12 that pebble accretion can lead directly from planetesimals to giant planet cores. Chambers found that, under certain conditions, pebble accretion may enhance the subsequent planetary growth and evolution if one already has reached the “classical” oligarchic growth stage in planet formation via accretion of planets from planetesimals Kokubo & Ida (1998, 2000), but did not address the evolution of a system in which pebbles are important from the very beginning.

In this paper we utilize the results of OK10 and LJ12 to investigate how planetary systems would evolve from an initial population of pebbles and planetesimals through the formation of giant planets. We begin in Section 2 by discussing the numerical planet formation code we use in this paper. We then present a short review of the physical processes behind pebble accretion in Section 3. We demonstrate the implementation of pebble accretion into our code in Section 4 and then present our fiducial model in Section 5. In Section 6 we briefly discuss the analytical reasons behind the behavior seen in this calculation. In Section 7 we validate the analytical discussion by demonstrating how the disk structure and pebble size impact our numerical results. In Section 8 we discuss the importance of the initial planetesimal population. In Section 9 we then turn to other ideas that have been suggested in the literature to aid core formation with pebbles: shielding, evaporation fronts, and planet traps. In Section 10 we discuss the caveats and implications of our results, and summarize our conclusions.

## 2. Description of our Numerical Model

For this work we have employed a modified version of LIPAD (Lagrangian Integrator for Planetary Accretion and Dynamics), a particle-based Lagrangian code that can follow the dynamical/collisional/accretional evolution of a large number of planetesimals through the entire growth process to become planets. For full details of the code and extensive test suites see Levison et al. (2012), but we summarize the most relevant attributes of the code here.

LIPAD is built on top of the N-body integrator SyMBA (Duncan et al., 1998). In order to handle the very large number of sub-kilometer objects required by many simulations, LIPAD utilizes tracer particles. Each tracer represents a large number of small bodies with roughly the same orbit and size, and is characterized by three numbers: the physical radius, the bulk density, and the total mass of the disk particles represented by the tracer. LIPAD employs statistical algorithms that follow the dynamical and collisional interactions between the tracers. When a tracer is determined to have been struck by another tracer, it is assigned a new radius according to the probabilistic outcome of the collision based on a fragmentation law (Benz & Asphaug, 1999, using the ice parameters , , and ). This way, the conglomeration of tracers represents the size distribution of the evolving planetesimal population. In this work, we do not allow our pebbles to either grow or fragment; therefore, particles below 1 km in size are not involved in the collisional cascade. LIPAD also includes statistical algorithms for viscous stirring, dynamical friction, and collisional damping among the tracers. The tracers mainly dynamically interact with the larger planetary-mass objects via the normal -body routines, which naturally follow changes in the trajectory of tracers due to the gravitational effects of the planets and vice versa. LIPAD is therefore unique in its ability to accurately handle the mixing and redistribution of material due to gravitational encounters, including gap opening, and resonant trapping, while also following the fragmentation and growth of bodies. Thus, it is well suited to follow the collisional cascade of a population of planetesimals while they gravitationally interact with a system of planets.

For this work we have made a number of optional additions to LIPAD. While in principle the basic LIPAD engine can treat the interaction and accretion of pebbles, calculating the trajectories of incoming pebbles as they interact with the gas around the (proto-)planets is not computationally practical in these large-scale calculations. In order to significantly increase the computational speed of the calculation, instead of directly integrating the trajectory of a pebble as it spirals toward the planet, if a pebble is within a planet’s Hill sphere () and satisfies the criteria found by OK10 (which will be described in Equation 2) then we allow the pebble to be accreted. Additionally, although the pebbles can be accreted by planetesimals, we ignore collisions between pebbles.

Another optional addition to LIPAD is gas turbulence. In the absence of turbulence, aerodynamic drag causes these small pebbles to rapidly concentrate at the midplane. Not only are there many mechanisms to drive turbulence in disks (e.g. Balbus, 2000; Lithwick, 2009; Vorobyov & Basu, 2008), but even in the absence of external forcing, as the pebbles settle, they can self-stir the gas disk either by vertical shearing (Goldreich & Ward, 1973; Cuzzi et al., 1993) or by streaming instabilities (Youdin & Goodman, 2005; Youdin & Johansen, 2007; Johansen & Youdin, 2007). We are not modeling the detailed behavior of the gas; therefore, we have added random kicks to the pebbles in order to simulate the effects of turbulence.

Given our resolution, we are unable to resolve any correlated behavior of pebbles due to turbulence; therefore, we use a simple prescription inspired by Youdin & Lithwick (2007) and give the pebbles random kicks of strength

Here is the desired particle scale height and is the Keplerian orbital frequency. This is equivalent to assuming that the correlation time for eddies is the same as our time step (). For our fiducial model we choose such that the pebble density in the midplane is equal to the gas density (consistent with the streaming instability). Following Youdin & Lithwick (2007), we can relate this to the common Shakura & Sunyaev (1973) parameter. For particles near , if the solid-to-gas ratio is 0.01, this corresponds to an . As described in Levison et al. (2012), we calculate the aerodynamic drag on the particles using the formalism of Adachi et al. (1976).

As in previous LIPAD runs, we have used the prescription in Inaba & Ikoma (2003) to include the enhanced capture cross section due to aerodynamic drag. Additionally, as we are interested in the gross evolution of a system after the formation of a potential giant planet core. Therefore, we have added in a simple optional prescription allowing cores to accrete gas envelopes. In order to accrete gas, the core size must be above a critical value, which depends on the mass accretion rate of solids onto the core. We follow Rafikov (2006) to determine whether the current core mass is above the critical value given the current mass accretion rate. If this criterion is met, then we grow the planet of mass by allowing gas to accrete on the Kelvin-Helmholtz timescale (), so that

We follow (Ida & Lin, 2008) and approximate this timescale as

We limit gas accretion to the Bondi accretion rate,

where is the gas density and is the local sound speed. We arbitrarily cut off gas accretion when a planet reaches one Jupiter mass.

Additionally, once a giant planet begins to grow we allow it to open up a gap in the gas disk. A planet opens a gap if

where (Crida et al., 2006), is the gas scale height, is the mass of the Sun, and is the Keplerian velocity. We follow Varnière et al. (2004) and assume that the half-width of this gap is

where is the distance from the Sun to the planet. Inside this region we turn off aerodynamic drag, type I migration, and eccentricity damping.

In most of the simulations presented here we do not allow the growing planets to migrate via type I migration (with the exception of the planet trap simulations in Section 9.3), although we do include type-I eccentricity damping (Papaloizou & Larwood, 2000). In Section 9.3 we present the details for our type I migration for those runs. In all runs we neglect type II migration.

## 3. Pebble Accretion Rates

The physics behind pebble accretion has been discussed in detail in OK10 and LJ12. We encourage the interested reader to find full details within those two papers, but, for convenience, we summarize the key points here.

In order for gas drag to significantly enhance the probability that a particle is accreted by a target object, two criteria must be met. First, the particle’s stopping timescale must be long compared to the timescale for the particle to be deflected by the target object’s gravity, or else it is too well coupled to the gas to be accreted. Second, the particle’s stopping time must be short compared to the time for it to drift past the target, or else the gas drag is not important in the accretion process.

Relevant to the first criteria, OK10 found that, specifically, the gravitational encounter timescale must be shorter than four times the stopping time, i.e.

where is the impact parameter, is the gravitational constant, is the mass of the growing planetesimal, and is the speed at which the particle approaches the planetesimal. This means that in order for a pebble to be captured, it must approach at a distance of

(1) |

Relevant to the second criterion, in the two-body problem a particle is deflected by 90 deg if its impact parameter is

(this is called the Bondi radius in LJ2012, but we distinguish it from the canonical Bondi radius for gas accretion). Let us define the critical crossing timescale as the time for a particle to cross this 90 deg scattering radius

This is effectively the time that a particle spends strongly interacting with the target. If then the particle does not feel the gas drag, and the actual capture cross section is much smaller than the limit given by Equation (1). OK10 found that a pebble is efficiently accreted only if

(2) |

where . This provides quite a sharp cutoff in the mass accretion rate, so that small targets grow only as a result of direct hits, while larger ones readily grow.

In our particle-based numerical model (described in Section 2) we employ Equation (2) to determine when a pebble is accreted by a growing planetesimal or planet. However, we can also integrate over and calculate the expected pebble accretion rate

(3) |

where and are the local Cartesian coordinates in the radial and vertical directions, respectively, and , so that and . is the capture radius which depends on the size of the particles and the relative velocity of those particles.

Aerodynamic drag between pebbles and the pressure-supported disk gas causes pebbles to drift (Adachi et al., 1976; Weidenschilling, 1977a). The relative velocity of a pebble drifting by at a radial distance, , from a planetesimal can be approximated as a sum of this aerodynamic driven drift and Keplerian shear,

(4) |

Here is the planetesimal’s semimajor axis, is the Keplerian orbital frequency, and the parameters and parameterize the amount of pebble drift. We define these in the next two paragraphs.

Due to the pressure () gradient in the disk, the gas generally rotates at a slightly sub-Keplerian velocity of

where

A power-law disk with a surface density profile of

(where ) and a disk scale height of

has

For reference, in the minimum-mass solar nebula (MMSN, Weidenschilling, 1977b; Hayashi, 1981) .

The pebbles have a motion different from both the gas and a Keplerian orbit. It is convenient to parameterize the pebbles using the non-dimensional Stoke’s number The gas drag causes the pebbles to drift relative to a Keplerian orbit, in both the azimuthal and radial directions. The net relative velocity in Equation (4) can be found in the laminar case using

(LJ2012). For pebbles with , can be well approximated by unity.

If the relative velocity is given by Equation (4), the capture radius according to the first criteria can be found analytically by solving the cubic equation

(5) |

In the shear dominated regime (for larger planets), capture can only occur if

(6) |

while in the head wind dominated regime (for smaller planets), capture can only occur if

(7) |

However, for small planets (or planetesimals), the second criteria (and hence the exponential part of Equation (2)) is normally more restrictive.

For reference, in Figure 1 we show example capture radii as a function of planet mass in our fiducial disk model (described in Section 5). Note that for planetesimals larger than the effective capture radius in pebble accretion is substantially larger than the comparable capture radius for planetesimal accretion.

## 4. Initial Tests

To demonstrate the efficiency of pebble accretion to produce large objects quickly, we first place a single embryo at 4.5 AU in a sea of pebbles and allow it to accrete pebbles that drift inward from large heliocentric distances. We assume a surface density for the gas disk of , which is about 2.5 times the MMSN at 1 AU but falls off more shallowly. The gas disk has a scale height of . Pebbles are chosen to have Stoke’s numbers near unity at the location of our planetesimal. We assume solar composition so that the solids (assumed to be all incorporated into pebbles) are 1% of the mass of the gas (consistent with Lodders, 2003). We set the turbulence so that the initial density of pebbles at the midplane equals the midplane gas density.

In Figure 2 we show the growth of an individual planetesimal as a function of time. One can see that if an initial planetesimal is small then there is negligible growth over the time frame shown. However, for initially large planetesimals there is exceptionally rapid growth; a 500 km planetesimal can grow to 10 in 2000 yr! In order to roughly maintain the pebble-to-gas ratio in these simulations as pebbles drift inward, we add new pebbles to the outer edge of the computational domain at a rate of . The pebbles continue to enter the domain for the length of the simulation. In total, over 650 of pebbles drifted to the planetesimal during the simulations. This demonstrates a key aspect of pebble accretion: as long as enough pebbles are present in the disk, a single planetesimal (of sufficient size) can grow extraordinarily rapidly.

We show multiple runs at each initial planetesimal size so that the stochastic nature of pebble accretion in LIPAD is apparent. For most of the runs the growth rates are similar between planetesimals with the same initial mass. The most notable exception to this is one of the initial 300 km particles. In this run the planetesimal happens to accrete a large number of pebbles just before the 2000 yr mark, and is able to grow to 1 . This has to do with the very steep change in the accretion radius for planetesimals in this size range (see Figure 1). The planetesimal happened to grow large enough that its growth rate became significantly larger than other, slightly smaller planetesimals. While this stochastic nature arises from our use of tracer particles, we believe that it is likely representative of real, undoubtedly somewhat inhomogeneous, protoplanetary disks.

In Figure 3 we show the mass accretion rate onto the planetesimals as a function of mass for the above simulations. We compare this to what would be expected by the analytical model given in Equation (3) if the relative velocities are assumed to be those of a laminar disk (Equation (4)). The stochastic nature is apparent, but these simulations demonstrate that at a given planetesimal mass (and for our chosen level of turbulence) the mass accretion rate is well determined, and is dominated by the laminar, not the turbulent, components of the relative velocities. This means that while turbulence is very important for determining the particle scale height, it does not dominate the encounter velocity. Additionally, it is clear that LIPAD is able to successfully reproduce the behavior described in OK10.

## 5. Fiducial Model

As expected, pebble accretion works extremely well when one has an isolated planetesimal sweeping up pebbles. However, in this paper we are concerned with a population of planetesimals. Therefore, for our fiducial model we begin with an initial population of planetesimals between 4 and 10 AU. We draw these planetesimals from a distribution of with a minimum planetesimal radius of 100 km (), a maximum size of 1000 km (), and a total planetesimal mass of 20 . For code stability reasons, we also place a single mass planetesimal on a low-eccentricity orbit at 30 AU. While this particle was intended to not interact with the simulation, it turns out to provide an interesting proxy for the retention of a cold Kuiper Belt and we will return to it later. We use the same disk model as in our test case, assuming a surface density for the gas disk of and that a gas disk has a scale height of . For these simulations we assume that the gas disk remains constant for the lifetime of the simulation.

Following the results of Johansen et al. (2009), we place 20% of the mass in planetesimals, and the remaining is in pebbles. The pebbles have a size distribution of from 0.5 to 5 m. Pebbles are chosen to have Stoke’s numbers near unity, ranging from around to 25 in the inner regions and from to 10 in the outer regions of the disk. Pebbles are given random turbulent kicks consistent with the value that would keep the initial midplane pebble density equal to the gas density. This is similar to a disk with a Shakura & Sunyaev (1973) disk with . We generate new pebbles at 10 AU and allow them to drift inward owing to the aerodynamic drag at a rate of , a rate that maintains the surface density of pebbles in the outer regions. The surface density of pebble tracer particles is shown in Figure 4. For reference, in Table 1 we show the parameters for this fiducial run and for the subsequent runs in this paper.

t | pebbles | pebbles | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

Run | () | (AU) | (years) | pebble size | intermixed | from outside | evaporation | planet trap | ||

Fiducial | 4500 | 0.035 | 1.25 | -1.0 | 0.07 | 0.5–5m | ✓ | ✓ | ||

Fiducial – Run 2 | 4500 | 0.035 | 1.25 | -1.0 | 0.07 | 0.5–5m | ✓ | ✓ | ||

Fiducial Shielding | 4500 | 0.035 | 1.25 | -1.0 | 0.07 | 0.5–5m | ✓ | |||

Large | 17 400 | 0.035 | 1.35 | -1.9 | 0.07 | 0.5–5m | ✓ | ✓ | ||

Small | 1150 | 0.035 | 1.15 | -0.1 | 0.07 | 0.5–5m | ✓ | ✓ | ||

Small Pebbles | 4500 | 0.035 | 1.25 | -1.0 | 0.01 | 0.2–0.5m | ✓ | ✓ | ||

Large Pebbles | 4500 | 0.035 | 1.25 | -1.0 | 0.2 | 5–10m | ✓ | ✓ | ||

Large Shielding | 4500 | 0.035 | 1.25 | -1.0 | 0.2 | 5–10m | ✓ | |||

Snow line | 4500 | 0.035 | 1.25 | -1.0 | 0.07 | 0.5–5m | ✓ | ✓ | ||

Planet Trap | 4500 | 0.035 | -0.5 | -1.0 | 0.07 | 0.5–5m | ✓ | ✓ | ✓ |

We continue generating pebbles for 500 yr, leading to 65 of pebbles entering the domain. This rather high flux illustrates one of the peculiarities in this model. Because pebbles drift inward very rapidly, in order to maintain the pebble surface density at the initial level we must supply them at a high rate. This is a rather high flux of pebbles, and in an actual disk the mass flux of pebbles might be slower and spread out over a longer period of time. However, this would simply slow down the accretion of pebbles relative to the other growth processes going on. Given the fact that we are interested in testing whether pebble accretion can be the dominant mechanism we choose to use these large pebble flux rates.

In Figure 5 we show the initial cumulative distribution of particles in light gray and the distribution after 90, 300, 600, 3000, , and yr of evolution in subsequently darker curves (with the final two curves dashed and dotted, respectively). One can see that even after a very short amount of time the largest planetesimals have grown quite large. As time progresses, the sharp cutoff in mass at the largest size softens and shallows as a few of the planetesimals are lucky and accrete more pebbles than their siblings. However, this runaway growth is not enough to create only a few proto-cores. Instead, by 600 yr almost 200 planetesimals (now protoplanets) grow to above Mars mass in size and none grow to much above an Earth mass. By 3000 yr all of the pebbles have been lost and there are almost 50 protoplanets more massive than 1 , with the largest being around 4.5 .

In the subsequent million years there is some evolution of this profile in which a few of the largest embryos grow to around 10 by eating their siblings, but there are still planets more massive than Earth. Because some of the cores have grown large enough to accrete gas, by 5 Myr three of the cores have become gas giants. At 1 Myr there are still eight non-gas planets larger than 1 and another half-dozen objects between Mars and Earth mass.

Despite the growth of three gas giants, the plethora of Mars- and Earth-mass bodies generated in this fashion provide a significant problem. The reason for this can be seen by looking at the evolution of the embryos’ semimajor axes. In Figure 6 we show the semimajor axis and eccentricity distribution at five different times during the simulation; the initial state and the states after 6000, , , and yr. Early on, the embryos grew extremely rapidly by sweeping up the pebbles. Accretion was so fast that they did not have time to dynamically interact with one another. In the subsequent million years a few of these Earth-mass protoplanets collide and grow, however the vast majority of them scatter and spread instead. This means that even after the few lucky cores that are able to accrete gas form, there are many Earth-mass planets in orbits that are farther from the Sun. We continue the integration of these remaining bodies for a few hundred Myr in a gas-free environment. Eventually these planets become gravitationally unstable, crossing orbits. Some collide and merge, while others are ejected from the system. The initial low-mass planetesimal at 30 AU is completely removed, confirming that this process has catastrophically disrupted the region where today we would expect the Kuiper Belt to have formed. This scenario is clearly incompatible with our own solar system.

Additionally, to speed up the computations, we removed objects that had semimajor axis below 2 AU. To more clearly see what has been lost by this assumption, in Figure 7 we show the semimajor axis evolution of the protoplanets within the computation. At various time steps we indicate Mars-mass, Earth-mass, 10 , and 100 planets with progressively larger points. By looking at the semimajor axis evolution of the protoplanets we can see that there were eight objects greater than 1 that were lost because they crossed into the inner solar system. We would expect many of these (water-rich) bodies to have survived in this region, disrupting any formation of terrestrial planets.

Furthermore, we may be overly optimistic in identifying these final objects as gas giants. Remember that our model for giant planet formation is simple; once the core mass is greater than the critical core mass as a function of its current mass accretion rate we allow it to accrete gas. Subsequent impacts cannot reduce the planet’s gaseous envelopes. This means that despite the fact that the giant planets are impacted with high-velocity collisions with massive embryos, we do not allow erosion of the atmosphere as would be expected (Broeg & Benz, 2012).

Of course, the fact that the oligarchs reached around an Earth mass was a function of the amount of pebbles we injected into the simulations. If we had kept adding pebbles to the disk, then the planetesimals would eventually reach 10 without relying on mutual collisions. However, the relative growth rates would be the same. With no clear protoplanet running away, instead we would have perhaps one hundred 10 cores competing to become gas giant planets, clearly unlike anything within our own solar system. On the other hand, if we had added fewer pebbles then we would have been left with a population of lower mass embryos. This would look very similar to the traditional planet formation scenario except for the fact that most of the mass in the system is in the large oligarchs, so the growth timescale would be even longer than the classical case.

To test the robustness of the result of our fiducial run, we ran another calculation with the same disk and pebble properties, but with different randomly placed planetesimals as our initial conditions. In this second run we saw that the early growth of the planetesimals was very similar to the first calculation (Figure 8). However, the giant planets formed a little closer to each other and scattered each other more strongly. This caused one giant planet to move outward, causing the collisions of other protoplanets (Figure 9). Because we keep the gas surface density constant throughout this simulation, even these later-forming massive cores can accrete gas and become giant planets. In the end this left a system of six relatively low eccentricity Jupiter-mass planets (Figure 10). However, there still are several Earth-mass objects that drifted inward out of the computational domain, and the Mars-mass bodies were scattered throughout the outer solar system. So while this system may look tantalizingly like some known exoplanet systems (such as HR 8799, Marois et al., 2008) it is clearly not what happened in our own solar system.

## 6. Semianalytic Growth Rates

In the fiducial runs it appears that the growth of the planetesimals resembles the standard so-called oligarchic regime in the planet formation process (Kokubo & Ida, 1998, 2000). In this regime the mass doubling time for the largest object becomes long compared to smaller objects, allowing the smaller ones to grow relatively faster and “catch up” to the larger objects. This is in contrast to the “runaway” regime in which a slightly larger object has a shorter mass doubling time, allowing it to grow significantly (and increasingly) faster than its smaller neighbors.

To understand the observed evolution better, we take a closer look at the growth rate of planetesimals due to pebble accretion. Let us define

(8) |

If , then we are in a regime in which the doubling time for the largest (proto-)planets is the shortest, allowing the biggest object to run away from its smaller neighbors. This causes the cumulative size distribution to become shallower. However, if then the largest objects have mass doubling times that are longer than the smaller ones. This means that the smaller objects catch up to the larger ones, causing the size distribution to steepen. At the slope of the size distribution remains unchanged.

By investigating Equation (3), we can semianalytically calculate . In Figure 11 we show as a function of planet mass for the pebble sizes and disk properties in our fiducial model. In these plots we include not only the pebble accretion cross section described in Section 3 but also gravitational focusing as described as the hyperbolic regime in Ormel & Klahr (2010) (their Equation (28)). This is only relevant for the smallest planetesimals. We can see that for small planetesimals , indicating runaway-type growth. However, at around the curve dips below unity, indicating that the growth has become oligarchic-like. This is consistent with Figure 5, in which we see an initial shallowing of the size distribution early in the simulation, but as the pebbles continue to accrete, there is a steepening of the size distribution at high masses as a significant number of objects grow together, as oligarchs.

One way to avoid the problem of having too many Mars- and Earth-mass objects at the time of giant planet formation is to maintain runaway accretion all the way to . To investigate whether this is possible, we can approximate the growth rate as , where is the effective capture cross section and is the density of pebbles at the midplane. Therefore, is determined by a simple combination of the two effects,

Let us first look at the capture cross section. If the scale height of the pebbles () is much larger than the capture radius (), then accretion is fully three-dimensional (3D) and the effective capture cross section is well approximated by . However, if is much larger than , then accretion is effectively two-dimensional (2D) and . This means that

where in the 2D case and in the 3D accretion.

To understand the behavior of as a function of we can look at the two extreme solutions to Equation (2), the shear-dominated regime (Equation (6)) or the head-wind-dominated regime (Equation (7)). In these two extremes,

(9) |

where in the shear-dominated regime and in the drift regime . Since we are interested in the growth rates of larger planets, so long as we limit ourselves to large planets in which

(10) |

then the second term in Equation (9) is negligible. Additionally, so long as , this is the same condition for being in the shear-dominated regime. Therefore, when looking at large planets, we can conclude that .

To understand the behavior of as a function of , we can look at the relative velocity from Equation (4) and take . This leads to

Therefore,

The bracketed part of this equation approaches for small planets (when ), and in the other extreme it approaches . Therefore we can easily conclude that, for a monodisperse pebble distribution (and probably for other distributions as well), at high masses (so that ) the only way to avoid oligarchic growth is to be in the 3D accretion regime (). It is extremely hard to keep pebbles in the 3D accretion regime at high masses. For example, to keep the particle scale height of pebbles larger than the accretion capture radius of a 1 planet in a disk with an , there must be (following Youdin & Lithwick, 2007). So for more realistic turbulence high excitation requires small pebbles. And even in those most favorable cases, where the pebbles are lofted high enough to keep the growth in the 3D regime, we only expect . This is a marginal case, not a strong runaway.

By using Equation 3 we can calculate the parameter as a function of for different pebble sizes and disk structure parameters. In Figure 11 we compare the fiducial model to other disk models (whose disk parameters are shown in Table 1). In red we show the results for a strongly pressure-supported disk (large in Table 1), resulting in a faster drift speed with . Note that at the lowest masses the physical radius of the planetesimal determines the cross section. In blue we show the results for a flatter disk model with (small ) .

We cannot avoid the “oligarchic dilemma” through moderate reasonable changes to the disk parameters. In order to maintain runaway accretion until an embryo reaches 10 , for particles near , we must have . Not only is this value an order of magnitude higher than expected in typical models (such as the MMSN), but in a disk that is so strongly pressure supported the pebbles would drift inward an order of magnitude faster, strongly decreasing the capture probability (as discussed subsequently and shown in Figure 13).

The only way to change the high-mass asymptote of the curve is to keep the growth in the 3D regime. To demonstrate this, in Figure 11 we compare two models with different turbulent strengths such that one is in the 2D regime at large masses, while the other is in the 3D regime. The dashed green curve shows the results for a disk with small pebbles when the disk turbulence is set so that the particle scale height (and thus midplane density) is the same as in the fiducial model. Compare this to the dashed black curve, which shows the result for the same fiducial disk model and small pebbles (parameters in Table 1) with a higher pebble scale height such that protoplanets remain in the 3D accretion regime to high masses. The two curves diverge at high mass. The fact that is slightly less than unity at high mass is due to the small spread in pebble sizes. However, even if had been identically one (as we would expect for a monodisperse pebble population), this would only preserve the size distribution. We cannot preserve runaway growth to high masses, and with a spread of pebble sizes, we cannot even reach the marginal case.

Additionally, in Figure 12 we show that this nearly marginal growth rate comes at a price. The mass accretion rate (for the same pebble surface density) is more than an order of magnitude lower for the excited small pebbles model compared to the fiducial model. In order to not fall deeply into oligarchic growth, the accretion rate is necessarily low. This low accretion rate is problematic because pebbles have a limited life span in the disk. The particles drift inward at a rate of

where

(Weidenschilling, 1977a). We can roughly estimate the probability that a particle is accreted by looking at the ratio of

Figure 13 shows for the various disk models. Although large planets (above 1 ) have fairly substantial capture probabilities, the probability for a pebble to be accreted is quite small for small planetesimals. This means that it takes quite a lot of pebbles to initially grow a planetesimal from, say, to . This is not insurmountable. However, if one looks at a disk with small excited pebbles, then the problem becomes much more acute. Several orders of magnitude more mass must have drifted past the planetesimal without being accreted in order to have the same growth rate as the larger, more settled pebbles. This means that to even just get that almost marginal growth rate, we have to sacrifice the efficiency that makes pebble accretion attractive.

## 7. Basic Model Variations

### 7.1. Disk Structure

In the previous section we can see that simple analytic estimates suggest that the oligarchic growth problem described above cannot be solved by changing the structure of the disk. However, to confirm the validity of the analytic approximations and test whether other physics in the complete LIPAD simulations could possibly lead to a solution, we run a series of simulations where we vary the disk structure.

In Figure 14 we show the evolution of the size distribution for a disk with a larger value (, , which leads to at 5 AU). We have adjusted the disk profile so that the midplane density is the same at 4.5 AU as our fiducial model. Note that this leads to an extremely high surface density at 1 AU (). Thus, this disk is probably not physical and should be considered an extreme case.

Compared to the fiducial model, we see that, during the few thousand years in which the pebbles persist, there is a notably shallower slope in the size distribution between and . This indicates that the planetesimals had a stronger runaway phase that extended to larger planet masses. However the effect is still not strong enough to avoid producing a population of Earth-mass planets that are not incorporated into giant planets (as seen in Figure 15).

Despite the fact that the total amount of solids that enter the computational domain in the above run is the same as in the fiducial case, the total amount of solids retained in the domain after the pebbles disappear is significantly smaller. As we show in Figure 16, the total amount of solids converted into planetesimals varies significantly as a function of the parameter in the disk. In this large model (red curve) only around 40 of pebbles (less than 30% of the pebbles) are accreted by the planetesimals. This compares with the 80% retention in the fiducial model. This agrees with what we would expect from the analytic estimates in Figure 13.

When we look at a less pressure-supported disk (, , , , at 5 AU) with a slower particle drift speed, we see that, except for the very smallest planetesimals, there is little evidence of runaway growth (Figure 17). Indeed, in this model no giant planets are able to form within 5 Myr despite the formation of a couple of 10 cores. Recall that in order to accrete gas the core needs to become larger than the critical core mass, but that the value of the critical core mass is a function of the accretion rate of solids (Mizuno, 1980; Stevenson, 1982). If the solid accretion rate is too high, then the core’s atmosphere is too hot to accrete more gas. In this run, these cores are actually still subcritical owing to the high solid accretion rates, so they do not become gas giants. In Figure 18 we show three snapshots of the disk, demonstrating the large number of excited planetary-mass objects. This shows an additional constraint, not only do we need large embryos, but we also need to shut off solid accretion.

As can be seen in Figure 16, in the low- model (the blue curve) nearly 95% of the pebbles that enter into the domain are accreted by the planetesimals. This high degree of retention is what we expected; the timescale to accrete material is much shorter than the timescale for it to drift, even when the planetesimals are quite small. Therefore, the large- and small- runs confirm that we can get either high retention and many oligarchs, or low retention and (even for an unrealistic disk) an only marginally more favorable size distribution.

All together, the results are as we would expect from the simple analytic model. We conclude that changes to the disk structure cannot, on their own, produce results consistent with our own solar system.

### 7.2. Pebble Size

In Section 6 we argued that for small, more excited particles (i.e., particles that are stirred by outside turbulence rather than self-stirring) it is possible to maintain a growth exponent close to one even at high masses. This is due to the fact that the scale height of the particle remains large compared to the capture cross section of these small particles, keeping the planet in the 3D accretion regime.

To test this effect, we performed a small pebble run where we include pebbles with radii of 20–50 cm and placed them in our fiducial disk, but forced the initial midplane pebble-to-gas ratio to be 0.03 instead of the usual 1.0. This lower gas-to-dust ratio is equivalent to keeping the turbulent parameter the same as in our fiducial run. However, we note that this higher level of turbulence is not directly consistent with efficient planetesimal formation via the streaming instability. In Figure 19 we show the cumulative distribution of planetesimals at different times. The pebbles have all disappeared by 6000 yr (either by being accreted or by drifting inward), and no planetesimal has grown to more than an Earth mass. In the next yr the size distribution does not change notably. The slight steepening of the high mass slope between 600 and 3000 yr indicates that we have still transitioned into the oligarchic growth regime even in this more favorable case. The total amount of growth was much smaller in this run as a result of the low accretion efficiency. In Figure 20 we show that a much smaller fraction of the pebbles are accreted in this simulation than in the fiducial one.

## 8. Initial planetesimal population

Another, perhaps even larger, unknown in these models is the mass and size distribution of the initial planetesimal population. In the previously presented simulations we have assumed that planetesimals form within a size range of 100–1000 km, have a slope of (where ), and have a total mass of 20 . This size distribution was motivated by studies that suggest that the high-mass end of the existing small-body reservoirs in our solar system may be indicative of the initial populations (Nesvorný et al., 2010; Morbidelli et al., 2009). Observationally, is roughly consistent with the present-day slope for objects larger than 50 km (Trujillo et al., 2001; Bernstein et al., 2004; Fuentes & Holman, 2008; Fraser et al., 2008, 2010). The total initial amount of solids in our computational domain was chosen to be 100 , a round number roughly comparable to the total amount of solids currently in the giant planets in our solar system (Hayashi, 1981). We set the percentage of mass in planetesimals to be 20%, which is roughly consistent with streaming instability calculations (Johansen et al., 2009). This mass is also comparable to the mass needed for the postulated so-called Nice disk, the planetesimal disk that was necessary to move the giant planets into their current configurations (Gomes et al., 2005; Tsiganis et al., 2005).

Despite these motivations, very little is actually known about the size distribution of planetesimals expected in the formation process. While theoretical models have shown that it is possible to form planetesimals even as large as via the streaming instability in only a few local orbital periods (Johansen et al., 2007; Chiang & Youdin, 2010), these models do not yet have predictions of the expected size distribution of planetesimals. If planetesimals form with a significantly different size distribution than assumed in this study, one could get very different results. For example, LJ2012 suggested that planetesimals may be produced in such a way that only a few (four or five in our solar system) are large enough to accrete pebbles efficiently, while the rest are too small. In this case the oligarchic dilemma we have described will never arise because there are not enough planetesimals available. However, as we address in detail in the remainder of this section, the LJ2012 model requires a very special set of initial conditions, which we believe may be difficult to achieve.

To aid in our discussion of the initial conditions required by LJ2012, in Figure 21 we show the growth curves for planets in the fiducial disk as a function of amount of material drifting past them (allowing the cores to grow to arbitrarily large sizes). If we want a 1000 km planetesimal to grow to , then around of pebbles must have drifted by it. During that time, any nearby planetesimal much larger than 300 km would also have grown significantly. The tendency toward oligarchic growth means that if a planetesimal is large enough to have grown, then it also will grow to approximately the same mass as all of the objects larger than it. It is important to note that there are no planetesimals with zero growth rates as a planetesimal can always accrete owing to direct impacts, but the timescale for a planetesimal to grow becomes significantly longer as one looks at subsequently smaller initial planetesimals.

Figure 21 shows that it is not simple to define a given mass for which there will be no significant growth; given enough time and pebbles, even small planetesimals will eventually grow. Under the LJ2012 scenario, objects must grow to around , which means that at this location . Additionally, no other objects can grow significantly in the presence of this many pebbles, and thus they all must have initially been smaller than a few times .

If we are interested in an equation to estimate this transition between growth and no growth, we can leverage the power of the exponential cutoff in Equation (2). This suggests a critical cutoff of the form

(11) |

If we choose , this corresponds to the mass at which the mass accretion rate in Equation (2) is reduced by a factor of (shown as a dotted line in Figure 21). The dashed line shows the mass at which the mass accretion rate is reduced by a factor of (). Given that we require little growth for objects below the cutoff when , seems to be a reasonable choice as a critical curve.

In Figure 22 we show the critical curves for various pebble sizes in the fiducial disk model. For each pebble size we show our preferred critical curve (Equation (11) with ). For comparison, these curves are overplotted on the initial planetesimals in our fiducial simulation. There are a number of notable aspects of these curves. First, the critical curve depends strongly on the size of the particle (or more explicitly, ). Particles with large or small stopping times have critical growth masses orders of magnitudes smaller than particles with . This means that if we have a range of pebble sizes then it is actually the extreme pebble sizes that contribute the most to the growth of the smaller planetesimals. It is only once these planetesimals have grown that all of the pebbles are accreted efficiently. Another interesting aspect to note about these critical curves is that, unless there is a transition from the Stokes to Epstein regime (which is the case for m, 1 m), the critical mass increases with increasing distance from the star. We will return to this second observation in a moment.

With these curves in hand, we can evaluate the likelihood of LJ2012’s idea that only four objects will grow. In particular, we constructed Monte Carlo simulations where we placed a population of planetesimals in our fiducial disk and used Equation (11) to determine which are larger than . We first randomly draw planetesimals from various size distributions where we vary the total mass () and the slope (). Then we calculate the number of planetesimals greater than one of these critical curves. Since we have a range of pebble sizes and the critical curves are different for each pebble mass, we define the relevant curve as being the lowest critical mass for all the pebbles in the distribution. We then calculate the probability that a given planetesimal population will have precisely three to five planetesimals greater than this critical curve. In Figure 23 we show the percentage of randomly drawn populations that match the above criteria, if we assume that planetesimals are only formed within a range of 100–1000 km in size. There is a relatively small region of parameter space for which there is a probability of successfully producing a handful of planetesimals that can grow. While it is possible that our solar system was simply a lucky one that happened to be in this range, we find that constraint a little uncomfortable.

Perhaps more importantly, for all reasonable values of the planetesimal disk must be extremely low mass () in order to avoid creating several tens (or more) of large planets. As a comparison, the Kuiper Belt today is estimated to have a mass of 0.04–0.1 (Gladman et al., 2001), and the Oort Cloud is expected to have a mass of at least an Earth mass (Dones et al., 2004). This makes the low mass seem to be a challenging requirement. Therefore, without specific arguments for how a disk would self-regulate to produce the proper mass of planetesimals such that there are only a few planetesimals larger than the critical mass (which depends on both the disk structure and the amount of pebbles drifting through the disk) while still being consistent with the observed small-body reservoirs, we conclude that the idea that only four planetesimals could grow to become giant planets is highly unlikely, although of course we cannot rule it out.

## 9. Additional Model Variations

In the literature there have been several suggestions to enhance giant planet formation by pebble accretion. In the next subsections we address a few of these ideas. We note that we have certainly not exhausted all possibilities related to these ideas, but as the reader will see, none of them have yet obviously presented a solution to the difficulties presented in this paper.

### 9.1. Shielding

As noted by Morbidelli & Nesvorny (2012), when planets get large, they can accrete a substantial fraction of the pebbles that cross their orbits. As a result, they suggested that such a planet might be able to effectively starve planets interior to them. If the larger planetesimals could prevent the growth of their interior neighbors, then this would break down the assumption used in our analytic derivation that all planetesimals have access to the same pebble swarm and potentially lead to a system where only a small number of embryos grow.

Of course, this shielding is natively included in the LIPAD simulations. However, because we begin the simulations with a substantial amount of pebbles intermixed with the planetesimals, it is difficult to ascertain to what extent shielding is occurring. Therefore, in order to specifically investigate shielding, we construct a new set of runs without initial intermixed pebbles and only allow pebbles to drift in from beyond 10 AU.

Our first shielding run is identical to the fiducial model except for the lack of initial pebbles.
In Figure 24(b) we show the distribution of pebbles as a function of time in our new shielding run.
This figure shows how pebbles drift inward from beyond 10 AU.
For ease of comparison, we also show the pebble surface density for the same pebble production, but without any planetesimals in Figure 24(a). It is clear that
some pebbles are accreted as they drift inward, but a substantial fraction of the pebbles still reach the inner regions of the disk.
In the yellow curve in Figure 25 we show the fraction of the pebbles that have been accreted by the planetesimals as a function of time.
By the end, only 67% of the pebbles have been accreted.
A comparison with our fiducial run (black curve) shows that *more* pebbles are accreted there than in the shielding run. This is because the planetesimals in the fiducial run have access to more pebbles in total so they grow larger, and as they grow larger they become more efficient at accreting pebbles.
But by the time the capture efficiency has gotten large, there already is a population of massive oligarchs, so shielding is not effective at alleviating the oligarchic problem in either the fiducial run or our shielding run.

However, as shown in Figure 13 and discussed by Morbidelli & Nesvorny (2012), the efficiencies of pebble accretion depend strongly on the pebble properties. By focusing on the 2D accretion regime, Morbidelli & Nesvorny found that the efficiency of pebble accretion increases for particles with slow radial drift speeds; particles with either relatively large () or small () Stoke’s numbers. If one allows that particles with may be in the 3D accretion regime because they are easily excited to high scale heights in a moderately turbulent disk, then we would expect the capture efficiency to be the highest for large pebbles. This is related to Figure 13; in our simple approximation the growth rate of large planets, when accreting large pebbles, is set by the rate at which those pebbles can be supplied to the growing planet. While in a realistic simulation the capture efficiency is likely not actually unity, we would expect large planets to let very few (if any) large pebbles drift by without being accreted.

To investigate shielding in a more favorable environment, we created runs in the fiducial disk with large “pebbles” (radii =5–10 m). As in the previous shielding experiment, we feed pebbles in only from the outside. Figure 26 shows the distribution of the pebbles in this simulation. As compared to Figure 24, it is clear that fewer pebbles are making it into the inner regions of the disk. We can quantify this by comparing the number of pebbles that have crossed a given annulus with a test run in which we produce the same number of pebbles at 10 AU but have no planetesimals to accrete them. Indeed, over 90% of the pebbles are accreted when there are large pebbles. As seen in Figure 25, the fraction is quite large both in a run in which large pebbles are intermixed (cyan curve) and in the shielding run in which large pebbles only drift in from beyond 10 AU (magenta curve). For clarity we will discuss only this later run from this point forward.

Despite the high pebble accretion efficiencies in the large pebble shielding run, there is no sign of one (or a few) planetesimals running away to become giant planets in the cumulative size distribution (Figure 27). Indeed, the ineffectiveness of shielding is not that surprising if one looks at the strong mass dependence of the filtering fraction as seen in Morbidelli & Nesvorny (2012). Initially, when the planetesimals are small, the filtering fraction of an individual planetesimal is not very large. It is only when the bodies get relatively large (greater than 1 ) that they sweep up a substantial fraction of the pebbles that cross their orbits. In this simulation, no one planet becomes large enough to really “starve” their neighbors before they had a chance to grow.

However, from Figure 24 it is clear that even if no individual planetesimal is capable of preventing the growth of neighbors, as an aggregate population the outer planetesimals do in fact accrete most of the pebbles, preventing the inner ones from growing. Thus, we might expect another possible mechanism by which shielding could aid in the growth of a single core; if all of the growing proto-cores are in the outer annulus of the disk, perhaps they could collide and merge into a single object. However, as the cumulative size distribution clearly shows, this did not happen in our calculation. By looking at the semimajor axis evolution in Figure 28, we can see that though all of the planets did indeed form in the outer region of the disk, these planets scattered each other and spread. This indicates that a high aggregate shielding rate is not sufficient to form a single core; for shielding to be effective, it must be for an individual object. Even under these most favorable conditions the filtering fraction of individual planetesimals is not high enough to stop hundreds of Mars-mass and tens of Earth-mass planets from growing, and spreading, in this disk. It appears that despite the higher accretion efficiency, shielding cannot provide a solution to this oligarchic problem.

### 9.2. Presence of a Snow Line

Another possible way to enhance growth in a localized region is by use of a snow line. As icy pebbles drift inward, they will reach a point in the disk where they will begin to sublimate. The volatiles will then diffuse both inward and outward, but inward-diffusing vapor can only be lost on the gas diffusion timescale, a timescale long compared to the pebble drift timescale while outward-diffusing vapor can be deposited on grains and return to the snow line on the short drift timescale. This can significantly enhance the solid surface density on the outside of the snow line (Stevenson & Lunine, 1988; Cuzzi & Zahnle, 2004). Recently, Ros & Johansen (2013) demonstrated that this effect is particularly effective at enhancing the growth of pebble-sized objects. They demonstrated that condensation fronts can produce a significant amount of pebbles in a local region. They suggested that this mechanism could potentially produce an excess of pebbles, sufficient to allow pebble accretion to dominate in a narrow annulus.

To investigate the above idea, we have run a variant of our fiducial model in which we place an evaporation front at 4 AU. For simplicity, instead of reproducing the detailed pebble growth rates of Ros & Johansen, we include a simple model to roughly mimic their conclusions. We assume that 99% of the pebbles that cross the evaporation front are recycled and reformed at a random distance between 4 and 5 AU. The other 1% of the mass is assumed to have remained in the gas phase and diffused into the inner solar system. In this simulation we do not have pebbles coming in from beyond the computational domain, but we do start with pebbles intermixed between our planetesimals. In Figure 29 we show the number of pebble tracer particles as a function of location and time in the disk. One can see the enhancement of pebbles in the 4–5 AU region. Despite this high retention efficiency, all of the pebbles end up leaving the domain in just over 3000 yr, so overall the pebbles are not around significantly longer than in our standard simulations in which we add pebbles from the distant disk.

In Figures 30 and 31 we see the evolution of the disk with the evaporation front in the same form as in prior figures. The pebbles initially in the outer region drain inward very quickly, causing most of the growth to occur on the inner edge of the domain. However, there is not just one core growing in the inner region; instead, a number of large cores compete to accrete the pebbles around the snow line. These cores mutually excite each other, gaining high eccentricities. This is related to the behavior that we saw in the shielding run; by forcing all of the pebble accretion to occur in a relatively narrow annulus, we do not prevent multiple planets from forming, but we do increase the violence of the dynamical instabilities once the planets have time to interact with one another. In this simulation eventually three giant planets do grow, but the eccentricities of all these planets are above 0.1. This again is not consistent with our own solar system.

### 9.3. Presence of a Planet Trap

In the previous calculations we have neglected type I migration, the migration of embedded planets in a gaseous disk due to tidal interactions with the disk (Ward, 1997). This migration poses a significant challenge to giant planet formation via core accretion. Using the “classical” rate (such as that given in Tanaka et al. (2002)), the migration timescale for a 1 core at 5 AU is less than yr in a MMSN, and decreases with increasing planet and disk mass. This is shorter than the disk lifetime and thus has proved a substantial challenge for core formation. One promising way to avoid losing the cores of the giant planets is to modify the disk structure to slow or even halt type I migration; such a location has been known in the literature as a planet trap (Masset et al., 2006). There are a number of theoretical reasons why planet traps might exist, particularly around the snow line. For example, they could arise as a result of changes in the disk opacity (Baruteau & Masset, 2008; Paardekooper & Papaloizou, 2008; Bitsch & Kley, 2011; Bitsch et al., 2013) or changes in the disk viscosity (Kretke & Lin, 2007, 2012, Bitsch et al. 2014, submitted).

We might expect that including type I migration and a planet trap would aid the core formation process by forcing our proto-cores close together, thus allowing them to merge and grow at a specific location. Therefore, we performed simulations where, in addition to our snow line, we have added type I migration and a planet trap slightly interior to the snow line. We use the migration rate from Papaloizou & Larwood (2000) as described Levison et al. (2012), setting the scaling factor in that paper’s Equation (28), . We place a planet trap at 3.5 AU, interior to the snow line at 4 AU. We have a very simple model of the trap: exterior to the trap type one migration is directed inward as normal, but for a few disk scale heights (0.2 AU) interior to the planet trap we multiply the fiducial migration rate by , so that the planet moves outward. Within 3.3 AU we turn off migration. This is not intended to be a high-fidelity reproduction of any specific planet trap model, but instead is a simple, computationally convenient model for preliminary studies. Additionally, this represents an optimistic case because both planets of all masses are trapped and planets are not forced to continue migration inward if they are interior to the “trap.”

Nevertheless, when we look at the cumulative distribution of the planet trap model (Figure 32), we see that planets never grow beyond a few Earth masses and there is significant loss of material within 5 Myr. After 6000 yr tens of objects greater than 1 have formed, but over the next couple million years of evolution all but two have been lost. The reason for this loss becomes evident when one looks at the semimajor-axis evolution of the particles (Figure 33). Pebbles cause the planetesimals to grow to several Earth masses in the region between the snow line and 7 AU, and these all migrate inward as one would expect. However, once at the planet trap these embryos do not in general merge to form a single planet; instead, they scatter each other. Planets that are scattered outward migrate back toward the trap; however, inwardly scattered planets can be delivered beyond the width of the trap and into the terrestrial planet region where there is no migration. In this inner region planets continue to scatter off of each other (and the planets in the trap) allowing some planets to continue to spread inward. All the while, planets that are scattered outward always return to the trap. As a result, more and more icy planets are delivered to the realm of the terrestrial planets, which is clearly inconsistent with our solar system. As we have placed the inner edge of the computational domain at 1 AU, we lose around 70 of planets/planetesimals. Therefore, while we cannot be sure of the final fate of these objects, it is clear that the presence of a planet trap does not force all of the embryos to merge into a single object in this model.

One may speculate that, if we had allowed type I migration to continue in the inner region, many of these ice-rich planets may have migrated on to the star, leaving a dry terrestrial planet region. While this is possible, it would require two conditions to be met. First, the terrestrial planet building blocks need to not migrate; therefore, they have to be small enough so that type I migration does not remove them and large enough not to drift inward from aerodynamic drag. Second, after the icy planets drifted through exciting their eccentricities, they would need to damp in order to form planets. Therefore, they have to be either large enough to experience type I eccentricity damping or small enough to experience significant gas drag. It is not obvious that there must exist a range of protoplanet sizes that can satisfy both of these criteria and, if there is, that it will likely be small. Therefore, we believe that the disk would have to be very fine-tuned for the terrestrial planets to survive this type of migration.

As a side note, the “leakiness” of planet traps appears likely to be a recurring problem for any model of assembling giant planet cores at planet traps as long as the “building block” components are large enough to scatter each other outside of the trapping region. For traps with widths a couple times the disk scale height, this corresponds to objects with escape velocities the local Keplerian velocity (which is only roughly a Pluto mass at 5 AU).

## 10. Discussion and Conclusions

Given a population of planetesimals and pebbles (objects with stopping times comparable to their orbital times), pebble accretion is potentially an extremely rapid mechanism to grow planets. However, despite the rapid growth rates, we find that invoking pebble accretion as the primary mechanism to form the cores of the giants planets is challenging.

As a test of the viability of pebble accretion for giant planet core formation in our solar system we preform a series of numerical simulations. We follow the dynamical evolution of an initial population of of planetesimals and 80 of pebbles spread between 4 and 10 AU as they are allowed to collide, grow, fragment, and accrete pebbles. We find that pebble accretion efficiently converts planetesimals with initial radii between 300 and 1000 km into Mars-, Earth-, or even larger-mass cores with only a relatively modest amount of pebbles drifting by. However, pebble accretion tends to convert too many planetesimals into these large embryos.

The crux of the issue is that although small planetesimals grow in a runaway manner (in which the mass doubling timescale decreases with planet mass), once protoplanets have entered the regime in which they are most efficiently accreting pebbles, their growth rate looks much more like the oligarchic stage in the standard planet formation scenarios (Kokubo & Ida, 1998, 2000), in which all of the planets grow together. This is due to the fact that their effective capture radius only grows with the Hill sphere, . Additionally, once this large cross section is combined with the fact that pebbles tend to settle toward the midplane on a much shorter timescale than they drift inward (Garaud et al., 2004), accretion in the Hill regime tends to be two-dimensional. Even with the fact that the Keplerian shear brings more material into the capture area of larger planets, the total growth rate only scales as . If the oligarchs reach several Earth masses, then it is likely that a few will collide and allow some to grow large enough to accrete gas, but in addition to these “lucky” cores there will remain a large population of leftover “oligarchs” that violently dynamically interact with each other, spreading throughout the solar system and thereby polluting the terrestrial planet region with too much water-rich material and destroying any possibility of a cold population of planetesimals such as we see in our Kuiper Belt.

Changes to the disk structure cannot easily overcome this problem. As one alters the disk structure in order to push the transition to oligarchic growth to higher masses, one must by necessity decrease the capture efficiency, driving up the needed mass of pebbles to unreasonably high levels. Even if the initial disk was massive enough to create a sufficient quality of pebbles, with this method we require an extraordinarily steep disk profile (, which requires for standard temperature profiles), much steeper than the MMSN or observed disks (Andrews et al., 2010). Moving to smaller, more excited pebbles to keep the growth closer to the runaway regime suffers a similar dilemma in that the efficiency of growth goes down significantly, requiring a very massive pebble reservoir, and still it is not enough to allow a small number of cores to run away to .

Additionally, it has been suggested that perhaps only four objects in our solar system were large enough to efficiently accrete pebbles, meaning that the oligarchic-type growth would not be problematic. We find that a population of planetesimals consistent with the observed small-body reservoirs is very unlikely to have produced only a handful of planetesimals capable of growth (Section 8). Furthermore, changing the parameters to increase the efficiency of shielding (Section 9.1), adding sublimation fronts (Section 9.2), and even adding planet traps (Section 9.3) do not appear to substantially aid the growth of a single (or a few) planetesimals at the expense of their neighbors. Situations in which growth is concentrated in a relatively small annulus still have oligarchic growth within that annulus, and the proto-cores scatter each other out of the preferred region. This demonstrates that the scattering and spreading that we see in the fiducial simulation are quite robust processes and hard to escape. Of course, the calculations in the paper have not exhausted parameter space, and particularly our snow-line model and planet trap model are fairly simplistic. Still, these runs show that it is not trivial to claim that a preferred region of planet formation automatically allows a single (or a few) giant planet cores to form. The fundamental oligarchic nature of pebble accretion tends to win out.

While pebble accretion appears unavoidable if one creates planetesimals via streaming instability, we find it to be challenging to use the sweeping up of leftover pebbles as the primary mechanism to form our own solar system giant planets. Thus the solution to giant planet core formation remains elusive.

## References

### References

- Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756
- Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241
- Balbus, S. A. 2000, ApJ, 534, 420
- Baruteau, C., & Masset, F. 2008, ApJ, 672, 1054
- Benz, W., & Asphaug, E. 1999, Icarus, 142, 5
- Bernstein, G. M., Trilling, D. E., Allen, R. L., Brown, M. E., Holman, M., & Malhotra, R. 2004, The Astronomical Journal, 128, 1364
- Bitsch, B., Boley, A., & Kley, W. 2013, A&A, 550, A52
- Bitsch, B., & Kley, W. 2011, A&A, 536, A77
- Blum, J., & Wurm, G. 2000, Icarus, 143, 138
- —. 2008, ARA&A, 46, 21
- Broeg, C. H., & Benz, W. 2012, A&A, 538, A90
- Chambers, J. E. 2014, Icarus, 233, 83
- Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493
- Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
- Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
- Cuzzi, J. N., & Zahnle, K. J. 2004, ApJ, 614, 490
- Dones, L., Weissman, P. R., Levison, H. F., & Duncan, M. J. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 323, Star Formation in the Interstellar Medium: In Honor of David Hollenbach, ed. D.~Johnstone, F.~C.~Adams, D.~N.~C.~Lin, D.~A.~Neufeeld, & E.~C.~Ostriker, 371–+
- Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067
- Fraser, W. C., Brown, M. E., & Schwamb, M. E. 2010, Icarus, 210, 944
- Fraser, W. C., Kavelaars, J. J., Holman, M. J., Pritchet, C. J., Gladman, B. J., Grav, T., Jones, R. L., Macwilliams, J., & Petit, J.-M. 2008, Icarus, 195, 827
- Fuentes, C. I., & Holman, M. J. 2008, The Astronomical Journal, 136, 83
- Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292
- Gladman, B., Kavelaars, J. J., Petit, J.-M., Morbidelli, A., Holman, M. J., & Loredo, T. 2001, AJ, 122, 1051
- Goldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A, 42, 549
- Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051
- Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466
- Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153
- Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
- Hernández, J., Hartmann, L., Calvet, N., Jeffries, R. D., Gutermuth, R., Muzerolle, J., & Stauffer, J. 2008, ApJ, 686, 1195
- Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419
- Ida, S., & Lin, D. N. C. 2008, ApJ, 673, 487
- Inaba, S., & Ikoma, M. 2003, A&A, 410, 711
- Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475
- Johansen, A., Oishi, J. S., Low, M.-M. M., Klahr, H., Henning, T., & Youdin, A. 2007, Nature, 448, 1022
- Johansen, A., & Youdin, A. 2007, ApJ, 662, 627
- Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75
- Kokubo, E., & Ida, S. 1998, Icarus, 131, 171
- Kokubo, E., & Ida, S. 2000, Icarus, 143, 15
- Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55
- Kretke, K. A., & Lin, D. N. C. 2012, ApJ, 755, 74
- Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32
- Levison, H. F., Duncan, M. J., & Thommes, E. 2012, AJ, 144, 119
- Levison, H. F., Thommes, E., & Duncan, M. J. 2010, AJ, 139, 1297
- Lithwick, Y. 2009, ApJ, 693, 85
- Lodders, K. 2003, ApJ, 591, 1220
- Marois, C., Macintosh, B., Barman, T., Zuckerman, B., Song, I., Patience, J., Lafrenière, D., & Doyon, R. 2008, Science, 322, 1348
- Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478
- Mizuno, H. 1980, Progress of Theoretical Physics, 64, 544
- Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558
- Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18
- Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785
- Ngo, H., Duncan, M., Levison, H., & Thommes, E. 2012, in AAS/Division of Dynamical Astronomy Meeting, Vol. 43, AAS/Division of Dynamical Astronomy Meeting, #05.04
- Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43
- Paardekooper, S.-J., & Papaloizou, J. C. B. 2008, A&A, 485, 877
- Papaloizou, J. C. B., & Larwood, J. D. 2000, Astronomy, 833
- Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62
- Rafikov, R. R. 2004, AJ, 128, 1348
- —. 2006, ApJ, 648, 666
- Ros, K., & Johansen, A. 2013, A&A, 552, A137
- Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Stevenson, D. J. 1982, Planet. Space Sci., 30, 755
- Stevenson, D. J., & Lunine, J. I. 1988, Icarus, 75, 146
- Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257
- Trujillo, C. A., Jewitt, D. C., & Luu, J. X. 2001, AJ, 122, 457
- Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459
- Varnière, P., Quillen, A. C., & Frank, A. 2004, ApJ, 612, 1152
- Vorobyov, E. I., & Basu, S. 2008, ApJ, 676, L139
- Ward, W. R. 1997, Icarus, 126, 261
- Weidenschilling, S. J. 1977a, MNRAS, 180, 57
- —. 1977b, Ap&SS, 51, 153
- Williams, J. P. 2012, Meteoritics and Planetary Science, 47, 1915
- Youdin, A., & Johansen, A. 2007, ApJ, 662, 613
- Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
- Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
- Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57