On the growth efficiency of pebble accreting planetesimals

(1)

On the growth efficiency of pebble

accreting planetesimals

A bachelor thesis presented to the faculty of science of the university of Amsterdam by:

Rico Visser rico.visser@student.uva.nl

Amsterdam, the Netherlands, 28 of June, 2015

Daily supervisor : Dr. Chris Ormel first examiner : Prof. dr. Carsten Dominik

Second examiner : Dr. Anna Watts

(2)

Nederlandse samenvatting

Je kent vast wel dat moment van schoolkamp vroeger of op vakantie, dat je ’s avonds even naar de heldere sterrenhemel keek en jezelf vragen stelde als; wat is daar nou eigenlijk verder, behalve de aarde? Welk wonder heeft ervoor gezorgd dat wij hier als mens terecht zijn gekomen? Tenminste, dat vroeg ik mij wel vaak genoeg af. Dit zijn vragen die van nature voortkomen uit onze nieuwsgierigheid en men heeft deze vragen door de historische tijdlijn heen proberen te beantwoorden. Als we op-timistisch zijn hebben wij het ook al best ver geschopt. Het zonnestelsel is aardig in kaart gebracht, we hebben velen planeten en sterren ontdekt buiten ons zonnes-telsen en we weten zonder al teveel moeite een karretje op mars te zetten. Toch is dit niet bevredigend voor de echte wetenschapper. Deze wil alles tot in de puntjes begrijpen. Dit begint bij het allerkleinste, de deeltjes die fungeren als bouwstenen van het heelal zoals wij die kennen, tot de gigantische constructies die hieruit voort zijn gekomen zoals hemellichamen, zonnestelsels en sterrenstelsels. Door de jaren heen zijn voor al deze vakgebieden modellen en theorieen opgesteld om de fysische processen te begrijpen en blijkt dat er veel vragen nog onbeantwoord blijven.

Een van de vragen waar men zich al jaren lang over buigt is hoe planeten nou vormen. Het is niet zo makkelijk om deze vraag te beantwoorden, omdat het labo-ratorium waar deze processen zich afspelen niet binnen handbereik is. In dit thesis wordt gekeken naar een mogelijk scenario om planeetformatie te verklaren. Hiervoor is het wel zo handig voor de lezer om te weten hoe ik te werk ga. Om te begin-nen is het belangrijk om te weten hoe wij denken dat een ster ontstaat. Sterren zijn fascinerende objecten die geboren worden uit gaswolken die niet meer in even-wicht kunnen blijven. Alle deeltjes in de gaswolk bij elkaar worden als het ware naar het centrum van de wolk getrokken door de zwaartekracht. Hierdoor wint de zwaartekracht van de druk die naar buiten staat door het gas omdat de wolk nog re-latief koel is en stort de wolk in. Hierbij komt een enorme hoeveelheid energie kijken wat natuurlijk een gigantische explosie veroorzaakt. Een ster is geboren! Veel gas en stof wordt bij de explosie weggeblazen maar veel zal ook in de nabije omgeving blijven. Dit zal voornamelijk gaan ronddraaien in een platte schijf door behoud van hoekmoment. Dit wordt de protoplanetaire schijf wordt genoemd. Deze zijn ook obervationeel waargenomen, zie hiervoor figuur1. Vanaf hier begint de vorming van een zonnestelsel.

(3)

Figure 1: Dit is een daadwerkelijke foto gemaakt door de ALMA telescoop van een pro-toplanetaire schijf. De donkere ringen zijn mogelijk regio’s waar jonge planeten al stof opgeveegd hebben.

Er zijn verschillende scenario’s denkbaar om te verklaren hoe planeten hieruit vormen. Een mogelijk model is zwaartekrachtsinstabiliteit. Een moeilijk woord voor hetzelfde als wat er met de eerder besproken gaswolk is gebeurd. Doordat er in de schijf dichte klonten van gas ontstaan die onder zwaartekracht instorten, worden planeten gevormd. Vooral de vorming van gasreuzen zoals Jupiter en Sat-urnus kan hiermee goed verklaard worden. De theorie wordt echter niet onderste-und door waarnemingen omdat de chemische samenstelling van deze gasreuzen niet overeenkomt met de verwachtingen.

Een ander model dat hetzelfde proces probeert te beschrijven heet kern-opveging. In dit model gaat men ervan uit dat hele kleine stofdeeltjes aan elkaar plakken door krachtinteracties, een beetje wat er gebeurt met stof als je je kamer niet regelmatig stofzuigt! Deze zullen dan weer met grotere deeltjes aan elkaar plakken, enzovoort. Uiteindelijk gaat zwaartekracht het overnemen als de massa groter wordt, en zullen de grotere jongens snel groeien door botsingen met elkaar, zie figuur2. Een mogelijk probleem met dit model is dat het groeiproces te langzaam gaat om groot genoeg te worden binnen de levensduur van de protoplanetaire schijf.

(4)

Figure 2: Een impressie van het kern-opveging model waarbij objecten blijven botsen tot de schijf is uitgedund en een paar grote jongens overblijven.

Een alternatief scenario genoemd kiezel-opveging gaat uit van een object die een kilometer tot duizend kilometer groot is en in de schijf helemaal alleen groter probeert te worden door het opvegen van milimeter tot decimeter grote deeltjes. In dit thesis onderzoek ik hoe efficient de groei van deze zogeheten planetesimalen is. Planetesimalen zijn eigenlijke hele kleine planeten. Het constante gevecht tussen de zwaartekrachtsaantrekking van de planetesimalen en de wrijvingskracht van het gas waar de planetesimaal doorheen gaat bepaalt de groeitijdschaal. Met andere woorden, worden de kiezels meegenomen met het gas om de planetesimaal heen, of weet de planetesimaal de kiezels uit het gas te trekken door de zwaartekracht? Lees vooral mijn thesis verder voor het antwoord!

(5)

Acknowledgements

Ten years ago I was a problem child in youth detention without a goal and without the will to accomplish anything in life. Now I am contributing to science in the

surrounding environment of my sweet girlfriend, loving friends and the kind, helping people and endless resources at the university of Amsterdam. I want to

thank everybody who have made it possible to be the person who I am today despite the person that I appeared yesterday. I want to thank my parents and my sister for standing by me every step of the way even after the hell I made them go through in my early years. I want to thank my loving friends Maikel, Dennis and Bart for being there when I needed them. Also, I want to thank Sjoerd, Eric and Samantha for the moments of unique insights and friendship along the way of my bachelor.I of course want to thank Chris Ormel for being my supervisor. You were

always available for discussions and help. Your quick reasoning and remarkable insights kept me sharp throughout the whole project and I learned a lot from you

the past couple of months. It was truly a privilege having you as a supervisor! I want to thank my sweetheart Tess for standing by me every step of the way, you

(6)

Abstract

