Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Hughes, Anna
Master of Science - MSc
University of British Columbia
The growth and migration of planetesimals in a young protoplanetary disk is fundamental to the planet formation process. However, in our modeling of early growth, there are a several processes that can inhibit smaller grains from growing to larger sizes, making growth beyond size scales of centimeters difficult. The observational data which are available ( e.g., relics from asteroids in our own solar system as well as gas lifetimes in other systems) suggest that early growth must be rapid. If a small number of 100-km-sized planetesimals do manage to form by some method such as streaming instability, then gas drag effects would enable such a body to efficiently accrete smaller solids from beyond its Hill sphere. This enhanced accretion cross-section, paired with dense gas and large populations of small solids enables a planet to grow at much faster rates. As the planetesimals accrete pebbles, they experience an additional angular momentum exchange, which could cause slow inward drift and a consequent back-reaction on growth rates. We present self-consistent hydrodynamic simulations with direct particle integration and gas-drag coupling to estimate the rate of planetesimal growth due to pebble accretion. We explore a range of particle sizes and disk conditions using a wind tunnel simulation. We also perform numerical analyses of planetesimal growth and drift rates for a range of distances from the star. The results of our models indicate that rapid growth of planeteismals under our assumed model must be at orbital distances inwards of 1 AU, and that at such distances centimeter-sized pebbles and larger are required for maximized accretion. We find that growth beyond 1 AU is possible under certain limited, optimized conditions.
Planetesimal Growth through theAccretion of Small SolidsbyAnna HughesB.Sc., Rensselaer Polytechnic Institute, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2016c© Anna Hughes 2016AbstractThe growth and migration of planetesimals in a young protoplanetary disk is fundamental to theplanet formation process. However, in our modeling of early growth, there are a several processesthat can inhibit smaller grains from growing to larger sizes, making growth beyond size scales ofcentimeters difficult. The observational data which are available ( e.g., relics from asteroids inour own solar system as well as gas lifetimes in other systems) suggest that early growth mustbe rapid. If a small number of 100-km-sized planetesimals do manage to form by some methodsuch as streaming instability, then gas drag effects would enable such a body to efficiently accretesmaller solids from beyond its Hill sphere. This enhanced accretion cross-section, paired with densegas and large populations of small solids enables a planet to grow at much faster rates. As theplanetesimals accrete pebbles, they experience an additional angular momentum exchange, whichcould cause slow inward drift and a consequent back-reaction on growth rates.We present self-consistent hydrodynamic simulations with direct particle integration and gas-drag coupling to estimate the rate of planetesimal growth due to pebble accretion. We explorea range of particle sizes and disk conditions using a wind tunnel simulation. We also performnumerical analyses of planetesimal growth and drift rates for a range of distances from the star.The results of our models indicate that rapid growth of planeteismals under our assumed modelmust be at orbital distances inwards of 1 AU, and that at such distances centimeter-sized pebblesand larger are required for maximized accretion. We find that growth beyond 1 AU is possibleunder certain limited, optimized conditions.iiPrefaceThis Master’s Thesis includes my research conducted from September 2014 through May 2016 inthe University of British Columbia Department of Physics and Astronomy. It presents simulationsand numerical analysis on the conditions of planet formation via pebble accretion. A paper iscurrently in preparation for submission for publication in a journal of astrophysics.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Grain Growth and Spatial Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Going Beyond the Meter Barrier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.1 Concentrating Solids and Direct Gravitational Collapse . . . . . . . . . . . . 41.3 From Pebbles to Planetesimals: Core Accretion and Timescale Issues . . . . . . . . 42 Theory and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Boxzy and Basic Simulation Design . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2.1 Extrapolated Mass Growth Estimates . . . . . . . . . . . . . . . . . . . . . . 163 Idealized Models and Accretion Times . . . . . . . . . . . . . . . . . . . . . . . . . 183.1 Results from Analytical Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27ivList of Tables2.1 Disk and object properties used in the hydrodynamic simulations, following the flared diskprofile outlined in equations (2.13) through (2.15). While most disk conditions change dra-matically at different stellar separations, the relative velocity between the gas and the plan-etesimal stays roughly constant. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Most efficiently accreted pebble sizes in the Epstein regime at each modeled stellar separation,calculated using equation (2.18). Ds is the orbital distance of the planetesimal. . . . . . . . 92.3 Estimated accretion rates for a 100-km planetesimal at densities 1.5 gcm3 (second column)and 3 gcm3 (third column), as well as for a 150-km planetesimal at density 3gcm3 only (lastcolumn) using the planet’s gravitational cross-section as the collisional cross-section. . . . . 142.4 Accretion rates for a 100-km and a 150-km planetesimal at different distances (Ds) rangingfrom 0.1 AU to 10 AU. The 100-km planetesimal results are shown with planetesimal densities1.5 gcm3 as M˙1.5 in column 3, 3gcm3 as M˙3 in column 4, and accretion rates for a 150-kmplanetesimal (only 3 gcm3 ) in column 5. The most efficient accretion rates are bolded wheneverthey are significantly higher for a certain particle size. . . . . . . . . . . . . . . . . . . . . 152.5 Mass growth times for a 100 km sized planetesimal at ρ = 1.5 gcm3 and ρ = 3gcm3 . Thesubscripts on time indicate the planetesimal density (1.5 or 3.0) as well as whether or notthis is the lower estimate (L) on accretion rates (and therefore longer growth times), or theupper estimate (H) on accretion rates (and therefore shorter growth times). . . . . . . . . . 172.6 Mass growth times for a 150 km sized planetesimal of ρ = 3 gcm3 . . . . . . . . . . . . . . . . 173.1 Analytic growth and drift for a 100-km-sized planetesimal of density 1.5 gcm3 and 3gcm3 . Theless dense planetesimal has an initial mass of 6.3× 1021g, while the denser planetesimal hasinitial mass 1.3 × 1022g. The time column gives the time to reach runaway accretion or50 Myr, whichever occurs first. Note that the tabulated masses are found at various timesduring their runaway. Since runaway growth is rapid, this capture a range of final massvalues within the runaway. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.2 Analytic growth and drift for a 150-km-sized planetesimal of density 3 gcm3 . This planetesimalhas an initial mass of 4.241× 1022g. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3 Analytic Growth and Drift Rates for a 1000-km Planetesimal of Density 3 gcm3 . This plan-etesimal has an initial mass of 1.257× 1025g. . . . . . . . . . . . . . . . . . . . . . . . . 21vList of Figures2.1 The Epstein and Stokes regimes for particle size and mean free path as a function of distancefrom the star. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 Distribution of particles on the grid - the higher population at larger y-values takes account ofthe cylindrical geometry. A clear and physical wake has formed to the left of the planetesimaldue to particle and gas flow around it. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Two traced 0.1 cm particles at 10 AU, separated by a distance dy, 10 km. The particle ontop is accreted, whereas the bottom one flows past the planetesimal and is reflected upwards,taking account of particles that would come in from the negative half of the y-axis. . . . . . 122.4 Traced particle tracks at 0.1 AU (left panel) and 1 AU (right panel) for various differentparticle sizes, showing their movement as they are accreted onto or deflected around theplanetesimal. The largest pebble size, 10 cm, comes in from beyond the geometric radius andis efficiently accreted. The most dramatic looping pattern is seen with the 3 cm particles at0.1 AU, which accrete onto the planetesimal even after they have moved past it significantly.At smaller sizes, accretion becomes much less efficient. The 0.0001 cm particles mostlystream past the planetesimal, and are accreted from a smaller area than even its geometriccross-section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.1 Runaway accretion for a 100-km planetesimal at 1 AU in the disk. When the planetesimalreaches runaway accretion, the mass rapidly increases until reaching Miso. In all of the tablespresented in this section, planetesimal growth is recorded as the planetesimal first approachesits runaway or 50 Myr, whichever time comes first. . . . . . . . . . . . . . . . . . . . . . 213.2 Radial growth of the same planetesimal as it hit runaway accretion. The radius goes with3√M , making radial growth much more gradual. . . . . . . . . . . . . . . . . . . . . . . . 223.3 Inward migration is slow and roughly linear. Migration as the planetesimal reaches runawaygrowth deviates slightly from linearity, but does not have the rapid turnoff seen with massand radial growth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22viAcknowledgementsI would like to express my gratitude to my supervisor, Aaron Boley, for his help and guidancethroughout the completion of this project. I would additionally like to thank my professors andcolleagues in the Physics and Astronomy department for their support and encouragement.viiDedicationI would like to dedicate this thesis to my parents - Kathy and Peter Hughes, who have supportedme throughout my academic career.viiiChapter 1IntroductionThe planet formation process is critically dependent on the coagulation of dust grains into largerbodies. During the collapse of a nebular cloud into a star, a circumstellar disk composed of gas andmicrometer-sized grains also forms. It is within this protoplanetary disk that dust grains must, bysome mechanism, conjoin with one another to reach planetary sizes in order to build a planetarysystem. The timescale of growth to planet sizes is thought to be a few million years. This limit ismotivated by several lines of empirical evidence, such as the incidence of gas giant planets, inferredchronologies from meteorites, and infrared observations of protoplanetary disks in nascent systems.In the core-nucleated instability model of planet formation (Pollack et al. 1996), gas giantplanet formation begins with the development of a large solid embryo. As the embryo grows in sizeand mass, a substantial atmosphere is accreted onto it from gas in the disk. If the atmosphere con-tinues to grow before the surrounding gaseous disk dissipates, it can become unstable and undergoa runaway accretion of gas, at which point it can grow into a gas giant. This model implies thatin order to form a gas giant planet, the disk must first develop planetary embryos large enoughto begin accreting atmospheres. In the solar neighborhood, gas giant planets are found aroundroughly 10% of stars within the orbital space between 2 and 2000 days (Cumming et al. 2008). Forcomparison, Jupiter has an orbital period of over 4000 days at an orbital distance of roughly 5 AU.Consequently, the inferred frequency should be thought of as a lower bound for the total occurrencerate of gas giants, as exoplanet detection methods are biased against planets with orbital periodsof Jupiter and longer. The high occurrence of gas giant planets indicates that the mechanism thatforms the embryos must be fairly commonplace, and that the formation time of gas giant planetsalso serves as a reasonable proxy for planet formation timescales in general. As such, a viableplanet formation mechanism must explain growth from micron grains to planets large enough toaccrete a gaseous envelope before the gas dissipates. Observations suggest that the typical gaseousdisk is dispersed after 3 Myr (Haisch et al. 2001, Mamajek 2009). Consequently, significant embryogrowth must be possible within this timescale.Without any detailed observations of circumstellar disks in the beginning phases of planet build-ing, we look to relics within our solar system itself for further constraints. Meteorites and theirparent bodies are remnants of planet formation, containing information about the nature of solidsin the disk and disk conditions during the formation epoch. Some parent bodies underwent signif-icant changes, including differentiation after reaching high internal temperatures; i.e., complete orpartial melting followed by gravitational settling. Meteorite parent bodies that did not differenti-ate preserve information about the conditions of the young solar system during planet formation.These primitive meteorites, called chondrites, provide more pristine information on the primordialcomposition and magnetic field, as well as the distributions of the protoplanetary disk conditions.Moreover, chondrites contain chondrules, which are 0.1-1mm igneous spherules (Desch et al. 2012).As many chondrules have remained largely unaltered since their formation, isotopic dating can beused to infer formation timescales as well as the sequence of formation events. The latter is onlypossible if radiogenic material was homogeneously distributed in the disk, which appears to be the1Chapter 1. Introductioncase (Villeneuve et al. 2009).Isotopic dating further reveals that the oldest chondritic solids in the solar system would haveformed within less than 0.3 Myr of calcium-aluminum-rich inclusions (CAIs) in the protostellardisk, potentially giving the solar system a wealth of solids in the early stages of planet formation(Scott & Krot 2003). CAIs are the oldest known solids in the solar system, dating to 4568.2± 0.2Myr (Bouvier & Wadhawa 2010). Beyond just dating, Scott & Krot (2003) find one chondrite inparticular, Kaidun, may provide unique insight into the planet formation process. Kaidun containsonly millimeter-sized solids, indicating that it may have formed after chondrules were dominant,but before all of the gas in the disk had dissipated. Low abundances of igneous rock in this andsimilar objects could be evidence that these solids first formed inwards of 2 AU and were laterscattered outwards in the disk.Other data indicate that some large planetesimals could have formed much earlier than mostchondrules. Schersten et al. (2006) find that iron meteorites have a narrow range in 182W , and areconsistent with the inferred initial solar system 182W values based on CAIs. This indicates thatthe parent bodies of these iron meteorites must have formed rapidly, within ≤ 1.5 Myr of CAIs.Such rapid planetesimal formation is consistent with recent observations of young protoplanetarydisks.As noted above, infrared observations show that most systems lose their disks within 3-6 millionyears (Hasich et al. 2001), defining a characteristic limit to gas giant formation time. Detailed ob-servations of the protoplanetary disk systems HL Tauri (HL Tau) and TW Hydrae (TW Hya) castlight on the timescales and conditions of planet formation. ALMA (Atacama Large Millimeter/sub-millimeter Array) observations of both systems reveal dark rings, most likely indicative of regionscleared of millimeter grains by large planets. These systems are especially noteworthy because oftheir age; TW Hya is constrained to be between 5 and 10 million years old and HL Tau likely lessthan just 1 million years old (Rhee 2007, Boss et al. 2011). We discuss each of these systems below.TW Hya is a young 0.8-M star with observed dark rings, reminiscent of the effects of planetclearing (Andrews et al. 2016). At a distance of 54 parsecs, it is the closest T Tauri variable star toEarth, enabling high-resolution observations. Its massive circumstellar disk extends to nearly 230AU (Weinberger et al. 2002), and the spatial distribution of solids compared with the gas suggeststhat there may have been significant migration towards the center of the system. Recent ALMAobservations with a spatial resolution of 1 AU reveal bright gaps and bands at sub-millimeter wave-lengths (Andrews et al. 2016). While zonal flows due to magnetization within the disk could beresponsible for the thinner gaps in the disk, it seems unlikely that they could account for the innerring-like structure. The dynamical explanation of these gaps is that planets as small as a few Earthmasses clear material near their orbits. If these gaps are indeed being cleared by planets, this showsplanet formation extends to 8 million years in the case of TW Hya and maybe other systems (Rhee2007).Possibly at the other extreme, HL Tau is a bright young star first identified with far-infraredobservations as hosting a protoplanetary disk (Cohen 1975). HL Tau is the best resolved planetarydisk observed by ALMA, showing very prominent ring structure (ALMA Partnership, 2015). Theeffective temperature and luminosity of HL Tau constrain its age to less than 1 million years, mak-ing it potentially the youngest observed forming planetary system. HL Tau in particular indicatesthat there must be some very rapid mechanism of planet formation (Stephens et. al 2014).21.1. Grain Growth and Spatial EvolutionAll of the above suggests that planet formation can occur within a few million years of the onsetof the disk. Such rapid growth seems to be at odds with classical dust coagulation models. In par-ticular, there is no accepted solution to two major barriers that limit the growth of planetesimals:radial drift and bouncing/fragmentation, discussed next.1.1 Grain Growth and Spatial EvolutionThe early stages of planet formation rely on the growth of large solids from micron to sub-millimeter grains. Even as grains settle into the midplane and stay there on long timescales,collisions between grains will not necessarily result in sizes larger than mm-cm. The outcome ofcollisions is dependent on the properties of the grains themselves. Coagulation through collisionalgrowth can promote the formation of sub-millimeter/millimeter grains, but larger grains have theirgrowth frustrated by bouncing at low collisional speeds and fragmentation at high speeds (Blum& Warm 2008, Testi et al. 2014). Moreover, these larger grains are no longer perfectly coupled tothe gas, and will exhibit non-trivial dynamical evolution according to gravity and gas drag. Forexample, grains will tend to settle towards the midplane (D’alessio et al. 2001), but the likelihoodof settling into the midplane and the settling velocity are dependent on the gas density and par-ticle size (Testi et al. 2014). Midplane settling is opposed by turbulent gas mixing in the verticaldirection of the disk, where small variations over one disk scale height can cause much larger oscilla-tions, bringing dust fragments from collisions into the upper layers of the disk (Dubrulle et al. 1995).Radial drift adds an additional complication to grain evolution. As will be discussed in moredetail in Chapter 2, the gas in the disk is supported in part by a radially outward pressure gra-dient under most disk conditions. This pressure gradient causes the gas to orbit at sub-Keplerianspeeds. In contrast, solids are not supported by gas pressure, and must orbit at Keplerian speeds.The velocity difference between solids and gas, for most conditions, causes particles to encountera headwind and thus gas drag. Subjected to the drag force, particles lose angular momentum andinspiral on timescales proportional to gas density and particle size (Adachi et al. 1976, Weiden-schilling 1977). The resulting inward drift is especially efficient at intermediate object sizes, causingmeter-sized particles to inspiral on the order of a few hundred years, which is faster than they arebelieved to grow through collisions (Birnstiel et al. 2010). This rapid inward drift is known asthe “meter barrier”, as it seemingly frustrates planet formation. Altogether, radial drift, bouncing,and fragmentation make the formation of planetesimals difficult.1.2 Going Beyond the Meter BarrierThe relatively short lifetime of gas in protoplanetary disks means that there must be some rapidmethod by which solids grow to gas giant core sizes within 6 million years. Chondrules also provideevidence that formation must be rapid, and the size distribution of solids in meteorites indicatesthat parent bodies may have preferentially accreted sub-millimeter/millimeter grains. There areseveral mechanisms that could cause the collapse of such smaller solids to form a planetesimal,introduced below.31.3. From Pebbles to Planetesimals: Core Accretion and Timescale Issues1.2.1 Concentrating Solids and Direct Gravitational CollapseIf solids can be concentrated to a sufficiently high density regardless of the size of the solids,a planetesimal could form by direct gravitational collapse, owing to the self-gravity of the solidconcentration (Goldreach & Ward 1973). There are a number of processes that can enhance thesolid density in different areas throughout the disk. One is direct settling of solids, but this iscountered by turbulence, including turbulence produced by shear and particle-gas feedback (Jo-hansen et al. 2009). Turbulence is also able to sort solids by size and concentrate them further.Non-axisymmetric spiral arms can produce pockets of high pressure and density, and reversals inthe pressure gradient draw solids into these pressure and density maxima (Lodato & Rice 2004,Haghighipour & Boss 2003). At first, this would rely on spiral structure being generated withouta planet, such as by gravitational instabilities in the gas (Darisen et al. 2007). In older disks, gapscaused by orbiting giant planets could sort solids into rings of high particle density. Once solids aretrapped, increased particle density causes more collisions between them, speeding up the growth ofsmall solids in the region (Rice et al. 2014).An interesting mechanism that is self-starting from low-amplitude noise in the disk is the stream-ing instability. Streaming instability models rely on solid-gas feedback changing the flow of the gasin the disk. If enough solids are concentrated into local regions due to, e.g., local pressure max-ima, then the solids’ collective behavior can reduce or effectively eliminate drag on the aggregate(Youdin & Goodman 2005). This process can allow inward-drifting solids to be collected in thelocal overdense region until it becomes gravitationally unstable. This could provide the disk witha population of 100-km sized objects, such as those seen in the Kuiper Belt and Trans-Neptunianpopulations. However, the growth of larger planetesimals or the planets themselves cannot be ad-dressed by this alone (Pollack et al. 1996). Moreover, the instability appears to be very sensitiveto metallicity even though this is not clearly seen in the observed exoplanet distribution for planetssmaller than approximately Neptune size (Buchhave et al. 2012).The above discussion is not meant to be exhaustive. Rather, it highlights some potential waysa few planetesimals could form from a large population of micron grains. If a few rare planetesi-mals are able to form in the disk despite impediments, preferential accretion of small solids couldfacilitate rapid growth of even larger planetesimals.1.3 From Pebbles to Planetesimals: Core Accretion andTimescale IssuesThe most commonly accepted model of planet formation is core accretion by pairwise collision.In core accretion models, planets grow directly from collisions between the population of largeplanetesimals. This model assumes that a swarm of meter- to kilometer-sized objects form in thedisk along with a few larger, heavier cores. The larger cores grow by accreting the smaller solidsin the disk, at a rate ofM˙ = ρsvrelR2Fgrav, (1.1)41.3. From Pebbles to Planetesimals: Core Accretion and Timescale Issueswhere vrel is the relative velocity between the solids, R is the core radius, and Fgrav is the gravita-tional focusing factor,Fgrav = 1 +(vescvrel)2. (1.2)The gravitational focusing factor increases the effective cross-section of the planetesimal by a factorof Fgrav due to the planetesimal’s gravity.In the case of low relative velocity, the planetesimal is able to grow very rapidly through run-away growth (Kokuba & Ida 1999). Embryos eventually grow to reach their isolation mass,Miso ≈ 5× 1020((rAU)2 Σkgm−2) 32(M M?) 12kg (1.3)where Σ is the local surface density of solids, and beyond which growth by runaway accretionstops due to insufficient material in the disk. For a roughly minimum-mass solar nebula at 1 AU,Weidenschilling (1977) estimates Σtot = 42, 000kgm2for the total surface density (dust and gas), orΣ ∼ 420 kgm2if the dust to gas ratio is 0.001. If there is significant mixing in the disk, however,through perturbations smaller objects may migrate into the feeding zone of the embryo, enablingit to grow beyond this limit.The pairwise collisional model faces a few setbacks to reaching significant growth, however.Collisions between intermediate-sized objects are more likely to result in fragmentation than coag-ulation into a larger object (Leinhardt & Stewart 2011). Inward drift of fragmented pieces can causea serious depletion of solids, giving larger planetesimals nothing to accrete and making growth tosuch sizes unlikely. Beyond that, the timescales of building planetesimals in this scenario can beexceedingly long (which we highlight in Chapter 3) seeming to outlast the lifetime of gas in thedisk. This would prohibit the formation of gas giant planets. At large distances from the star, thetime it takes to build a core for a gas giant planet can be as long as 10-20 million years (Teiser &Dodson-Robinson 2013; Fortier et al. 2007). This is inconsistent with the observed and predictedtimescales of gaseous disks.The pebble accretion model of planet formation is an alternate to the pairwise collisional model.It relies on some mechanism such as the streaming instability to rapidly form an initial distributionof rare core planetesimals. Instead of growing planets through collisions between many large plan-etesimals, pebble accretion relies on these few planetesimals in the disk to sweep up the potentiallyabundant small, co-orbiting solids. Drag forces between the small solids, or pebbles, and the gasin the disk can change the trajectories of the pebbles as they encounter the planetesimal, causingthem to accrete from beyond the planetesimal’s gas-free gravitational cross-section. This enablesrapid growth of the planetesimal in certain regimes.The focus of this work is to investigate the conditions under which pebble accretion is optimizedand by directly simulating accretion under different conditions, as well as investigating the long termtimescales of growth and inward drift of planetesimals. First we examine in detail the theoreticalframework behind pebble accretion Chapter 2. We then present the results of our hydrodynamicsimulations of pebble accretion under various conditions, and briefly discuss some of the implicationsof such results. Finally in Chapter 3 we use numerical analysis to estimate the rates of growth andinward drift of planetesimals over long timescales, as well as to assess the the effects of gravitationalfocusing with and without gas on the time to reach isolation mass.5Chapter 2Theory and MethodsGas in protoplanetary disks is expected to generally orbit at sub-Keplerian speeds. The orbitalmotion of the gas is mainly set by the potential of the star, which would entail a circular orbitat Keplerian speeds for orbiting solids. However, the gas pressure gradient causes a non-negligibleperturbation to the orbital motion, leading the gas to behave as if the effective gravity wereg′ = −GM?r2− 1ρgdPdr, (2.1)where r is the stellar separation, M? is the stellar mass, ρg is the gas density, the gas pressureis P, and dPdr is the pressure gradient of the disk. This equation is a consequence of hydrostaticequilibrium. Assuming the gas orbit is circular, the effective gravity must be balanced by centrifugalacceleration,g′ = −r(vgas)2, (2.2)− GM?r2− 1ρgdPdr= −r(vgas)2. (2.3)The pressure-supported circular speed is thusvgas =√GM?r+rρgdPdr. (2.4)In contrast to the gas, solids do not have additional pressure support and must orbit at Keplerianspeeds. For dPdr < 0, solids will encounter a headwind, creating an interaction between the solidsand gas through drag. Small solids quickly become entrenched with the gas flow, while large plan-etesimals carry enough momentum that gas drag only affects their motion on very long timescales.The acceleration of a particle due to gas drag can be writtenv˙ = − 1τf(v − u), (2.5)where v is the particle velocity, u is the gas velocity, and τf is the friction time (Whipple 1972).The friction time is determined by the particle size and local gas conditions. There are two differentregimes for small solids that characterize their interactions with the gas. A pebble is in the Epsteinregime when its radius (assuming a spherical shape) is much smaller than the mean free path ofthe gas, and in the Stokes regime (i.e., the fluid limit) when the radius is larger than the mean freepath. In the Epstein regime, the friction time isτf =sρscsρ, (2.6)where s is the pebble radius, ρs is the internal density, cs is the adiabatic sound speed, and ρ is thegas density. In the fluid limit, the equation used for the stopping time depends on the Reynoldsnumber,Re =2sρvη, (2.7)6Chapter 2. Theory and Methodswhere η is the gas viscosity (Weidenschilling 1977). For Re < 1, which is valid for many conditionsin planetesimal formation, the Stokes drag law is applicable withτf =49s2ρscsρgλ(2.8)where λ is the mean free path of the gas (Whipple 1972). The friction time takes additional formsfor Re > 1, but for simplicity we do not discuss these cases further. The mean free path can befound usingλ =m¯ρσ(2.9)where m¯ is the mean particle mass, and σ is the cross-section of molecular hydrogen. The divisionbetween particle sizes for the Epstein and Stokes regimes is indirectly dependent on the stellocen-tric separation due to the pressure and temperature dependence, which on average, decrease withstellocentric separation. For the average molecule radius, we assume a value corresponding to H2,σ = 2.0 × 10−15cm2 (Rafikov 2004). The pressure can be found assuming ideal gas conditions,P =Rgµ Tρ. For 1 AU, P = 1.08 Pa and the mean free path is λ = 1.82 cm. Figure 2.1 showsthe separation between regimes as a function of stellocentric separations following our assumptionsabout the gas and disk profile.Figure 2.1: The Epstein and Stokes regimes for particle size and mean free path as a function of distancefrom the star.Small particles rapidly lose angular momentum to the gas, but because τf is so low for them,radial drift becomes negligible. Large solids have such high τf values that they only lose angu-lar momentum on very long timescales. In between these sizes, when τfΩ is approximately unity(where Ω is the angular frequency), solids lose angular momentum efficiently. For these sizes, the7Chapter 2. Theory and Methodsstopping time is long enough (but not too long) that significant radial drift of the particles can occur.Assuming the planetesimal moves in a circular, Keplerian orbit, it has velocityvkep =√GMr. (2.10)Since the gas is sub-Keplerian, the planetesimal will experience a relative wind speedvrel =√GMr−√GMr+rρdPdr. (2.11)To evaluate different positions in the disk, we assume that the gas distribution is similar to typicalprofiles for the protostellar disk used.For midplane conditions, we use the following surface density, gas density, temperature, andpressure profiles,Σ = Σ0(rAU)−1.5. (2.12)ρ = ρ0(rAU)−2.5, (2.13)T = T0(rAU)−0.5, (2.14)P =RgµTρ, (2.15)where r is the stellar separation measured in AU, ρ0 = 1×10−9 gcm3 and T0 = 300 K are gas densityand temperature values at 1 AU. The gas conditions at various stellar separations are shown inTable 2.1, along with the corresponding relative wind speeds. For our assumptions, the wind speedvrel is roughly constant.Stellar Distance (AU) Gas Density gcm3Temperature (K) Relative Velocity ms10 3.2× 10−12 95 54.83 6.4× 10−11 170 54.71 1.0× 10−9 230 54.70.3 2.0× 10−8 550 54.60.1 3.1× 10−7 950 54.6Table 2.1: Disk and object properties used in the hydrodynamic simulations, following the flared disk profileoutlined in equations (2.13) through (2.15). While most disk conditions change dramatically at differentstellar separations, the relative velocity between the gas and the planetesimal stays roughly constant.Particle accretion or deflection is dependent on the gas properties as well as the particle size.In the absence of any gas, the particles would move unimpeded and accrete onto the planetesimal,based on the geometric area enhanced by the gravitational focusing factor defined in equation (1.2).The presence of gas in the disk means that the pebbles will feel a drag force as their trajectories arealtered by the encounter with the planetesimal. Here the pebble size is important for determiningits behavior in the gas, for which larger particles are minimally influenced by the gas and smaller8Chapter 2. Theory and Methodsones are strongly influenced.To quantify a pebble’s response to gas drag as it is diverted around the planetesimal, we intro-duce the stopping distance, which is estimated using the friction timescale and the relative windspeed, i.e.,dstop = vrelρsρgscs, (2.16)where cs is the adiabatic sound speed. Accretion is optimized when the stopping distance iscomparable to the size of the Bondi radius of the object onto which the pebbles are accreted. TheBondi radius is defined here asRb =Gmv2rel. (2.17)The distance Rb is roughly where the kinetic energy of the incoming pebbles becomes comparableto the gravitational potential due to the planetesimal. If a significant fraction of a pebble’s kineticenergy is dissipated while it traverses the volume defined by Rb, then the pebble can be accretedonto the planetesimal. Efficient accretion occurs when dstop ∼ Rb, as the pebble is not perfectlycoupled to the gas and still experiences drag during its encounter with the planetesimal. The mostefficient size can be estimated by equating dstop with Rb for the Epstein regime, which yieldss =ρgGmρsvrel3√γRgTµ, (2.18)where γ = 1.4 is the adiabatic index of the gas and µ = 2.3 gmol is the mean weight of the gas.In practice, the Bondi radius cannot exceed the Hill radius of the planetesimal. This condi-tion is not revealed in the simulations explored here, as the planetesimal is never large enough toapproach this limit. The Hill radius defines the region of gravitational influence in a three-bodyconfiguration. It is the effective area of a smaller body in orbit around a host body, as defined fora massless particle.Using equation (2.18), Table 2.2 shows the ideal particle size for each stellar separation, for ourassumed disk.Ds (AU) Pebble Radius (cm)10 0.00013 0.0061 0.120.3 3.20.1 65Table 2.2: Most efficiently accreted pebble sizes in the Epstein regime at each modeled stellar separation,calculated using equation (2.18). Ds is the orbital distance of the planetesimal.In our research, we investigate the accretion rates of a range of pebble sizes at all of the stellarseparations in Table 2.2. In order to model planetesimal growth at different distances, we use aCartesian grid with particles streaming past a centered planetesimal. The movement of gas past92.1. Boxzy and Basic Simulation Designthe planetesimal is modeled by a wind tunnel, where the speed of oncoming gas is set by the rel-ative velocity between the planetesimal and gas at the specified stellar separation. Particles areintroduced onto the grid at velocities equal to the relative gas speeds, adapting to friction withthe gas. All simulations are run over short timescales until the accretion rates are stabilized torecover instantaneous accretion rates in a number of conditions and pebble sizes. Over such shorttimescales, we do not consider inward radial drift of the planetesimal or inward drift of the pebbles.2.1 Boxzy and Basic Simulation DesignWe ran our simulations using Boxzy Hydro (Boley et al. 2013), a Cartesian, self-consistentfinite volume hydrodynamics code. Boxzy includes a detailed equation of state with rotational andvibrational H2 states, self-gravity, direct particle integration, and drag terms with feedback on thegas. Boxzy solves the hydrodynamics continuity equations in the following form:∂ρ∂t+ ~5 · (~vρ) = 0 (Mass Continuity),∂ρ~v∂t+ ~5(~vρ~v) = −~5P − ρ~5Φ (Momentum Equation),∂E∂t+ ~5 · (~v(E + P )) = −ρ~v · ~5Φ (Energy Equation),where E = 12ρ|v|2 + , for a gas density ρ, velocity v, and internal energy density .The code uses a Reimann solver for shock capturing. The hydrodynamic state is advanced byfirst solving the hydrodynamics equations using the current state to predict the state at the end ofthe time step, and then uses that prediction to update the final state. The inclusion of gravity issolved using a kick-drift-kick-like method.The gas is evolved adiabatically, comprised of a mixture of H2 (mass fraction 0.73), He (massfraction 0.25), and metals (mass fraction 0.02) such that µ = 2.3 gmol . The rotational and vibra-tional modes of H2 are included in the Equation of State tables presented in Boley et al. (2007).A fixed 3:1 ortho-parahydrogen mixture is used. The full gas conditions using the above equationsare shown below in Table 2.1, for various stellar separations.Boxzy can be run in 2D, 3D, or 2D-3D. We use the latter configuration in this study. The2D-3D simulations use a cylindrical cross-section to capture 3D effects on a 2D grid, creating awind tunnel around the planetesimal, set in the frame of the planetesimal. The wind is the relativevelocity between the gas and the planetesimal.Our simulations use a grid of 256 × 128 cells, where each cell has the dimensions of 10 km × 10km or 15 km × 15 km, depending on the planetesimal size, so that the radius is always resolved by10 cells. We use a 100-km- and a 150-km-sized planetesimal, horizontally centered on the grid asshown in Figure 2.2. In the 2D-3D setup, we exploit azimuthal symmetry about the x-axis, makingthe y-axis the radial direction.102.1. Boxzy and Basic Simulation DesignParticles are coupled to the gas using a cloud-in-cell approximation and a kick-drift-kick inte-gration algorithm. The gas drag and gravity kicks are applied while sourcing the gas. All particlesthat hit the x-axis are reflected back upwards with the same velocity, taking account of pebblesthat would enter the grid from below the planetesimal. This effect can be seen in the bottom panelof Figure 2.3.Pebbles are modeled by “superparticles”; i.e., each particle on the grid represents a swarm ofreal particles. The grid is considered a thin wedge of a cylindrical volume surrounding the planetes-imal, such that the placement of the particles must take this geometry into account. Each wedgecontains more particles at higher y-values, as seen in Figure 2.2. Consequently, the particles areplaced on the grid with a bias towards higher y-values. This is done by placing particles according to,Y = (ny)× (dy)×√P ′, (2.19)where ny is the number of cells in the y-direction, dy is the length of each cell, and P ′ is somerandom variate between 0 and 1. When projecting onto a rectangular grid, this simply places agreater number of particles towards larger y-values (or larger radial coordinates) onto the grid.Figure 2.2: Distribution of particles on the grid - the higher population at larger y-values takes account ofthe cylindrical geometry. A clear and physical wake has formed to the left of the planetesimal due to particleand gas flow around it.Of the 100,000 particles, Boxzy traces in detail the positions and velocities of 60 particles,stacked in the y-direction from 0 to beyond the Bondi Radius. This allows us to plot the motionsof certain individual particles and observe if they are accreted (lower panel of Fig. 2.3, to theleft of the planetesimal) or if they bypass the object (lower panel of Fig. 2.3, to the right of the112.2. Simulation Resultsplanetesimal).Figure 2.3: Two traced 0.1 cm particles at 10 AU, separated by a distance dy, 10 km. The particle on topis accreted, whereas the bottom one flows past the planetesimal and is reflected upwards, taking account ofparticles that would come in from the negative half of the y-axis.Whenever any particle (traced or not) makes contact with the planetesimal, it is accreted,recorded, and reset to enter along with the incoming wind according to equation (2.19). Distribut-ing the particles as above means that each accretion event already takes account of the cylindricalgeometry, allowing the actual accretion rate to be directly calculated from the simulation’s countof particles accreted per unit time. The particle-in-cell algorithm effectively spreads a particle overthe nearest 4 cells (in 2D-3D mode), with a distance-like weighting. Whenever any part of the cloudcontacts the planetesimal, the entire superparticle is accreted. This slightly enhances the physicalsize of the planetesimal, but is taken into account in our analysis.From the output of total accreted particles per unit time, ∆N∆t , it is straightforward to use thegrid’s cylindrical geometry to find the accretion rate, m˙, as measured in kgs . This requires that eachsuperparticle represents a fixed amount of mass in solids. The mass of each superparticle is thenMs = ρfpiLy2LxNpart, (2.20)where ρ is the unperturbed gas density, f is the dust-to-gas ratio, here assumed to be 0.01, Ly isthe simulation length in the y-direction (radial), and Lx is the length along the symmetry axis.The accretion rate is then simply m˙ = ∆N∆t Ms.We can use this to further estimate growth timescales. In particular the mass growth time, orthe time needed to increase the planetesimal mass by a factor of e, assuming a fixed accretion ratetm =moriginalm˙. (2.21)2.2 Simulation ResultsWe ran a total of 132 simulations using the 2D-3D Boxzy configuration at distances rangingfrom 0.1 AU to 10 AU and particle sizes ranging from 1 µm to 10 cm, for a 100-km planetesimal122.2. Simulation Resultswith density 1.5 gcm3, a 100-km planetesimal with density 3 gcm3, and a 150-km planetesimal withdensity 3 gcm3. In all cases, particles are assumed to have an internal density of 3 gcm3, which isneeded for the drag calculations.The traced particles (tracks shown in Figure 2.4) show the path taken by the innermost pebbleson their trajectory towards a planetesimal for select cases. At 1 AU, the gas density is ρ =1× 10−9 gcm3, and the most efficiently accreted pebble size is ∼0.3 cm. Pebbles much smaller thanthe most efficient size are largely deflected for the low dust-to-gas ratio assumed, with an effectiveimpact radius that is smaller than the gravitational cross-section.Figure 2.4: Traced particle tracks at 0.1 AU (left panel) and 1 AU (right panel) for various different particlesizes, showing their movement as they are accreted onto or deflected around the planetesimal. The largestpebble size, 10 cm, comes in from beyond the geometric radius and is efficiently accreted. The most dramaticlooping pattern is seen with the 3 cm particles at 0.1 AU, which accrete onto the planetesimal even afterthey have moved past it significantly. At smaller sizes, accretion becomes much less efficient. The 0.0001 cmparticles mostly stream past the planetesimal, and are accreted from a smaller area than even its geometriccross-section.Table 2.4 shows accretion rates from simulations assuming a 100-km planetesimal with densities1.5 gcm3and 3 gcm3and a 150-km planetesimal with density 3 gcm3. The gravitational focusing factorchanges between cases, as well as the analytical accretion rates assuming accretion only from thegravitational cross-section. The gravitational focusing factor, defined in equation (1.6), depends on132.2. Simulation Resultsthe escape velocity of the planetesimal and the relative velocity of incoming material. For constantdensity, equation (1.6) can be written,Fg = 1 +83piGR2ρs(vrel)2. (2.22)Assuming a relative velocity of approximately 54ms , Fg = 3.88 for ρ = 1.5gcm3, and Fg = 6.75 forρ = 3 gcm3.The analytic accretion estimates are found usingm˙ = fρgpiR2vrelFg. (2.23)As noted, Boxzy uses the cloud-in-cell method, in which individual particles are integrated ontocontinuous phase space. In our grid, each cell represents 10 kilometers, so that the planetesimal’seffective radius in the simulations includes an additional 10 km cell, or 15 km for 150-km planetes-imal. This gives accretion rates in Table 2.3. For large particle sizes, the simulation accretion ratesand the gravitational focusing rates converge, as expected. This also shows that our simulationboxes are large enough to capture the accretion trajectories.Ds (AU) M1.5kgs M3.0kgs M150kgs10 250 440 20003 5100 8900 4.1× 1041 8.0× 104 14× 104 6.4× 1050.3 1.6× 106 2.8× 106 13× 1060.1 210× 105 420× 105 2.0× 108Table 2.3: Estimated accretion rates for a 100-km planetesimal at densities 1.5 gcm3 (second column) and3 gcm3 (third column), as well as for a 150-km planetesimal at density 3gcm3 only (last column) using theplanet’s gravitational cross-section as the collisional cross-section.We run simulations for each stellar distance and a range of particle sizes, tabulated below for100-km sized planetesimals of density 1.5 gcm3and 3 gcm3as well as a 150-km sized planetesimal ofdensity 3 gcm3. The larger planetesimal more dramatically highlights enhancement of growth atefficient pebble sizes, making the effects of pebble accretion more pronounced.142.2. Simulation ResultsTable 2.4: Accretion rates for a 100-km and a 150-km planetesimal at different distances (Ds) ranging from0.1 AU to 10 AU. The 100-km planetesimal results are shown with planetesimal densities 1.5 gcm3 as M˙1.5 incolumn 3, 3 gcm3 as M˙3 in column 4, and accretion rates for a 150-km planetesimal (only 3gcm3 ) in column5. The most efficient accretion rates are bolded whenever they are significantly higher for a certain particlesize.Ds(AU) Particle Size (cm) M˙1.5(kgs ) M˙3(kgs ) M˙150(kgs )10 3 250 430 19001 250 430 19000.3 250 430 19000.1 250 450 20000.03 250 470 35000.01 250 460 66000.003 240 820 84000.0001 120 250 940Ds (AU) Particle Size (cm) M˙1.5(kgs ) M˙3(kgs ) M˙150(×104 kgs )3 10 5100 8700 3.93 5200 9100 4.11 5100 9500 6.80.3 5000 9400 130.1 4800 15000 170.03 6600 20000 120.01 5600 13000 5.70.003 2500 5500 0.440.0001 150 250 0.084Ds (AU) Particle Size (cm) M˙1.5(×104 kgs ) M˙3(×104 kgs ) M˙150(×105 kgs )1 10 7.9 15 173 7.7 20 261 8.8 31 230.3 10 25 120.1 5.7 12 4.90.03 2.1 4.5 1.60.01 0.79 1.6 0.610.0001 0.13 0.13 0.028Ds (AU) Particle Size (cm) M˙1.5(×104 kgs ) M˙3(×104 kgs ) M˙150(×105 kgs )0.3 100 160 280 12010 160 360 5003 200 490 2301 54 110 420.3 15 31 110.1 6.0 12 4.20.03 3.2 5.2 1.40.01 2.6 2.9 0.680.0001 2.5 2.5 0.42Ds (AU) Particle Size (cm) M˙1.5(×105 kgs ) M˙3(×105 kgs ) M˙150(×106 kgs )0.1 100 260 470 25030 250 470 56010 240 830 7903 290 640 2801 63 130 460.3 8.1 13 4.70.1 3.8 3.9 6.20.03 3.7 3.7 0.640.01 3.6 3.7 0.620.0001 3.6 3.6 0.62 152.2. Simulation ResultsFor most distances, there is at least one pebble size with accretion rates exceeding those pre-dicted using gravitational focusing alone (shown in table 2.3). This enhancement is most efficient atsmall stellar separations, and least efficient at large distances. At 10 AU, only the 3 gcm3planetesi-mals accrete from significantly beyond the gravitational cross-section. At closer distances, accretionbecomes much more efficient and many particle sizes have accretion rates enhanced by larger thana factor of two. The smallest planetesimal (of density 1.5 gcm3) at 10 AU accretes all but micrometerpebbles sizes at roughly the same rate, so that accretion onto such planetesimals would reflect thesize distribution of particles in the disk, as no one size is preferentially accreted. At closer distances,accretion rates are optimized for particles of a few centimeters in size, consistent with expectationsfrom Table 2.2. Much smaller particles have considerably lower accretion rates, even a few ordersof magnitude lower than estimates from gravitational focusing. This indicates that particles mustundergo growth beyond sub-micrometer regimes before being effectively accreted onto a growingplanetesimal in the presence of gas.Solids much larger than the preferred size have accretion rates approaching what would be ex-pected from gravitational focusing alone. This is due to very weak coupling between large particlesand the gas during the encounter with the planetesimal. Numerically, this further provides animportant check on the simulation volume, which the results demonstrate to be sufficient for thisstudy. For example, if the accretion rates for decoupled particles were notably smaller than theexpected gravitational focusing rates, then the simulation volume would need to be lengthened.Using the results for the larger particles at 10 AU, all planetesimals are have accretion rates within5% of the expected rates.The largest planetesimal shows more exaggerated effects of pebble accretion. Simulation resultsfor the 150-km planetesimal have enhanced accretion rates for most particle sizes, and larger factorsof enhancement. For the most efficient particle sizes at 0.1 AU, the accretion rate enhancement isa factor of 4 × larger than that due to gravitational focusing. If the initial planetesimal starts outlarger, then pebble accretion is able to produce much higher instantaneous accretion rates.Accretion rates at 0.1 AU tend to exceed those at 10 AU by up to 5 orders of magnitude, whichis mainly due to the disk gas density law. This highlights the perceived difficulty in planetesimalformation at large stellar separations and is emphasized further by extrapolating our accretion ratesinto mass growth estimates.2.2.1 Extrapolated Mass Growth EstimatesThe mass growth timescale, tm, given in equation (2.23), uses the simulation results to estimatethe time needed for a planetesimal to increase in mass by a factor of e, assuming a constant ac-cretion rate. We calculate mass growth timescales for all three simulated planetesimals: 100 km at1.5 gcm3and 3 gcm3, and 150 km at 3 gcm3for all distances. For each orbital distance and planetesimal,we calculate an upper and lower limit on the mass growth timescale by selecting the highest andlowest m˙ values found in our simulations. The tabulated mass growth times are presented below,first in the case of a 100-km planetesimal and then for a 150-km planetesimal.162.2. Simulation ResultsTable 2.5: Mass growth times for a 100 km sized planetesimal at ρ = 1.5 gcm3 and ρ = 3gcm3 . The subscriptson time indicate the planetesimal density (1.5 or 3.0) as well as whether or not this is the lower estimate (L)on accretion rates (and therefore longer growth times), or the upper estimate (H) on accretion rates (andtherefore shorter growth times).Ds (AU) t1.5;H (years) t1.5;L (years) t3.0;H (years) t3.0;L (years)10 2.4× 108 1.7× 109 4.9× 108 1.6× 1093 3.6× 107 1.4× 109 2.0× 107 1.6× 1081 2.8× 106 1.5× 108 1.3× 106 3.0× 1080.3 1× 105 8.0× 106 8.2× 104 1.6× 1070.1 7800 5.5× 105 4800 1.1× 106Table 2.6: Mass growth times for a 150 km sized planetesimal of ρ = 3 gcm3 .Ds (AU) tH (years) tL (years)10 1.6× 108 1.4× 1093 7.9× 106 1.6× 1091 5.2× 105 4.8× 1080.3 2.7× 104 3.2× 1070.1 1700 2.2× 106At large stellar separations, even the high estimates of these timescales exceed the lifetime of gasin the disk. At short distances, however, both low and high timescales show rapid growth. There isa clear cutoff at distances greater than 3 AU, where planetesimal growth slows dramatically. Thee-folding times at such large distances indicate that planetesimal growth at far stellar separationsis unlikely. Instead, planetesimal growth through pebble accretion requires close stellar separationsin order to reach embryo sizes within the expected timescales.17Chapter 3Idealized Models and Accretion TimesOur simulations from Chapter 2 run over short timescales, providing only instantaneous accre-tion rates. These rates can be used to project long-term growth of the planetesimal as done insection 2.2.1 for assumptions regarding the typical particle size. However, these estimates do nottake into account the changing mass of the planetesimal and consequently, the evolving accretionrates due to a larger cross-section and increased gravitational and drag focusing. Simplified nu-merical analysis can provide insight into the long-term evolution of the planetesimal in the disk,evolving the accretion rates as the mass increases. Moreover, the analysis can be used to evaluateinward migration of the planetesimal. This migration is due to at least two effects. The planetesi-mal moving through the sub-Keplerian gas experiences a drag-induced torque, causing inward drift.The solids accreted by the planetesimal are also sub-Keplerian, and thus cause the planetesimal topreferentially accrete slightly low-angular momentum material. The latter is most effective whenthe stopping distance of the incoming material is larger than the planetesimal’s Bondi radius. Forthe analysis here, we include these two torques and track the long-term migration and growth ofplanetesimals, for planetesimals starting out in the given size ranges.First consider the torque on the planetesimal due to drag. We assume to good approximationthat at any given moment the instantaneous mass is fixed. The drag torque is then given by~Γd = ~Fd × ~r, (3.1)~Γd = −12cDρpiR2v2relreˆz, (3.2)where Fd is the drag force, cD is the drag coefficient (assumed to be ≈ 0.5 for a sphere in the plan-etesimal limit), R is the planetesimal radius, and r is the stellar separation. We use the conventionthat the headwind is negative for vrel. We further restrict motions to be within the midplane ofthe disk, which is perpendicular to the eˆz direction, as noted in equation (3.2). Hereafter, we dropthe explicit vector notation, with it being understood that we are exploring radial motions in themidplane, and that the relevant torques and angular momenta are in the ±eˆz direction.Given a constant mass, the time-derivative of the orbital angular momentum is,ddt(m(GM?r)12 )|d = 12mvkepr˙|d = L˙d (3.3)where m is the instantaneous planetesimal mass, which again is held constant for gas drag alone,and vkep is the circular Keplerian velocity. We have further assumed that the planetesimal remainson a circular orbit. Since torque is just the change in angular momentum,L˙d = Γd, (3.4)12mvkepr˙|d = −12cDρpiR2v2relr. (3.5)18Chapter 3. Idealized Models and Accretion TimesSolving for the drift rate due to drag, r˙d,r˙d = −34cDρρsrRv2relvkep, (3.6)where the planetesimal (solid) density is demarcated ρs.The drag torque is complemented by an accretion torque, Γa, directly due to the accretion ofpebbles. In this case,Γa = rm˙(vrel + vkep). (3.7)Here, we can no longer hold the mass constant; thus,ddt(m(GM?r)12 )|a = m˙vkepr + m2vkepr˙|a. (3.8)Again equating this expression with the accretion torque and solving for r˙|a,r˙|a = 2m˙mrvkepvrel (3.9)in the case of an accretion-only process. When vrel = 0, there is no migration due to accretionbecause the specific angular momentum does not change.To evaluate the growth of the planetesimal, we require an estimate for m˙. To move forwardwith this simplified model, we consider only the case of gravitational focusing. We will, however,return to this issue and explore how pebble accretion under the Bondi limit would enhance thegrowth rates. For gravitational focusing only,m˙ = −fρpiR2vrelFgrav = −fρpiR2vrel(1 +v2escv2rel), (3.10)where Fgrav is the gravitational focusing factor, f is the dust-to-gas ratio (0.01), and vesc is theescape velocityv2esc =2GmR=8piGρsR23. (3.11)Substituting Equation (3.11) into (3.10) yields,m˙ = −fρpiR2vrel(1 + 8piR2ρs3v2rel), (3.12)which can be used to arrive at the accretion drift,r˙|a = 32fρρsrRvrel[1− vrelvkep](1 +8piR2ρs3v2rel). (3.13)In this framework, the total drift rate r˙ is the sum of the drift due to drag and accretion, r˙|d + r˙|a.We use a fourth order Runge-Kutta integrator to solve for the mass growth (under equation 3.3) andthe migration rates, using timescales of up to 50 Myr to see long-term growth and drift. Growth anddrift are co-dependent, which is why we choose to solve the growth of the planetesimal numerically.As we show next, the inward migration is small, so we can solve the growth equation analyticallyfor two limits. We first give the results of our numerical model because it is more flexible, and can193.1. Results from Analytical Workbe easily altered to fit a range of disk conditions and planetesimal sizes.The Runge-Kutta integrator requires a disk model for input. We assume the same disk param-eters as in the wind tunnel simulations, as given in Table 2.1. We follow planetesimals with aninitial size of 100 km and densities of 1.5 gcm3and 3 gcm3. We also explore planetesimal sizes of 150km and 1000 km, each with an internal density of 3 gcm3. In all cases, we assume a dust-to-gas ratioof 0.013.1 Results from Analytical WorkTables 3.1 and 3.2 show the growth and migration results for a 100-km and a 150-km planetes-imal respectively, for five different locations in the disk. For reference, the detailed growth as afunction of time is shown in Figures 3.1-3.3. In this case of a 100-km planetesimal with density1.5 gcm3at 1 AU, mass growth accelerates considerably after 8 Myr and has clearly runaway by 10Myr. Planetesimals only reach runaway accretion within a few Myr at short stellar separations (<1 AU). Recall that this is calculated using gravitational focusing only, and thus further emphasizesthe need for enhanced accretion rates.Table 3.1: Analytic growth and drift for a 100-km-sized planetesimal of density 1.5 gcm3 and 3gcm3 . Theless dense planetesimal has an initial mass of 6.3 × 1021g, while the denser planetesimal has initial mass1.3× 1022g. The time column gives the time to reach runaway accretion or 50 Myr, whichever occurs first.Note that the tabulated masses are found at various times during their runaway. Since runaway growth israpid, this capture a range of final mass values within the runaway.Initial Distance Time (years) Final Mass (×1022 g) % Decrease in Distance(AU) 1.5 gcm33 gcm31.5 gcm33 gcm31.5 gcm33 gcm310 5 ×107 5 ×107 0.62 1.3 0.47 0.253 5 ×107 5 ×107 2.0 3.6 4.4 2.51 1 ×107 5 ×106 730 1.7 6.7 2.00.3 5 ×105 5 ×105 1.1× 104 600 3.8 2.00.1 3.5×104 3× 104 1.0× 106 210 3.0 1.2Table 3.2: Analytic growth and drift for a 150-km-sized planetesimal of density 3 gcm3 . This planetesimalhas an initial mass of 4.241× 1022g.Initial Distance (AU) Time (years) Final Mass ×1023(g) % Decrease in Distance10 5× 107 0.45 0.23 5× 107 0.20 2.01 5× 106 8.3 2.00.3 1× 105 1.1 0.40.1 1× 104 20 0.4We repeat these calculations for much larger initial planetesimal in order to see exaggeratedeffects on growth and drift, using the same gas conditions. In this case it should be possible to seethe planetesimal reach 10 M⊕, the size of Jupiter’s solid core, within an appropriate timescale.203.1. Results from Analytical WorkTable 3.3: Analytic Growth and Drift Rates for a 1000-km Planetesimal of Density 3 gcm3 . This planetesimalhas an initial mass of 1.257× 1025g.Initial Distance (AU) Time (years) Final Mass ×1025 (g) % Decrease in Distance10 1× 107 1.4 0.13 1× 107 12 1.51 1× 106 180 0.430.3 5× 104 210 1.00.1 3800 4× 104 1.2For a planetesimal of an initial size of 1000 kilometers, rapid growth to gas giant core sizes ispossible within 3800 years at the closest orbital distance, and on the order of a few tens of thousandsof years at 0.3 AU. Beyond 1 AU, accretion slows again, restricting even the larger planetesimalfrom reaching gas giant size within 1 Myr. For most simulations, the planetesimal undergoes themost rapid growth and migration at close orbital distances.Figure 3.1: Runaway accretion for a 100-km planetesimal at 1 AU in the disk. When the planetesimalreaches runaway accretion, the mass rapidly increases until reaching Miso. In all of the tables presented inthis section, planetesimal growth is recorded as the planetesimal first approaches its runaway or 50 Myr,whichever time comes first.213.1. Results from Analytical WorkFigure 3.2: Radial growth of the same planetesimal as it hit runaway accretion. The radius goes with 3√M ,making radial growth much more gradual.Figure 3.3: Inward migration is slow and roughly linear. Migration as the planetesimal reaches runawaygrowth deviates slightly from linearity, but does not have the rapid turnoff seen with mass and radial growth.The point at which the planetesimal reaches a runaway can be thought of in terms of growthof the planetesimal’s radius. This limit can be found using the m˙ in equation (3.12), using R′ asthe effective radius, yielding,m˙ = fρvrelpiR′2. (3.14)For a planetesimal of constant internal density ρs, m˙ can be written 4piR2R˙ρs, where R is thegeometric radius of the planetesimal. Rearranging equation (3.14) in terms of the planetesimal223.1. Results from Analytical Workradius,R˙ =1RfρρsvrelR′2R2. (3.15)This expression has two limits for R′: the gravitational focusing limit and the Bondi limit. Inthe gravitational focusing regime, R′ is defined as the gravitationally enhanced radius only, notconsidering any outside effects such as gas drag. The Bondi radius, however, takes into accountboth gravitational focusing as well as gas drag effects. As we saw in the wind tunnel simulations,accretion is highly dependent on the size of the particles being accreted. First consider the grav-itational focusing limit, as was explored in the above numerical integration. We will discuss theBondi radius next.In the gravitational focusing limit, R′ = R√1 + v2escv2rel. Using the expression for the escapevelocity in equation (3.11), R˙|G is the planetesimal growth for gravitational focusing,R˙|G = 14fρρsvrel +2pifρGR23vrel. (3.16)The first term in equation (3.16) is the contribution from the geometric component. As gravita-tional focusing becomes more important, the second term dominates while the first term becomesnegligible. Equation (3.15) is of the form R˙ = AR2 +B. The solution isR(t)|G = vrel√38piGρstan[fρ√piG6ρst + arctan(R0vrel√8piGρs3)]. (3.17)Runaway occurs when the tangent argument approaches pi2 ; thus,t|G = 1fρ√6ρspiG[pi2− arctan(R0vrel√8piGρs3)], (3.18)where t|G is the runaway time in the gravitational focusing limit. Since migration is small, theanalytic growth is consistent with our numerical integrations. Some important scalings are notimmediately clear from this form, so we return to equation (3.16), but ignore the geometric term,solving the equation in the limit of large gravitational focusing. In that case, the approximatesolution to R˙(t)|G is,R˙(t)|G ≈ 3vrelR03vrel − 2pifρGR0t , (3.19)which has a clear runaway when t = 3vrel2pifρGR0 .Thus, as planetesimals approach the strong focusing limit, the runaway time is inversely pro-portional to the initial size of the planetesimal itself.The second limit is set by the Bondi radius rather than the gravitational focusing radius,R′ = Gmv2rel. Again using equation (3.15),R˙|B = 4pi2fρρsR4G29vrel, (3.20)233.1. Results from Analytical Workwhere the subscript B denotes the Bondi regime. The underlying assumption is that all particlesthat enter the Bondi region will be accreted due to a combination of gravity and drag, which onlyapplies to a particular evolving size range. Again integrating the R˙ expression, the planetesimalradius as a function of time isR(t)|B =(3v3relR303v3rel − 4pi2G2fρρsR30t) 13. (3.21)As in the gravitational focusing limit, runaway accretion happens as the denominator approacheszero:t|B = 3v3rel4pi2G2fρρsR30. (3.22)Notice that this timescale depends inversely on R30. For each regime, the time it takes for theplanetesimal to reach runaway accretion can be solved for taking a particular stellar distance andinitial planetesimal size. Taking the distance to be 1 AU and using corresponding ρ and vrel valuesfrom Table 2.1, and a planetesimal with initial radius R0 = 100 km and ρs = 3gcm3, t|B = 2.8× 106years and in the gravitational focusing regime, t|G = 12.2× 106years.Recall that these times are for the runaway in radius. Since the mass has an R3 dependence, theeffective mass runaway times are shorter, as highlighted in Figures 3.1 and 3.2. Enhancement dueto gas drag decreases the time required for a planetesimal to reach runaway growth in some regimesand the initial size of the planetesimal is extremely important for determining the magnitude ofenhanced accretion due to drag effects. For example, at 150 km for R0, the ratio becomest|Gt|B ∼ 10,all other things being equal.24Chapter 4DiscussionIn this work, we have investigated pebble accretion as a method of rapid planetesimal growthin protoplanetary disks. Our hydrodynamic simulations show that for nearly every position in thedisk and corresponding disk conditions, there is a preferentially accreted pebble size. The opti-mized accretion rates of such sizes are bolded in the comprehensive results shown in Table 2.4.Centimeter-sized pebbles are most efficiently accreted at 1 AU and inwards, while sub-millimetersizes are preferred at large stellar separations. Chondrule-sized objects in particular are prefer-entially accreted between 1 AU and 3 AU, which overlaps with the location of many observedmeteorite parent bodies. The size of preferentially accreted pebbles decreases with increasing stel-locentric separation, as is expected and exhibited in Table 2.2. By 10 AU, the low-density 100-kmplanetesimal shows no preference in accreted pebble size and only a minimal aversion to micronpebbles. We thus expect small planetesimals that formed at large stellar separations to roughlypreserve the size distribution of solids at the time of formation in its accreted material. However, asthe initial mass of the host planetesimal increases, sub-millimeter particles are strongly preferencedat 10 AU, as can be seen in the accretion rates of both the high-density 100-km planetesimal aswell as the 150-km planetesimal. Based on the correlation between distance and accretion ratesseen in our simulations, we anticipate that even large initial planetesimals accrete inefficiently atlarge stellar separations.Solids much larger than the preferred size have accretion rates approaching what would be ex-pected from gravitational focusing alone due to weak coupling with the gas, while small solids arepredominately excluded from accretion due to strong coupling with the gas as it streams around theplanetesimal. This effect can be dramatic in some regimes, and can significantly reduce the effectiveaccreting planetesimal area to be much smaller than its geometric cross-section. This transition tolow accretion rates can be rapid as well, with efficiency sharply dropping off with decreasing particlesize. If planetesimals do grow via pebble accretion, then the particle size distribution that can berecovered from an undifferentiated planetesimal will exhibit the following: (1) an optimal pebblesize for an initial mass and formation location, (2) a decrease in pebble filtering for pebbles largerthan the optimal size, and (3) a low presence of pebbles smaller than the optimal size. This couldpotentially be seen in samples retrieved by the upcoming asteroid-bound OSIRIS-REx (Origins,Spectral Interpretation, Resource Identification, Security - Regolith explorer) mission.Our simulations and numerical analysis both allow us to explore planetesimal growth timescales.As a consequence of decreasing accretion rates with increasing stellocentric separation, planetesimalgrowth timescales are exceedingly long outwards of 1 AU. The mass growth timescales tabulated insection 2.7 exceed the gas lifetime at distances beyond 0.3 AU for both of the 100-km planetesimals,and beyond 1 AU for the 150-km planetesimal, taking only the high estimates on accretion rates.All but one of the low accretion estimates exceeds the gas lifetime, indicating that even at closeseparations growth by pebble accretion requires particles at optimal sizes. These values do assumeconstant accretion rates, which may overestimate growth timescales. The numerical analysis up-dates accretion rates with growing mass, taking account of long term growth and feedback with25Chapter 4. Discussioninward drift of the planetesimal.The results of our numerical analysis confirm that while growth at close stellar separations issubstantially faster than at large separations, there are a few unlikely venues by which a planetesi-mal could grow substantially even at 10 AU. Under the gravitational focusing limit, planetesimalsat 0.1 AU to 1 AU can reach runaway growth within 1 Myr, while growth remains slow and steadyat 10 AU, even for a 1000-km planetesimal. These timescales can be drastically shortened in theBondi limit, where the time to runaway scales with 1R3, as shown in equation (3.22). The Bonditimescale to reach runaway growth is 4.4 × faster than the gravitational focusing timescale for a100-km planetesimal, 10 × faster for a 150-km planetesimal, and over 300 × faster for a 1000-kmplanetesimal. This expands the possibility of runaway growth to include 3 AU for the 150-kmplanetesimal, and up to 10 AU for the 1000-km planetesimal. Along with the simulation results,this shows how the addition of gas drag can significantly raise accretion rates and lower growthtimescales, facilitating potential for growth under a wider range of conditions.Both the simulations and numerical analysis show that pebble accretion functions best at closestellar separations, and that pebble accretion is very limited, but in principle possible, at largeseparations. Each stellocentric distance has a preferred pebble size, which accretes with significantlyhigher rates than expected from gravitational focusing alone. Growth by pebble accretion requires alarge population of solids of efficient sizes in order to capitalize on gas drag. In optimized conditions,pebble accretion can be very rapid and cause the planetesimal to reach runaway growth well withinthe time constraint set by gas giants, meteorite chronologies, and observed protoplanetary disks.We make a number of assumptions that could be altered to fit a range of disk models. The dust-to-gas ratio of 0.01 ascribed for the midplane may be a low estimate, while larger values due togreater midplane settling or higher metallicity could enhance accretion. The assumed disk profilescould additionally be tweaked to fit a number of models. Our simulations do not take inward driftof planetesimals into account, which would lead to a slight non-axisymmetric accretion of pebbles,with a preference towards the outward-facing side of the planetesimal. These deviations may alteraccretion rates slightly or impart an additional torque force on the planetesimal due to preferentialaccretion, but the distribution of thick gas at close stellar separations ensures that the planetesimalgrows most efficiently close to the star, with limited flexibility.26BibliographyAdachi I., Hayashi C., & Nakazawa K. 1976, PThPh, 56, 1756Andrews S. M., Wilner D.J., Zhu Z., et al. 2016 ApJ 820, 2Bai X., & Stone J. 2010 ApJ 722.2, 1437-459.Birnstiel T., Ricci L., Trotta, F., et al. 2010, A&A, 516, L14Blum J. & Wurm G., 2008. Annu. Rev. Astron. Astrophys. 46, 21-56.Boley A. C., Morris M. A.,Ford, E.B. 2014 ApJ 792, 2Boley A. C., Morris M. A., Desch, S. J. 2013 ApJ 776.2, 101Boss A. P. 2011 ApJ, 731, 74Bouvier A., & Wadhwa M. 2010 Nature Geoscience 3.9, 637-41Buchhave L. A., Latham D., Johansen A. et al. 2012 NatureCohen M. 1975 ApJ 343, 393Cumming A., Butler R. P., Marcy G. W., et al. 2008 PASP 120, 531D’Alessio P., Calvet N., & Hartmann L. 2001 ApJ 553, 321-334Desch S. J., Morris M. A., Connolly H. C., et al. 2.12 AP 47, 1139Dubrulle B. 1995 Icarus 114, 237-46Fortier A., Benvenuto O. G. , & Brunini A. 2007, A&A 473.1, 311-22Goldreich P. & Ward W. R. 1973 Astrophys. J. 183, 10511062Haisch K. E., Lada E. A., & Lada, C. J. 2001 ApJ 553.2Eunkyu H., Wang S. X., Wright J. et al. 2014 Publications of the Astronomical Society of thePacific 126.943, 827-37Haghighipour N. & Boss A. P. 2003 ApJ 583, 996100327BibliographyHayashi C., Adachi I & Nakazawa K. 1976 Progress of Theoretical Physics 55.3, 945-46Hayashi C. 1981 Fundamental Problems in the Theory of Stellar Evolution 113, 28Kokubo E. & Ida S. 1998 Icarus 131.1, 171-78Kothe S., Bloom J., Weidling, R. et al. 2013 Icarus 225.1, 75-85Leinhardt Z. & Stewart S. 2011 ApJ, 745.1, 79Lodato G., Rice W. K. M. 2004, MNRAS 351, 630Mamajek E. E. 2009, in American Institute of Physics Conference Series, Vol. 1158, AmericanInstitute of Physics Conference Series, ed. T. Usuda, M. Tamura, & M. Ishii, 310Morbidelli A., Bottke W., et al. 2009 Icarus 204.2, 558-73Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J. et al. 1996, Icarus 124, 62Rafikov, R. 2005 ApJ, 621.1Rhee J. H., Song I., Zuckerman B. et al. 2007 ApJ 660, 1556Rice W. K. M., Lodato G., Pringle J. E. et al. 2004 MNRAS 355,543Schersten A., Elliott T., Hawkesworth C. et al. 2006 Earth and Planetary Science Letters 241.3-4,530-42Scott E. & A. Krot. 2014 Treatise on Geochemistry, 65-137Stephens I. W., Looney L. W., Kwon K. et al. 2014 Nature 514.7524, 597-99Teiser J., & Dodson-Robinson S. E. 2013 A&A 555Testi L., Birnstiel T., Ricci L. et al. (2014) ArXiv e-prints.Villeneuve J., Chaussidon M., & Libourel G. 2009. Science 325, 985988Weidenschilling S. J. 1977 Astrophys & Space Sci 51.1, 153-58Weinberger A. J., et al. 2002 ApJ, 566, 409Whipple F. L. 1972 The Motion, Evolution of Orbits, and Origin of Comets 401-08Youdin A., & Goodman J. 2005 ApJ 620.1, 459-6928
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Planetesimal growth through the accretion of small solids
Text