One suggested method to build planets is through pebble accretion, where a plan-etesimal sweeps up dust particles (pebbles) in the protoplanetary disk surrounding the young star. In this thesis I study the growth stage of kilometer to thousand kilometer sized planetesimals by numerically integrating the equation of motion for variable size pebbles. These millimeter to centimeter size particles experience gas drag and interact gravitationally with the planetesimal, which follows an unper-turbed Keplerian orbit at 1 AU. Following a theoretical introduction to the disk structure and particle dynamics, the flow pattern is solved around the planetesimal using both viscid and inviscid flow solutions. Next, the equation of motion of the pebble is introduced theoretically including drag and the two body force. The equa-tion is integrated numerically using Runge Kutta Fehlberg variable step techniques. The results are obtained by obtaining the collisional cross-section and converting this to a growth time scale. I obtained a smooth finite transition between the flow dominated regime (Guillot 2014)[14] and the settling regime (Ormel and Klahr 2010) [7]. From this transition I conclude that a significant barrier is preventing growth of planetesimals sweeping up particles smaller than 1 centimeter. The minimum of this barrier is identified to be around 100 kilometer planetesimal size. For pebble sizes above 1 centimer, gravity takes over and planet formation is possible within the lifetime of the protoplanetary disk. The lifetime of the disk is taken to be 6 Myr.

(7)

(8)

1 Introduction

A natural curiosity has always triggered mankind to study the questions: where do we come from? Are we alone? What characterizes nature’s way of doing what she does so beautifully but destructively, elegantly but complexly? These are bigger questions following from sub questions that have to be answered first before we can dig deeper. If we cannot answer these fundamental sub questions, nature will main-tain a mystery to us. Understanding planet formation is one of those intermediate steps and is still not understood. We have to admit that it is a complex process with a variety of physical processes.

What we do however understand, is the bigger picture from the physics we know and observations. A cloud of gas becomes too heavy for its own gravity, collapses and forms a star. Due to conservation of angular momentum, gas and dust start to rotate around the star in a flat disk called the protoplanetary disk. Dynamics of the gas and dust in the protoplanetary disk need to be studied to obtain knowledge about planet formation. Stuart Weidenschilling, a pioneer in the study of planet formation, did this already in 1977 in his paper about dynamics in the protoplanetary disk and people are still considering this as the building blocks of their research.

1.1 Current planet formation models

Today there are two main models that try to explain planet formation. On one side we have the disk instability model or top-bottom model, on the other side we have core accretion.

In the first model we think that gas settles to the midplane of the disk because the disk wants to be flattened out due to angular momentum conservation and forms clumps. These clumps will eventually get so heavy that gravity makes them collapse and form planets. One problem with this model is inconsistency in expectation and observation of the metallicity ratio of gas giants (Matsuo et. al 2007)[5]

For the latter we expect small pebbles in the gas of the disk to grow by sticking together, forming bigger pebbles that stick together. Eventually you have one big chaotic disk where big gravitational bodies collide with each other and planets form. Bodies starting around 1 km until around 1000km are called planetesimals. One major problem with this model is that the efficiency in outer parts of the disk is too low to form planets within the life time of the disk and collisions for meter sized objects are not sticking anymore but whether destructive. This stagnates the growth process.(Piso and Youdin 2014)[1] This is why planetary formation is still such a hot topic nowadays and scientists all over the world are trying to come up with one consistent scenario.

(10)

1.2 This thesis

In this thesis I investigate a mechanism known as pebble accretion (M. Lambrechts and A. Johansen 2012) [12]. In this model it is assumed that planetesimals already formed in the disk and that it is accreting millimeter to decimeter sized test particles. I numerically integrate the equation of motion of a test particle including the drag force of the gas and the gravitational interactions of the planetesimal and star. This problem is also known as the restricted circular three body problem. The numerical integration is needed for a selection of test particles. From the numerical implementation a growth time scale is deduced for the planetesimal and I analyse how efficient this body can grow through the pebble accretion mechanism.

To obtain the numerical results I go through a number of intermediate steps. In the first section I briefly elaborate on the star formation process. In section two I introduce the protoplanetary disk model. In the third section I talk about the dynamics of pebbles in the disk. In section 4 I solve for analytical solutions for the gas around a spherical object, representing gas-particle interactions. In section 5 the three body problem of the pebble, the star and the planetesimal is solved including gas and gravity. I solve the equation of motion of one test particle numerically and extend this to a selection of non-interacting particles also, the simulation is tested thoroughly in this section for the specified initial conditions and parameters and I quantify my accretion model. Furthermore, the numerical results are discussed and the growth efficiency is determined in the last sections.

2 A brief overview of star formation

When a cloud of gas and dust in the interstellar medium becomes too heavy, gravity exceeds the outwards pressure. The critical mass for which this happens is called the Jeans mass. Assuming a spherical cloud, this can be written as (de Pater and Lissauer 2006)[13] : MJ ≈ 5k_bT G ¯m r 3 4πρ (1)

With ¯m the mean molecular weight, kb the Boltzmann constant and ρ the density

of the cloud. If the cloud is only supported outwards by pressure it will clump and eventually collapse when the mass exceeds this value. Depending on the mass and density of the cloud, it will collapse in a galaxy or a star. If I ignore the pressure term outwards, the collapse happens on the free fall time scale defined by (de Pater and Lissauer 2006)[13]: tf f = 3π 32Gρ 1/2 (2) A protostar forms as the cloud collapses further and eventually when temperatures become high enough, fusion starts and hydrostatic equilibrium is established again. Because the angular momentum initially contained in the cloud is too high for the newly born star, additional material will start rotating in a flattened disk around the star due to conservation of angular momentum. As I discussed earlier in the introduction, it is believed that planets form here. In figure3a schematic overview of

(11)

Figure 3: A schematic overview of how we think solar systems are being formed. Picture obtained from http://thehigherlearning.com

the discussed formation process is shown. From this point the process of planetary formation begins and in the next section I discuss the disk model used for my simulations.

3 The protoplanetary disk

3.1 A minimum mass solar nebula

For the protoplanetary disk model surrounding the young star, the minimum mass solar nebula model (MMSN) is used (Weidenschilling 1977b)[16]. In this model it is assumed that the mass of our eight solar system planets has been spread throughout their orbital planes, see figure 4. Although the protoplanetary disk is much more complex in real life, the proposed model gives a valid first order approximation. Some basic properties of the MMSN model need to be derived to proceed to the next section. The surface density profile follows directly from the MMSN model and at 1 AU this can be expressed as:

Σ(r) = Σ1

r 1AU

−3₂

(12)

Figure 4: An estimation of the surface density of the rocky planets and gas giants spread through their orbital zones starting from Mercury (left) and ending at Neptune (right). The point with the large error bars corresponds with the asteroid belt between Mars and Jupiter. Obtained from Weidenschilling (1977b)[16]

With r the distance from the solar center and Σ1 = 1700g cm2 (See figure 4).

Because the disk is optically thin for radiation from the star the temperature profile can be derived from Stefan-Boltzmann’s law:

T (r) = T1

r

1 AU −1₂

(4) I will assume that the temperature is constant in the vertical direction that is, T (r, z) = T (r). The density profile for the disk can be obtained by first considering the vertical direction of the disk. The gas needs to be balanced by gravity:

Fg =

1 ρ

dP

dz (5)

Because the radial distance is much larger than the vertical distance, r2 _z2_{, the}

vertical component of the stellar gravity can be simplified to:

Fg ≈ −Ω2z (6)

With Ω is the Kepler frequency and z the vertical axis. Using the constant temper-ature in the vertical direction, I can apply the isothermal equation of state for the gas pressure

(13)

Where the sound speed is given by: cs = r kbT ¯ m (8)

With kb the Boltzmann constant, ¯m the mean molecular weight and T the

temper-ature. This gives the differential equation: − Ω2_{z =} c2s

ρ dρ

dz (9)

The solution of this equation is a Gaussian: ρ(r, z) = ρ0(r)e− 1 2( z H) 2 (10) Where the scale height H is defined as:

H = cs

Ω (11)

I can integrate the density profile over the vertical distance to obtain the total surface density and to gain an expression for ρ0:

Σ(r) = ρ0(r) Z ∞ −∞ e−12( z H) 2 dz (12)

Which gives that:

ρ0 =

Σ(r)

H√2π (13)

This brings us to the complete solution for the density profile: ρ(r, z) = Σ(r) H√2πe −1 2( z H) 2 (14) To asses if the self-gravity of the disk can be neglected, I consider the Toomre stability criterion which says that disks with Q ≤ 1 are not stable. This is the case for instance in cool heavy disks. The Toomre parameter Q is given by:

Q = csΩ

πGΣ (15)

It can be seen from this equation that it compares centrifugal and pressure terms of the force balance with the gravitational term. Filling in the conditions at 1 AU stated in AppendixC, I obtain a value of: Q ≈ 55. The assumption that self-gravity can be neglected is therefore justified.

To determine the shape of the disk I consider the aspect ratio. This can be calculated by taking the ratio of the scaleheight given in equation (11) and the radial distance from the center. Taking into account also the r dependency of the scale height, I obtain:

H r ∝ r

1

4 (16)

Which implies that the disk is flared: the scale height is increasing with r. Figure5 gives a cartoon of how the disk structure looks like.

(14)

Figure 5: A rough sketch of the flared disk structure. The physics in this thesis mostly happens in the midplane of the disk.

3.2 Typical lifetime of the disk

To quantify the growth efficiency of the pebble accreting planetesimals, I need a disk lifetime to compare with. The evolution of a protoplanetary disk depends on many dynamic variables but in general there are two mechanisms responsible for depletion of gas from the disk:

1. Dissipation

2. Photo-evaporation

In case of dissipation, gas is accreted onto the star until only debris is left. For the latter, the gas is heated up significantly (thousands of Kelvin) because of exposure to UV-radiation from the star and is transferred away from the disk. In general it follows from observations and theory that the lifetime for a protoplanetary disk is in the range of (Ribas et. al 2015)[19]:

1 Myr ≤ tdisk ≤ 10 Myr (17)

For the minimum solar nebula I assume a typical lifetime somewhere in between (Lada and Kylafis 1999) [2]:

tM M SN ≈ 6 Myr (18)

3.3 Dynamics of the gas

When travelling with solid particles orbiting around a star in a Keplerian circular orbit, the force balance is given simply by:

Fgrav = Fcent (19)

Where Fgrav is the gravity exerted on the particle by the star and Fcent is the

(15)

also have to consider an additional pressure gradient: ∆g = 1

ρ dP

dr (20)

With ρ the gas density. This gives the total force balance for an ideal gas:

Fgrav = Fcent+ ∆g (21)

∆g can be expressed as a fraction of the centrifugal term: ∆g = −2ηv

2 k

r (22)

With vk the Keplerian velocity and r the distance from the star center. I used

a minus sign to make the fractional parameter η positive. The factor 2 becomes convenient for the Taylor approximation of the velocity. Writing out the radial force balance of the rotating gas gives :

v2 g r = v2_k r − 2ηv_k2 r (23)

With vg the velocity of the gas. This can be written as:

vg = vk

p

1 − 2η ≈ vg = vk− ηvk (24)

Where I used that η 1 to expand η to first order. Now, all I have to do is solve for η. For this I use the expression for the pressure of the gas given in equation 7. This can be written in the form:

P = P0

r 1 AU

−n

(25) With r the distance from the solar center. Differentiating gives the expression for the pressure gradient:

1 ρ dP dr = n c2 s r (26)

Tracing back the r dependency of P from equation 3 and equation 8, I obtain the pressure exponent:

n = 13

4 (27)

Equating expression22and 26gives the dimensionless η: η = nc 2 s 2v2 k (28) Inserting the constants gives that:

η ≈ 1.8 · 10−3 (29)

The total deviation from the Kepler velocity is given by:

ηvk ≡ vhw = 5.7 × 103cm s−1 (30)

This result brings us to the particle dynamics which will be discussed in the next section.

(16)

10-2 ₁₀-1 ₁₀0 ₁₀1 ₁₀2 ₁₀3 ₁₀4 ₁₀5 radius [cm] 101 102 103 104 105 106 107 108 109 1010 1011 1012 Stopping time [s]

Stopping time vs radius particle

Stokes regime

Epstein regime

1 year

Figure 6: The stopping time of particles in the gas ranging from particles to planetesimals in the three different regimes. Small particles live in the Epstein regime which slow down very fast, centimeter particles to kilometer sized planetesimals live in the Stokes regime and become more and more inert as is to be expected.

4 Particle dynamics in the disk

4.1 Stopping time and gas drag

Since the gas is rotating slower than the Kepler speed, a particle that is rotating with the Kepler speed in the gas will experience a drag. The stopping time, the time needed to slow down a pebble of mass m moving in the gas with velocity v with a factor e is given by (Weidenschilling 1977a)[15]

tst =

mv FD

(31) With FD the drag force exerted by the gas in the Epstein regime e.a. s < 9₄lmf p

(Weidenschilling 1977A)[15]:

FD =

4π 3 s

2_{ρv ¯}_ν ₍₃₂₎

Where lmf p is the mean free path of the molecules in the gas, ¯ν the mean molecular

velocity and s the radius of the particle. In this regime the drag exerted on the particle is the sum of collisions with individual molecules. In the Stokes regime, which is valid if s > 9₄lmf p, the gas should be considered as a flow because the

collection of molecules behaves as a fluid. The drag law becomes: FD = CDπs2ρ

v2

(17)

With CD a drag coefficient for a spherical object (Whipple 1972) [6] categorized in

three regimes (Weidenschilling 1977a)[15]:

CD =      24Re−1 for Re < 1 24Re−0.6 for 1 < Re < 800 0.44 for Re > 800 (34)

Where the Reynolds number is a dimensionless parameter comparing the viscous forces to the inertial forces in particle-gas interactions:

Re = 2sρv

ξ (35)

With ρ is the gas density and ξ the viscosity of the gas. This quantity will be discussed later in this section. For the region Re > 54 the stopping time is given by:

tst =

6ρss

ρv (36)

Which we refer to as the quadratic regime because the drag force in equation 33 goes quadratically with v. To summarize, the small particles are slowed down very easily. The planetesimals however are practically unperturbed by the gas and start around 105 _{cm. Thus I assume that in the relevant planetesimal range, the bodies}

are following a circular Keplerian orbit at 1 AU. The stopping time for the three different regimes as a function of particle radius is graphically represented in figure 6.

4.2 Particle speed relative to the gas

If a test particle follows an unperturbed circular orbit, the radial component of the velocity is zero. However, since I derived earlier that the gas rotates slower and the particle is moving in the gas, the particle looses angular momentum and compensates this by moving closer to the star. This spiraling inwards is called radial drift (Weidenschilling 1977A)[15]. The radial component of the particle velocity is not zero anymore and I have to solve for a azimuthal component as well as a radial component. The total radial force balance of the particle relative to the gas is:

u tst +(vg+ w) 2 r − v_k2 r = 0 (37)

Where u and w are the radial and azimuthal velocity components of the particle, respectively. Filling in the expression for vg and using that _vw

k 1 combined with

ηvk being small, gives an expression linear in w and u:

ru tst

+ 2vkw − 2ηvk = 0 (38)

To solve for u and w I need another equation. The torque N on the particle is related to the specific angular momentum loss h by:

N = dh

(18)

Figure 7: A cartoon of the defined speeds. Left: the inertial frame. Right: the Keplerian frame. The particle is spiralling inwards since u, the radial component is not zero anymore.

Where h is given by:

h = r(vg+ w) (40)

With w the relative transverse speed of the particle. Filling this expression in for h and using the product rule and the chain rule I obtain:

dh dt = dr dt(vg+ w) + r ∂ ∂r(vg + w) dr dt (41)

To simplify this expression I use that vg is given by equation 24 and r∂v_∂rk = −vk₂.

Neglecting small terms of the order ηvk, u and w, I obtain the expression:

rw tst

− uvk

2 = 0 (42)

Finally, I Solve for u and w, the radial and azimuthal speed of the particle relative to the gas respectively (figure 7):

u = 2ηvk τ 1 + τ2 (43) w = ηvk τ2 1 + τ2 (44) Where: τ ≡ tstΩk (45)

Where Ωk is the Kepler frequency. The total speed of the particle relative to the

gas is thus given by:

∆V = √u2_{+ w}2 ₍₄₆₎

In a Keplerian frame this can be explained as follows: The small particles are tightly coupled to the gas and move with the headwind velocity. In accordance with Wei-denschilling 1977A [15] the meter sized particles drift to the star at a very fast rate. When the particles get bigger and more inert the trajectories can be assumed approximately Keplerian, see figure 8. Now I found the speed components of

(19)

peb-10-2 ₁₀-1 ₁₀0 ₁₀1 ₁₀2 ₁₀3 ₁₀4 ₁₀5 Particle radius[cm] 10-5 10-4 10-3 10-2 10-1 100 101 102 103 104

Particle speed [cm/s]

_{Azimuthal Stokes}

Radial Stokes

Azimuthal Epstein

Radial Epstein

ηvk

Figure 8: Speed components in the discussed regimes, dust grains tightly coupled to the gas follow the head wind, for meter sized pebbles a maximum in radial velocity appears (in accordance with Weidenschilling 1977a [15]). For the big objects the speed can be approximated Keplerian.

bles moving in the disk, I can elaborate on the Reynolds number. It can be seen from equation 35 that the Reynolds number is directly dependent on the speed of the pebble. The speed of the pebble is dependent on the stopping time through equation 45. To trace back the s dependency of the Reynolds number I look at the expressions for the velocity components for the pebble given in equation43and 44. If I consider the Epstein regime and Stokes regime given in equation 32 and 33 respectively I obtain:

τ ∝ (

s Epstein

s2 Stokes (47)

From taking limits of τ I deduce that the magnitude of the relative speed is propor-tional to τ according to:

∆V ∝ (

τ0 if τ 1

τ if τ 1 (48)

From the s dependency of ∆V I obtain the s dependency of the Reynolds number in all four regimes, see figure 9:

∆V ∝          s [τ 1, Epstein] s2 _{[τ 1, Stokes]} s0 _{[τ 1, Stokes]} s0 _{[τ 1, Epstein]} ⇒ Re ∝          s2 s3 s s (49)

(20)

10-2 ₁₀-1 ₁₀0 ₁₀1 ₁₀2 ₁₀3 ₁₀4 ₁₀5 Particle radius (cm) 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101 102 103 104 Reynolds number

Stokes regime

Epstein regime

s=9₄ lmfp

Figure 9: The Reynolds number dependent on particle radius in the Stokes regime and in the Epstein regime.

The result determines whether a gas flow past the pebble is viscid (Re < 1) or inviscid (Re > 1). In the viscid (sticky) case, the flow is highly coupled perpendicular to the direction of movement and it will be difficult to deform single stream lines. In the inviscid case, the flow can be deformed easily and turbulence can play a role when passing a spherical object. These are effects that I am neglecting in this thesis. More about this in the next section.

(21)

Figure 10: Geometry of the system. A chosen point a distance r from the origin makes an angle θ with the incoming flow.

5 Gas-planetesimal interactions

5.1 Theoretical background of a flow passing a sphere

5.1.1 The potential solution

For Re 1 the viscosity of the gas is very small compared with the frictional forces and I use the potential solution for the gas flow to model an ideal inviscid flow. This requires solving the continuity equation:

∂ρ

∂t + ∇ · (ρv) = 0 (50)

The flow has azimuthal symmetry because there is no φ dependency so the flow is entirely described by (r, θ), see figure10. For the unperturbed flow, the components are given by:

v∞ = v_r vθ = v∞cos θ −v∞sin θ (51) I assume the flow is incompressible:

∇ · v∞= 0 (52)

It can be shown that the rotation of the field vanishes:

∇ × v∞ = 0 (53)

Which implies that the field is conservative on a simply connected domain. I can therefore write the potential flow as:

(22)

This requires that:

dΦ∞

dr = v∞cos θ (55)

dΦ∞

dθ = −v∞sin θ (56)

Integrating gives the final expression for the unperturbed potential:

Φ∞(r, θ) = v∞r cos θ (57)

Close to the sphere where the flow has velocity v, I use that for a steady flow the density is constant in time so the first term of equation 50 vanishes. Also, from Kelvin’s circulation theorem (Landau and Lifshitz 1959)[9], it follows that the inviscid flow stays incompressible and irrotational. This gives a more convenient equation:

∇ · (ρv) = 0 (58)

I obtain the Laplace equation:

∇2_{Φ = 0} ₍₅₉₎

To solve for the perturbed flow, I need to solve this differential equation using separation of variables:

Φ(θ, r) = f (r)g(θ) (60)

The Laplace equation in spherical coordinates is given by:

∇2_{Φ =}       ∂2_Φ ∂r2 + 2 r ∂Φ ∂r | {z } 1 r2 ∂ ∂r(r2 ∂Φ∂r) + 1 r2_{sin θ} ∂ ∂θ sin θ∂Φ ∂θ + 1 r2_sin2_θ ∂2_Φ ∂φ2       (61)

Filling in the proposed solution gives the following expression: r2 f f 00 + 2rf 0 f = − 1 g sin θ ∂ ∂θ(sin θ ∂g ∂θ) (62)

This requires both the left hand side as well as the right hand side to be a constant because they do not depend on each other. I define this constant to be n(n + 1). Cleaning up gives two separate equations:

− 1 g sin θ ∂ ∂θ(sin θ ∂g ∂θ) = n(n + 1) (63) r2 f f 00 + 2rf 0 f = n(n + 1) (64)

Equation (63) satisfies the Legendre differential equation given by: d

dx[(1 − x

2₎dy

(23)

Figure 11: Left: The potential flow pattern around a normalized sphere obtained by solv-ing Laplace’s equation. This flow is more applicable for high Reynolds number but I am neglecting turbulent wakes. Right: The Stokes flow pattern around a normalized sphere. The solution most applicable to flow with low Reynolds number.

With x = cos θ. The solutions are the Legendre polynomials: Pn(x) = 1 2n_n! dn dxn(x 2_{− 1)}n ₍₆₆₎

I fill in n = 1 for the reason that it has the fit the angular part of the potential solution and obtain:

g(θ) = P1 = cos θ (67)

For (64) I assume a power law solution of the form:

f (r) = A1r + A2r−2 (68)

I apply the boundary condition that Φ should reduce to Φ∞ when r → ∞ and I also

use that the radial velocity should equal zero at the boundary of the sphere to solve for the constants A1 and A2:

A1 = v∞, A2 =

R3v∞

2 (69)

This gives the final solution for the scalar potential to be: Φ = (v∞r +

R3_v ∞

2r2 ) cos θ (70)

The gradient of this field gives the velocity components: v = (1 −R 3 r3)v∞cos θ ˆer− (1 + R 3 2r3)v∞sin θ ˆeθ (71)

This solution is used to describe the non viscous flow of fluids around a spherical object e.a. for high Reynolds number.

(24)

5.1.2 The Stokes solution

For Re 1 I use the Stokes solution to model the gas flow. In this regime the viscosity is not negligible and the Stokes flow is more applicable. The Stokes flow solution is given by (Landau and Lifshitz 1959; originally found by George Gabriel Stokes) [9]: vr = v∞cos θ 1 + R 3 2r3 − 3R 2r (72) vθ = −v∞sin θ 1 − R 3 4r3 − 3R 4r (73) As can be seen from the given expressions, the Stokes flow contains extra terms correcting for the viscosity of the flow. These terms will have a long range contri-bution because they drop with 1_r. This solution is used to describe the very viscous flow of fluids around a spherical object e.a. for low Reynolds number (Beard et al. 1974) [18] and (Sellentin et al. 2013) [8]. The viscosity term in the solution couples streamlines in the horizontal direction which makes the influence of the physical object of the sphere drop down much slower than in the potential case. The flow patterns for both solutions are shown in figure11.

(25)

Figure 12: The global cylindrical coordinate system with from left to right the star, the planetesimal and the test particle. On the right the local Cartesian system I use for the numerical model.

6 Equation of motion of a pebble

6.1 The circular three body problem

Since I am considering the interaction of pebbles and gas with the planetesimal, it is natural to consider a frame which is rotating with the planetesimal with angular frequency Ω0 = Ωpltsml = ΩKepler. A local Cartesian coordinate system with the

planetesimal in the origin is therefore used for the planetesimal-pebble interactions. This is illustrated in figure12. The equation of motion of a pebble in the local frame is given by:

dv

dt = FCoriolis+ Fdrag+ Fpltsml+ Fcentrifugal+ FF (74)

Where I have used the units of acceleration and interaction between the pebbles is neglected. The expression for the Coriolis force is:

FCoriolis = −2Ω0× v (75)

Where v is the velocity of the pebble and Ω0 is the Kepler frequency. I am assuming

Ω0 points in the z-direction and that the motion is restricted to the xy-plane. This

easily leads to the expression:

FCoriolis = 2Ω0

v_y −vx

(76) The gravitational acceleration of the planetesimal is given by:

Fpltsml= −

GMpltsm

(x2_{+ y}2₎32

x (77)

With Mpltsml the mass of the planetesimal, G the gravitational constant and x, y

the horizontal and vertical distance to the pebble respectively. The gravity of the star is given by:

F_F = −GMF

(26)

With M_F the mass of the star and r the distance to the pebble. The centrifugal term following from the rotation is given by:

FCentrifugal= Ω20r (79)

However, since:

|F_F| ∼ |FCentrifugal| (80)

These forces can be combined:

F_F+ FCentrifugal = Ω20− Ω

2_r ₍₈₁₎

Using that the cylindrical coordinate system is related to the Cartesian coordinate system with the expression (figure12):

r = (r0+ x) ex+ yey (82)

This can be written as:

F_F+ FCentrifugal = Ω20 " 1 − r r0 −3# r (83)

From figure12 it also follows that x, y r0 and therefore I obtain:

r = q

(r0+ x)2+ y2 ≈ r0+ x (84)

Finally, the expression can be simplified to:

F_F+ FCentrifugal = 3xΩ20ex ≡ FTidal (85)

The gas drag can be expressed as: Fdrag =

(v − vgas)

tst

(86) With v the velocity of the particle. Note that the velocity of the gas is also influenced by the rotation of the coordinate system. Every object rotating around a fixed axis with angular frequency Ω satisfies:

dr(t)

dt = Ω × r (87)

For the gas velocity in the local frame I subtract the velocity of the rotating frame:

vgasφ= (1 − η) Ωreφ− Ω0reφ (88)

The same geometric argument used for equation (85) leads to the simplified expres-sion: vgas≈ −3 2Ω0x − vhw ey (89)

Where I found the Kepler shear corrected gas velocity. This gives the tools needed to proceed to the numerical implementation.

If I consider only the x-direction of the tidal and gravitational acceleration and I fill in y = 0 in equation (77) I obtain the Hill radius which I will need in section 7:

x = r0 3

s

Mpltsm

3M_F ≡ Rhill (90)

The Hill radius of the planetesimal indicates the sphere of gravitational influence. Within this region objects become gravitationally bound to the object.

(27)

Figure 13: A schematic overview of the numerical method. I integrate the streamlines of particles and look at the first and last hit to determine an impact parameter. The pebbles start far away from the planetesimal to ensure the initial conditions are unperturbed.

6.2 Numerical method

The equation of motion has to be numerically integrated to obtain the pebble trajec-tories in the gas. The used integrator for the problem is the Runge Kutta Fehlberg variable step method. For further details about this method I refer to appendix B. The pebbles start at an initial y-coordinate ys far enough from the planetesimal.

The separation ∆x between two pebble streamlines is given by the resolution: ∆x = Total number of integrations N

Total x range (91)

To determine the amount of particles being swept up by the planetesimal per unit time, it is useful to start with the impact parameter b:

b = x2− x1

2 (92)

Where x2 and x1 are the positions of the first hitting pebble and the last hitting

pebble respectively. The minimum number of pebble hits on the planetesimal is 15 to obtain a reasonable minimum impact parameter. The total flux through this region is just:

Fcoll = Av · ˆn (93)

With A the surface of the region, ˆn a vector perpendicular to this surface and v the velocity of the pebble. In the symmetrical case this reduces to:

(28)

For higher stopping time this leads to trouble because of the radial drift discussed in section 2 and in the gravitational case, higher interactions. Because the distance increases significantly in the x direction as well as in the y direction I need to correct for the Kepler shear. The z direction is not affected by the shear so I only correct for the x direction:

bef f = b

vy

vhw

(95) To quantify this area of influence I look at the collisional fraction. I define this to be the ratio of the collisional cross-section and the geometric cross-section of the planetesimal: fcoll = πb2 ef f πR2 pltsml = b Rpltsml 2 (96) A sketch of the accretion model is given in figure13.

(29)

7 Numerical results

7.1 Parameter space and initial conditions

Before I start with the accretion process, I have to specify the parameters and initial conditions. Many of the assumptions about the protoplanetary disk follow from the MMSN from Weidenschilling introduced in the first section. Furthermore, The planetesimal follows a Keplerian orbit at 1 AU from the solar like star and it has an assumed mass density of 1 g cm3_{. I am assuming it is unperturbed by the}

gas. I run the simulations for a planetesimal range of 1 to 1000 km. The pebbles considered go from 0.01 cm pebbles to pebbles of about 30 cm. The pebbles are assumed not to interact with each other. An overview of the full parameter space is given in appendix C. The results are obtained with the potential flow solution unless indicated otherwise.

7.2 No gravity case

To tackle the more complex problem of gravity, I first looked at the gas drag exerted on the pebbles encountering the planetesimal. It is useful to define the Stokes number at this point as:

St = tst tenc

(97) With tenc =

Rpltsml

vhw .This equation gives information about the coupling strength of

the pebbles to the gas. The two limits of this equation are simple to understand. For low Stokes number the pebbles are highly coupled to the gas simply because tenc < tst. This implies that the time to travel past the planetesimal is shorter

than the time the pebble is being slowed down. For high Stokes number tenc > tst.

In this case the particle is slowed down before full encounter, leading to a higher collision probability. To summarize, a higher stopping time implies a bigger pebble and this in turn implies that it is less perturbed by the gas. If the planetesimal gets physically bigger, the time to go past it gets higher. The obtained results agree with the hydrosimulation obtained by Sekiya and Takeda(2003)[11].

The pebble streamlines for St 1, St = 1 and St 1 are shown in figure14. For St = 0.01 the pebbles are strongly coupled to the gas and go past the planetesimal. For St = 1 the pebbles are affected by the gas but there are some pebbles that are not coupled enough to the gas and collide with the planetesimal. As for for St = 100 stream lines, the pebbles hardly feel the gas anymore and the collisional cross section equals approximately the geometric cross section of the planetesimal. The red trajectories are the hits and the gray trajectories the misses.

7.3 The full model

Although the above model gives rise to some interesting physics which needs to be understood first before I can move on, it is not the complete picture. Equation (74) gives the equation of motion that I have to use to develop the full accretion model for the planetesimal. The two body interaction of the planetesimal increases the influence of the planetesimal as can be seen in figure 15. This has been tested for

(30)

the same parameters as our gas drag study. There still is a coupling to the gas but the gravitational focusing of the two-body force pulls pebbles back to the planetesimal, which increases the region of influence. In this case only the Stokes number is not sufficient because gravity is also dependent on the radius of the planetesimal and the stopping time.

For St = 0.01 there are already numerous collisions due to gravitational influence. For St = 1 the collisional cross-section becomes a multiple of the planetesimal radius since the pebble radius is increased by a factor 100 and gravity becomes dominant. Also, because of the tidal force due to the star, the streamlines are starting to become asymmetrical. For St = 100 the scale of the figure is more global because of high gravitational interactions. Pebbles start far outside from the Hill sphere of the planetesimal (eq. 90) to ensure the pebbles are unperturbed at the starting point.

The collisions that I discussed so far can be categorized in three regimes: 1. The geometric regime;

2. The flow dominated regime (Guillot et. al 2014) [14]; 3. The settling regime (Ormel and Klahr 2010) [7].

The first regime shows accretion in which the impact parameter b is approxi-mately equal to the planetesimal cross section. In the second regime the pebbles hardly fall on the sphere because they are coupled to much to the gas, this is called the flow dominated regime. In the latter regime gravitational focusing in combina-tion with the gas drag increases the impact parameter. The gas drag is responsible for lowering the momentum of the incoming particle. This in turn makes a collision more probable with the planetesimal since energy is not conserved anymore due to friction. If only gravitation was present, conservation of energy would prevent point particles from colliding in close encounters. For gravity dominated interactions when including gas drag, the physical body can be shrunk down to a point particle and gravitational focusing will still be approximately equal. This is therefore also called the point particle regime or settling regime.

(31)

Figure 14: Left: pebble streamlines for St = 0.01. Right: pebble streamlines for St = 1. Below: pebble streamlines for St = 100. The red streamlines are the hits and the gray are the misses. The gas stream lines are not shown but resemble the potential solution given in 11.

(32)

Figure 15: Particle streamlines with the two body force on. Left above: pebble trajectories for St = 0.01. Right above: pebble trajectories for St = 1. Below: pebble trajectories for St = 100. The dashed circle represents the Hill sphere of the planetesimal given by equation 90.

(33)

Figure 16: Collisional fraction numerically obtained for both the Stokes flow solution and the potential flow solution.

7.4 Stokes flow vs. potential flow for 0.01 cm pebbles

For a pebble size of 0.01 cm a comparison has been made between the Stokes flow solution and the Potential flow solution for the collisional fraction fcoll, see figure

16.. The collisional fraction starts at approximately the geometrical cross-section for small planetesimals. Then it goes down significantly in the flow dominated (no gravity) regime discussed earlier until a minimum is reached around Rpltsml ≈

100km. After the minimum a smooth transition occurs to the settling regime in which gravity takes over. For the Stokes model as well as for the potential model there is a barrier. Next to the full model, the flow dominated curve (no gravity) as well as the point particle curve (flow neglected) are shown separately. The velocity of the gas stream lines increases near the sphere boundary while it goes to zero in the Stokes flow case. Because of this, the particles will be moving faster in the potential case close to the sphere, and the collisional probability goes down. This explains that the full model curve is starting below the settling curve for the potential model. When gravity increases, the full model curve merges with the settling curve again because the point particle regime is entered here. Also, as discussed in section 5, the Stokes flow is more viscous which can make the barrier more shallow because inertial forces are low and particle settle more efficiently.

7.5 Collisional fraction for full parameter space

The numerical results for the full parameter space are shown in figure17. These re-sults are obtained with the potential flow solution in accordance with high Reynolds number, as I discussed in section 5. The lowest three curves correspond to the right figure in figure 16. In all cases for 1 km planetesimals the collisional fraction starts at the geometric cross-section (geometric regime). The collisional fraction goes down significantly in the flow dominated regime since particles are coupled too much by the gas. Eventually the 10−2 and 0.2 cm curves have a minimum at

(34)

Figure 17: The collisional fraction as a function of planetesimal radius obtained with the potential flow solution.

around Rpltsml ≈ 100km. This shows that a finite barrier occurs when the particles

being accreted are too small. When the planetesimals become bigger, gravity dom-inates (settling regime) and the barrier disappears. For pebbles larger than 0.2 cm the cross-section exceeds the geometrical cross-section and there is no minimum anymore.

7.6 From collisional fraction to growth time scale

To compare the growth rate of the planetesimal with the lifetime of the disk, I need to derive a growth timescale. In the derivation I assume that collisions with the planetesimal are sticky. In this case, the rate of change of the mass is just what comes through the impact surface (eq: 94) multiplied by the density of pebbles:

˙

M = Fcollρpart (98)

This can be converted to a characteristic growth time: tgrowth =

M ˙

M (99)

Where tgrowth is defined as the time needed for the planetesimal to double its size.

(35)

Figure 18: The growth rate for planetesimals in the given range as a function of planetes-imal size and a pebble density of ρgas

100, following from the dust to gas ratio of the MMSN. The dashed lines represent the constant lines when filling in the corresponding stopping time in equation104.

the encounter time. Where the settling time is given by: tsett =

b vsett

(100) Where vsett= fgtst, this is also known as the terminal velocity. With fg the gravity

exerted on the pebble and b the impact parameter. This gives that: tsett=

b3

tstGMpltsml

(101) The extreme case is that tsett is equal to tenc:

b3

tstGMpltsml

= b

vhw

(102) Where I used that the encounter distance is equal to the impact parameter since this is the distance that the pebble is free from influence of the planetesimal. This can be written as:

b2 = 4tstGMpltsml vhw

(36)

The extra factor 4 was obtained in a more detailed analysis by Ormel and Klahr (2010) [7]. From equation (99) the following expression arises:

tgrowth ≈

1 4πGρparttst

(104) Where ρpart is the number of particles per unit volume. The value follows from the

gas to dust ratio in the MMSN and is taken to be 0.01ρgas. The expression given in

equation 104 is independent of the planetesimal properties and gives the constant the growth time must approach eventually. For test pebbles ranging from 0.01 cm to 30 cm and a planetesimal range of 1 km to 1000 km, tgrowth is obtained from the

numerical simulations. See figure 18. The smaller the test pebble size, the larger the grow barrier becomes. The barrier is identified at around 100 km planetesimal radius. It can be seen that for the largest three pebble sizes, curves are crossing each other. Since the stopping times are increasing, settling becomes harder for these pebbles because they are not being slowed down enough within the encounter time. However, when the gravity of the planetesimal becomes high enough, the process speeds up again and the pebbles settle eventually.

(37)

8 Conclusion and discussion

In this bachelor thesis I have investigated whether planetesimals in the range of 1-1000 km in an unperturbed Keplerian circular orbit at a distance of 1 AU from the star are able to grow by accretion of pebbles in the range of 0.01 cm to 30 cm. I conclude from the growth rate of the planetesimals that there is a growth barrier between the flow dominated regime and the settling regime until the test pebble size reaches around 1 cm. In this region pebble accretion kicks in and growth becomes easier. Using an optimistic typical disk lifetime of tdisk ≈ 6 Myr it appears that full

growth can be realized for Rpebble > 2 cm. The barrier is identified around 100 km

planetesimal size for each pebble size. While this barrier can be a real bottleneck for growth, the barrier is not infinite. A smooth transition has been found between the two regimes and gravity will take over at a certain point. One possible scenario could be that growth stops when the barrier has been reached and that another mechanism takes over from here that ensures tgrowth < tdisk and the model is still

viable.

Another implication of this research is that there are significant differences be-tween the gas flow solution in the Stokes regime, a viscosity dominated flow (Re 1) and the potential flow regime, the non-viscous flow (Re 1). In the first case the barrier is less deep than the latter because more pebbles stick to the planetesimal. This is probably due to the viscosity of the gas which overtakes the inertial forces.

The key assumptions I used in the pebble accretion model presented in this thesis are summarized as follows:

• Interaction between pebbles have been neglected • The planetesimal follows an unperturbed circular orbit • The flow pattern is assumed to be incompressible • Turbulent wakes have been neglected

• Pebble reservoir infinite in the disk

• Collisions with the planetesimal are sticky • Rotation of the planetesimal has been neglected • Migration of the planetesimal has been neglected

Although the model is highly simplified, the physics contained in it still gives a valid approximation for the global dynamics of solids and fluids in the disk. Interac-tions among pebbles are not taken into account and can alter the accretion process. Inside the Hill sphere however, I expect no significant differences as collisions lower momentum which may even increase the accretion of pebbles. For a more detailed analysis of the collisional physics of dust I refer to Ormel, Spaans and Tielens (2007) [3] and Krijt (2015) [17].

(38)

I also assumed that the gas flow pattern in the disk is incompressible. If the gas becomes compressible an atmosphere will form around the planetesimal which influences the accretion proces. Detailed verification of the incompressibility of the gas can be found in appendix A.

The assumption that the planetesimal follows an unperturbed orbit is justified for Rpltsml> 1 km because the decay rate is approaching the lifetime of the disk.

The effects of a turbulent wake behind the planetesimal due to low viscosity have also been neglected. This could alter the accretion process a lot since the analytical solution takes particles far from the sphere again, while in the turbulent case this could be different. Future research is needed to study these effects.

Furthermore, the amount of pebbles to sweep up in the disk is assumed infinite. Eventually however, there will be no pebbles left and the accretion process stops. This could be calculated more specifically for a finite dust to gas ratio.

The collision with the planetesimal is assumed sticky. This requires further re-search because small dust grains could bounce back elastically and go along with the gas stream pattern again. Research has been done on this subject at 1 AU. For the very small particles the collisions can be assumed sticky. For the bigger pebbles however, the collisions enter the red zone and fragmentation occurs (Testi et al. 2014) [10].

Including rotation of the planetesimal goes beyond the scope of this work. In this case the gas will spin up around the planetesimal and the rotation of the gas cannot be neglected anymore.

Future research building on the resources created within this bachelor project could be possible atmosphere creation by modeling a compressible gas. Another scientifically relevant aspect to this problem is rotation build up by the planetesimal. Varying parameters as gas density and orbital distance from the star are also an option.

(39)

APPENDICES

A

Is an incompressible flow justified?

In solving for the gas flow past the planetesimal I assume that the gas is incompress-ible. However at the point where the escape velocity of the planetesimal becomes bigger than the speed of sound in the gas, a significant atmosphere is being formed around it (A. A. Piso and A.N. Youdin 2014)[1]. Equating the speed of sound to the escape velocity of the planetesimal gives the Bondi radius:

Rb ≡

GM c2

s

(105) In the case that Rpltsml > Rb, the gas can be assumed incompressible. Since

the tidal interactions from the sun do not play a big role yet, there is spherical symmetry and hydrostatic balance between the gas and the planetesimal. Solving for the radius of the planetesimal equalling the Bondi radius as a function of disk radius, I obtain: R(r) = s 3kbT1 4πGµmHρ 1AU r 1₄ (106) If the planetesimal radius exceeds the Bondi radius, the gas can be considered incompressible. If the planetesimal radius is within the Bondi radius, the gas is compressible and the solution is not valid anymore, see figure19. In the planetesimal range I am considering, it is safe to assume the assumptions made are valid.

Figure 19: A logarithmic representation of the Bondi radius equalling the radius of the planetesimal as a function of disk radius. Below the line the gas is incompressible and above it, the flow is compressible.

(40)

B

The Runge Kutta Fehlberg variable step method

To integrate the equation of motion, I wrote an integrator based on the Runge Kutta Fehlberg variable step method (RKF45). This method is considered especially efficient and accurate for problems in which energy is not conserved. Since the velocity dependent gas drag is in our equation of motion, this is a suited integrator. Newton’s equation can be approached mathematically as an initial value problem. For the position dependent accelerations this can be written as:

¨

x = a(x, t) (107)

For the drag force I need to consider the velocity dependency and write it as: ˙

v = a(v, t) (108)

Numerically this differential can be written in array form as: ˙

Y = [ax, ay, vx, vy] (109)

With a the acceleration and v the speed. Given an array of initial conditions for the speed and position, which I will call Y0:

Y0 = [vx0, vy 0, x0, y0] (110)

These equations can be solved by iteration. For the iteration, Runge Kutta order four and five respectively are compared to determine the next step size h:

(

Y_{N +1}order 4 = Yn+ (a18k1 + a19k3+ a20k4− a21k5)

Y_{N +1}order 5 = Yn+ (a22k1 + a23k3+ a24k4− a25k5+ a26k6)

(111) The values of anand the prescription of kncan be found in a paper of Es-hagh (2005)

[4]. The difference between the two solutions gives an array of absolute errors: E = |Y_{N +1}order 4− Yorder 5

N +1 | (112)

If the error between the two approximations exceeds a user specified tolerance which I will call toluser, the step size is reduced with a fudge factor:

s = 0.84 Min [f E]

1

4 (113)

Where f is an array with components:

f = toluser[|v|, |v|, |x|, |x|] (114)

Used to convert toluser to a minimum absolute tolerance. Now I can relate this to

the error of the two solutions to obtain a relative error. A new step size will be used in the iteration provided that the error is smaller than the specified tolerance:

Max E f

< 1 (115)

If not the case, the step size will continue to be reduced until it becomes smaller than 1 and only then the iteration will continue to the next step.This provides an efficient but still accurate way to integrate the equation of motion given in eq (74).

(41)

C

Numerical values of used parameters

Table 1: Numerical values of parameters used in calculations

Parameter Definition Value

vhw Headwind velocity 5700 cm s−1

Rpltsml Planetesimal radius 1 - 1000 km

ρgas Gas density 10−9 g cm3

tst Pebble stopping time 102− 106 sec

Rpebble Radius of the pebbles 0.01 − 30 cm

ρparticle Amount of particles per volume ρ₁₀₀gas

ρpltsml Planetesimal density 1 g cm−3

r0 Orbital distance planetesimal 1 AU

Σ Gas surface density 1700 g cm2

vk Kepler speed 30 km s−1

¯

m Mean molecular weight 2.2 × 10−24g

cs Sound speed 105 cm s−1

ξ Gas viscosity 1.7 × 10−4cm2 s−1

¯

ν Thermal velocity molecules 1.87 × 104 cm s−1

T Temperature 300 K

lmf p Mean free path molecules 1.84 cm

(42)

References

[1] A.N.Youdin A.A.Piso. On the minimum core mass for giant planet formation at wide seperations. arXiv:1311.0011v2 [astro-ph.EP], ”2014”.

[2] C.J.Lada and N.D.Kylafis. The origin of stars and planetary systems. Springer, 1999.

[3] A.G.G.M.Tielens C.W.Ormel, M.Spaans. Dust coagulation in protoplanetary disks: porosity matters. Astronomy and Astrophysics, 461:215, 2007.

[4] M. Es-hagh. Step variable numerical orbit integration of a low earth orbiting satellite. Journal of the Earth Space Physics, 31:1–12, 2005.

[5] T.Matsuo et. al. Planetary formation scenarios revistied: Core-accretion versus disk instability. arXiv:astro-ph/0703237v1, 2007.

[6] F.L.Whipple. From plasma to planet. ed.A.Elvius,Wiley,London, 211, 1972. [7] C.W.Ormel H.H.Klahr. The effect of gas drag on the growth of protoplanets.

Astronomy and Astrophysics, 520:A43, 2010.

[8] C.P.Dullemond E.Sellentin J.P. Ramsey, F.Windmark. A quantification of hy-drodynamical effects on protoplanetary dust growth. Astronomy and Astro-physics, 560:A96, 2013.

[9] L.D. Landau E.M. Lifshitz. Fluid Mechanics ( Volume 6 of A Course of Theoretical Physics ). Pergamon Press, 1959.

[10] L.Ricci S.Andrews4 J.Blum J.Carpenter C.Dominik A.Isella A.Natta J.P.Williams D.J.Wilner L.Testim, T.Birnstiel. Dust evolution in protoplan-etary disks. arXiv:1402.1354v1 [astro-ph.SR], ”2014”.

[11] H.Takeda M. Sekiya. Were planetesimals formed by dust accretion in the solar nebula? Earth Planets Space, 55:263269, 2003.

[12] M.Lambrechts and A.Johansen. Rapid growth of gas-giant cores by pebble accretion. Astronomy and astrophysics, arXiv:1205.3030v2 [astro-ph.EP], 2012. [13] I.de Pater and J.J.Lissauer. Planetary sciences. Cambridge University Press,

2006.

[14] T. Guillot S. Ida, C.W.Ormel. On the filtering and processing of dust by planetesimals. Astronomy and astrophysics, 572:A72, 2014.

[15] S.J.Weidenschilling. Aerodynamics of solid bodies in the solar nebula. Mon. Not. Roy. Astron. Soc., 180:57–70, 1977A.

[16] S.J.Weidenschilling. The distribution of mass in the planetary system and the solar nebula. Astrophys. Space Sci., 51:153–158, 1977B.

(43)

[17] S.Krijt. From grains to planetesimals: the microphysics of dust coagulation. Ipskamp drukkers, 2015.

[18] K.V.Beard S.N. Grover. Numerical collision efficiencies for small raindrops colliding with micron size particles. Journal of the atmospheric sciences, 31:543– 550, 1974.

[19] B.Mern . Ribas, H.Bouy. Protoplanetary disk lifetimes vs stellar mass and possible implications for giant planet populations. ArXiv eprint, 2015.

On the growth efficiency of pebble accreting planetesimals