Pebble accretion
An analysis of the behaviour of non-accreting pebbles Niels van Hese
11646462 July 21, 2020
Report Bachelor Project Physics & Astronomy
15 EC
Conducted between 30-3-2020 & 21-7-2020
woo
Supervisors: Prof. Dr. Carsten Dominik, Dhr. Rico Visser Msc. Second Examiner: Dr. Jean-Michel Desert
w w w w w w w w w w w w w w
Anton Pannekoek Institute for Astronomy
Abstract
Pebble accretion is a promising new avenue for the formation of planets within observed disk lifetimes. Prior research has mostly covered pebbles directly accret-ing onto the planetesimal without consideraccret-ing collisions. In this thesis a numerical model for pebble accretion is constructed. With it, the effects of adding collisions, in particular for pebbles that do not directly accrete, on mass accretion and spin accumulation of planetesimals are investigated. The regimes of ballistic accretion and the onset of pebble accretion are considered. For the ballistic regime, the data implies that collisions have a significant effect. For a planetary radius of 100 km, adding collisions was found to increases the prograde spin gained via pebbles by a factor ten. For the accretion regime the data suggests that collisions have little to no impact on either mass accretion or spin accumulation. Suggestions for further research along the lines of this thesis are presented.
Populaire samenvatting
Wanneer een ster vormt, laat het een grote schijf met stof, gas en deeltjes achter. Hoe hier planeten uit ontstaan is een vraag waar nog altijd aan gewerkt wordt. Een van de nieuwste idee¨en is ’pebble accretion’, een model waarin kleine steentjes van een paar millimeter tot een paar decimeter groot ver in de buitenkant van de schijf vormen. Daar worden ze afgeremd door het aanwezige gas, en beginnen ze inwaards richting de ster te vallen. Als ze dan onderweg een planeet, of iets wat een planeet kan worden, tegenkomen, worden ze door diens zwaartekracht afgebogen, soms zo ver dat ze op de planeet terecht komen, zoals te zien in het plaatje hieronder. Daar is een dergelijke planeet in het midden te zien, met de paden van een aantal van die steentjes eromheen aangegeven, met in het rood degene die op de planeet eindigen, en in het grijs degene die afgebogen worden. Zoals ook te zien, kunnen deeltjes die afgebogen worden nog bij elkaar terecht komen en botsen. Eerder onderzoek heeft botsingen genegeerd, maar ik heb onderzocht of deze wel meenemen een effect heeft op de hoeveelheid van die steentjes die er op de planeet eindigen.
Acknowledgements
I would like to specifically thank a number of people very close to me for helping me through this project.
First of all, thank you Jurriaan and Steven! Without your expertise at C and keen insight, Rebound’s source code might well have remained a mystery to me. Thank you too for your invaluable assistance in getting the simulations running and gathering the data. You have saved me many a trip to stack exchange!
Second of all, thank you Annika! Without your eye for stylistic devices and experience at writing papers, this thesis would have looked much worse, and much less consistent. But thank you especially for the emotional support you gave me throughout this ordeal, without which I could have seen myself throw in the towel long ago.
Third and last, but certainly not least, thank you Maaike! Without your linguis-tic abilities my thesis would, without a shadow of a doubt, have remained riddled with slight mistakes and unsuitably informal language.
Contents
1 Introduction 1 2 The disk 2 2.1 Overview . . . 2 2.2 Underlying assumptions . . . 2 2.3 Gas dynamics . . . 4 2.4 Particle properties . . . 5 2.5 Particle movement . . . 6 2.6 Particle formation . . . 83 Three body system 10 3.1 Co-rotating frame . . . 11 3.1.1 Shear . . . 13 4 Numerical methods 13 4.1 Rebound . . . 14 4.2 Pebble injection . . . 15 4.2.1 Pebble flux . . . 15 4.2.2 Injection . . . 16 4.3 Numerical values . . . 17 5 Results 17 5.1 Model testing . . . 17 5.2 Initial conditions . . . 19 5.3 Ballistic regime . . . 21
5.4 Onset pebble accretion regime . . . 21
6 Discussion 24
7 Conclusion 25
8 Bibliography 26
A Weighted randomizer 28
1
Introduction
Planetary formation is a busy and ongoing field of study. There are a myriad of questions left unanswered as of yet. One of the most interesting of these, is the question of how the solar system formed. This question has not been answered to a satisfactory degree, as previous core accretion models have proven to be unlikely to be able to produce planets of the size our solar system contains over reasonable time scales (Haisch & Lada, 2001 [1]).
Thus, it might be time to look for a new model. In this thesis the pebble accretion model for planetary formation (Ormel and Klahr, 2010 [2], Lambrechts and Johansen, 2012 [3]) will be tested in a specific situation. The pebble accretion model assumes that there already exist some larger planetesimals (relatively small planetary bodies) within a protoplanetary disk, and that small particles (pebbles), between milli- and decimeter in radius, form in the outer reaches of that disk, driven inward by drag forces from the gas present. As they are driven inward, these peb-bles may encounter the planetesimals, be caught in their gravitational wells and subsequently accreted.
The pebble accretion model has been shown to explain the rapid formation of planetary cores, and can be used in conjunction with the other models to more accurately predict and/or retrace the formation of planetary systems much better than those other models could by themselves. Most research concerned with this new model has been focused on the accretion of the pebbles on planetesimals.
In this thesis however, I will be looking not at the pebbles that actually accrete directly onto the already formed planetesimals, but at the ones that narrowly miss, and potentially collide with pebbles that miss the other side of the planetesimal. The behaviour and potential contributions of these particles to the formation of planetary systems is a relatively new area of study, and potentially significant to the evolution of the disk as a whole. What I specifically want to know for this thesis is whether these colliding particles lose enough of their momentum in their collisions to then be accreted onto the planetesimal indirectly, thereby contributing to the overall accretion rate.
But before all that, I will first lay out the system we will be looking at, consisting of the different components like gas, dust and larger particles inside a protoplanetary disk. Then the interactions between these different components will be laid out as they pertain to the system that will be simulated, which will be a small box around a planetesimal as it orbits the central star. That leads naturally into how the system is modelled numerically in my simulation, and from there we can easily start taking a look at the results that the simulations have produced, and drawing conclusions from them.
2
The disk
2.1 Overview
As stars form from interstellar dust clouds collapsing in on themselves because of perturbations from nearby supernovae and likewise impactful events, so too do the protoplanetary disks surrounding said stars form (Prialnik, 2009 [4]). These disks contain all the gas and dust not (yet) accreted onto the central star. The dust, by some means other than pebble accretion, forms a number of planetesimals, and many, many pebbles. A possible avenue for the creation of these planetesimals is the streaming instability, which holds that the gas drag acting on small particles in a disk can cause spontaneous clumping, leading into a gravitational collapse to a planetesimal (Johansen & Jacquet, 2015 [5]).
In this chapter, the physics guiding the movement of and interactions between all these elements of the disk will be discussed. I will start by outlining the necessary assumptions, to then move on to the dynamics of the gas, and finally the effects those have on the movement of particles in the disk.
2.2 Underlying assumptions
The theory used in this thesis could not meaningfully be discussed without first dis-cussing the assumptions used about the system as a whole. Thus, these assumptions will be discussed first here, to then be referenced as they become relevant later on. The system discussed in this thesis will be a protoplanetary disk following the model of the Minimum Mass Solar Nebula (MNSN)(Weidenschilling, 1977b [6]; Hayashi et al., 1985 [7]; Nakagawa et al., 1986 [8]). In this model, a number of key assumptions are made. Most importantly the temperature and surface gas den-sity profiles, are both assumed to follow neat power-law distributions:
T (r) = T0 r au −1/2 , (1) and Σ(r) = Σ0 r au −3/2 , (2)
with r the distance from the central star. Further assumptions include the scale height H of the disk being
H = cs Ωk
= csr vk
, (3)
where cs is the speed of sound within the gas, and Ωkthe keplerian frequency, which
is defined as the keplerian velocity vk over r.
Last but not least, we assume the gas to follow the ideal gas law, meaning that the pressure is given by
P V = nkbT P = ρg ¯ mkbT P = ρgcs2, (4)
where ρg is the volumetric gas density, and ¯m is the mean molecular mass.
With these assumptions the volumetric gas density can be calculated, by inte-grating the surface density over the vertical direction. To do this, we first look at the force balance on a parcel of gas in the z-direction. We assume the disk is in a steady state, thus the forces should add up to zero. This gives
Ftotal= Fgz+ FPtop+ FPbottom
= sin(arctanz r )mg + PtopA − PbottomA Ftotal= 0 ⇒ sin(arctan z r )gρgA · dz + A · dPz= 0, (5)
which, given that z << r, can be simplified to z
rgρgdz + dPz= 0. (6)
Rewriting this and using equation 4 gives
1 ρg dρgcs2 dz = − gz r = − vk2z r2 = −Ωk 2z 1 ρg dρg dz = − z H2. (7)
Then by integrating both sides over dz
Z 1 ρg dρg = − 1 2 z H 2 ρg = Ce− 1 2( z H) 2 , (8)
where C is an integration constant resulting from the indefinite integration. Then by integrating this over all z, the gas surface density should result:
Z ∞ −∞ Ce−12( z H) 2 dz = Σ(r). (9)
By recognizing that this integral is symmetric around 0, and using the standard integralR0∞e−ax2dx = 12pπa, it can be solved:
2C Z ∞ 0 e−12( z H) 2 dz = Σ(r) √ 2πHC = Σ(r) C = √Σ(r) 2πH. (10)
Thus the volumetric gas density is given by ρg =
Σ(r) H√2πe
− z2
Definition Value Units
T0 Temperature in disk at r=au 170 K
Σ0 Gas surface density at r=au 13000 kg · m−2
M? Star mass 1.988 · 1030 kg
G Gravitational constant 6.674 · 10−11 m3 · kg−1 · s−2
kB Boltzmann’s constant 1.381 · 10−23 kg · m−2 · s−2 · K−1
Table 1: This table contains all the (assumed) constants used in this thesis, with T0
and Σ0 taken from Visser et al. 2020 [9].
Using this definition and the other assumptions, the pressure can be rewritten to the following power law:
P = P0
r au
−13/4
. (12)
Naturally, to produce numerical results, we will need to assume physical values for some constants, such as T0 and Σ0. The assumptions for these constants used in
this thesis are all outlined in table 1.
2.3 Gas dynamics
With the assumptions made clear, I can begin to describe the physics we will need to accurately model pebble accretion. To start then, let us look at the gas present in the disk.
Normally, without a gas present, any particle would simply orbit the central star with the keplerian velocity, given by
vk= r GM? r = √ gr, (13)
as it is merely affected by the gravity of the star, and the centripetal force caused by its own motion. A gas particle however, feels an extra force, namely, the pressure of all its fellow gas particles, pushing it along the pressure gradient (Weidenschilling, 1977A [10]). This pressure manifests as a correction to the gravitational force, equaling ∆g = 1 ρg dP dr = − 13 4 cs 2, (14)
and thus resulting in a new balance of forces, and thereby a new velocity too:
vg2 r = g − ∆g vg = vk s 1 −∆g g . (15)
Figure 1: Exaggerated representations of pebble (red) interactions with the gas in the Stokes regime on the left, where pebbles interact viscously with the gas, and the Epstein regime on the right, where pebbles interact on the level of individual gas particles.
This, in the limit where |∆gg | < 1, can be expanded in a Taylor series, and then by neglecting all terms of order higher than 1, be rewritten to:
vg = vk(1 −
∆g
2g) = vk(1 − η), (16)
where η is the fractional deviation of the gas velocity from the keplerian velocity, and ηvkis the actual difference between the gas and keplerian velocities, which will
from this point on be referred to as the headwind velocity.
It is because of this headwind velocity that the behaviour of particles in a proto-planetary disk where gas is present is so different from their behaviour in vacuum. How exactly it differs is what will be looked at next.
2.4 Particle properties
In this section the tools needed to desribe the movement of particles as it is influenced by the presence of the gas will be discussed, to then be put to use in the next section. Now, because of the headwind, any non-gas particles travelling through the gas will experience a drag force. There are two distinct regimes for this drag force (Whipple, 1972 [11]). They are equal to one another when the following condition holds:
s = 4
9λ, (17)
where s is the pebble radius, and λ is the mean free path of a gas molecule. When s is less than this, the pebble is in the Epstein regime. This means a drag force equaling FDrag = 4π 3 ρgs 2v sv,¯ (18)
where vsis the velocity of a particle, and ¯v is the thermal velocity of the gas particles.
When s is higher, the the drag force is equal to Fdrag=
1
2CDπρgs
2v
s2, (19)
where CD is a drag coefficient that depends only on the dimensionless Reynolds
number smaller than 1, which holds in all cases of interest to this thesis, the drag coefficient is then given by:
CD= 24Re−1, (20)
By plugging this back into the equation for the drag force above, we have the drag force in what is known as the Stokes regime, which looks like:
Fdrag= 6ηgπsvs. (21)
These expressions for the drag force can be used to calculate another useful quantity: the stopping time. The stopping time denotes the time it takes for the drag force to reduce the speed of a particle by a factor e, and is related to the drag force as follows:
ts =
mvs
Fdrag
, (22)
with m being the particle’s mass, given by 4π3 ρss3.
When this is worked out, the expressions for the stopping time in the Stokes and Epstein regimes read as
ts ( 2ρss2 9ηg Stokes Regime. ρss ρ¯v Epstein Regime. (23)
Now only one more definition needs to be established, that of the Stokes number. The Stokes number is a dimensionless number, given as
τs =
vkts
r , (24)
and characterizes the behaviour of our particles in a concise manner. Its use in our situation will be apparent in the next section.
2.5 Particle movement
Now that all the tools needed to discuss the behaviour of particles in the disk’s gas have been outlined, the actual discussion can begin.
To start, say a particle has a speed w relative to the gas in the transverse direction (thus a total speed of vg+ w in the transverse direction), and a speed of u in the
radial direction. The total speed relative to the gas then is v =√u2+ w2.
Of course, just naming these speeds is not enough. So with the derivations from the previous section and the help of the particular angular momentum
h = r(vg+ w), (25)
of the particles, we can create a linear system of equations that we can solve to get u and w.
In order to do this, we will have to look at two separate situations. First, with u and w defined, the drag force can be split in the different directions by scaling
it with wv and uv for the transverse and radial directions, respectively. The force balance in the radial direction then becomes
Fr= FDu mv + (vg+ w)2 r − g = 0 = FDu mv + vg2+ 2vgw + w2 r − g = FDu mv + (∆g + g)r + 2(vk− ηvk)w r − g = FDu mv + 2vkw t + ∆g, (26)
by neglecting terms of the order w2 or higher.
The other situation to look at is the change in h over time. This change is defined as the torque on a particle divided by the mass, thus given by
dh dt =
rFDw
mv . (27)
However, the change over time can also be calculated by taking the time derivative of h as follows: dh dt = d(r(vg+ w)) dt = dr dt(vg+ w) + r d(vg+ w) dr dr dt = dr dt[(vg+ w) + r d(vk− ηvk+ w) dr ] = u[(vg+ w) + r dvk dr ] = u[(vg+ w) + 1 2 √ gr] = vku 2 , (28)
which is possible in the limit where w ≈ ηvk.
Now, by equating these two versions of dhdt, we obtain the second part of the system of linear equations:
FDw
mv −
vku
2r = 0. (29)
Now that the system of linear equations is complete, we can get to solving it. Rewrit-ing equation 29 gives
w = vkteu
2r , (30)
which can be plugged into equation 26, thereby giving an equation defining u:
u = − r 2∆gt e r2+ v k2te2 . (31)
Plugging this back into equation 30 leads to the equation defining w: w = − ∆grvkte 2 2(r2+ v k2te2) . (32)
The last step here is rewriting these equations to their simplest form using the Stokes number τs from equation 24, which results in
u = − vhwτs (1 + τs2) w = − vhwτs 2 (1 + τs2) . (33)
To recap, these are the radial and transverse speeds respectively relative to the gas of a loose particle in the protoplanetary disk. To illustrate, they are shown in figure 2 as a function of the particle radius and in figure 3 as a function of the Stokes number.
Figure 2: In this plot the particle speed in the radial and transverse directions are shown as a function of the particle radius, thus explicitly for both the Epstein and Stokes regimes.
2.6 Particle formation
The next problem is figuring out how many pebbles will be flying through the disk. For this I look to Lambrechts & Johansen 2014 [12], where they go into detail about how to calculate this exact property.
Lambrechts and Johansen find that the mass flux through the disk surface is
˙ Mflux= 9.5 · 10−5 Σ 5 · 103 kg m−2 M? M 1/3 Z0 0.01 5/3 t 106 yr −1/3 ME yr−1, (34) where Z0 is the dust-to-gas ratio, and t is the time since the disk initially formed
Figure 3: In this plot the particle speed in the radial and transverse directions are shown as a function of the Stokes number, thus implicitly for both the Epstein and Stokes regimes.
dust to form pebbles that start drifting radially, and find that this depends on the distance to the star. Thus, pebbles will form in a ring around the star at a radius rp, which becomes larger as time goes on, and is given by
rp =
3 16
1/3
(GM?)1/3(dZ0)2/3t2/3, (35)
where d encapsulates two growth parameters of the pebbles, namely the dust
stick-ing efficiency, and the natural logarithm of the ratio between the radius at which the pebbles start growing and the radius at which they start drifting.
The total mass formed into pebbles will then be
Mflux= πrp2· Σd, (36)
where Σd is the surface density of the dust, which is taken to be Z0· Σ, simply the
gas surface density multiplied by the dust-to-gas ratio. To get the mass flux from this equation, the derivative with respect to time needs to be taken:
˙ Mflux= π d(rp2) dt · Σd = πd(rp 2) drp drp dt · Σd = 2πrp drp dt · Σd. (37)
Then by filling in this equation until only base constants and the time dependency remain and properly normalizing it, equation 34 is obtained.
However, this equation, as mentioned, describes only the mass flux through the surface of the disk, while in this thesis we are looking at a 3D system. The conversion to 3D is not very difficult however, as by taking the gas density ρg (see equation
11) for z = 0 instead of the surface density, and by multiplying the equation by the pebble scale height Hp, one obtains the mass flux through the 3D disk, as opposed
to the 2D disk as
˙
Mflux= 2πrpHp
drp
dt · ρd, (38)
where Hp is the scale height of the pebble layer, defined as
Hp= H
r αt
αt+ τs
, (39)
with αt the turbulence factor as described in C. Ormel’s contribution to the book
’Formation, Evolution, and Dynamics of Young Solar Systems’ [13]. Note that we have here implicitly assumed that ρd is constant up to the pebble scale height. Also
note that this equation denotes the mass flux through a surface ring at the distance from the star at which pebbles form. To know the pebble flux through a surface ring at an arbitrary distance from the star, an understanding of the elements in equation 40 is necessary. In that equation, there are the factors 2πrpHp, drpdt , and ρd, which
are the surface at which pebbles form in the disk, the speed with which that surface moves (thus, inversely, also the speed at which pebbles move through that surface), and the dust density, respectively. By scaling this to the distance from the star r which we want to look at, we arrive at the mass flux through the part of the disk we want to be looking at:
˙ Mflux= 2πrpHp drp dt · ρd rdrdt rpdrpdt = 2πrHp dr dt · ρd. (40)
3
Three body system
With the individual parts of the protoplanetary disk now laid out, it is time to com-bine some of those parts into the system we will actually be looking at in this thesis. This system is a three body system consisting of the central star, a planetesimal and a pebble, as illustrated in figure 4.
If we look at this system, which will be described using cylindrical coordinates, from the stellar frame of reference, there are three forces that act on the pebble, with the first being the stellar gravity
FG,star= − GM?mp r02 ~ r0 r0, (41)
where r0 is the distance between the pebble and the star, as shown in figure 4, whose components in the φ and z directions are small enough to be neglected, giving
The second is the planetary gravity FG,planet = GM mp x2+ y2+ z2 ~r + ~φ + ~z p x2+ y2+ z2 = GM mp (x2+ y2+ z2)3/2(x ˆr + y ˆφ + z ˆz). (43)
The third is the gas drag as derived in section 2.4, which in this case can be worked out as FD= − m(~v − ~vg) ts = −vrr + (vˆ φ− vg) ˆφ + vzzˆ ts . (44)
The total equation of motion for the system is then given by
~a(ˆr, ˆφ, ˆz) = −(r + x)Ω2− GM x (x2+y2+z2)3/2 −vrts − GM y (x2+y2+z2)3/2− vφ−vg ts − GM z (x2+y2+z2)3/2 −vzts . (45) 3.1 Co-rotating frame
Now, as stated earlier, equation 45 is the equation of motion in the stellar frame of reference. However, for the simulation it will be much more convenient to look at the situation through a frame that rotates along with the planetesimal. This does however complicate the equation of motion somewhat.
First though, we switch from cylindrical coordinates to Cartesian ones, which make the simulations easier. It is relatively easy to accomplish, since the ˆr and ˆφ directions directly map onto the ˆx and ˆy directions.
Figure 4: Diagram showing the individual parts of the three body system of star, planetesimal and pebble, along with their orbital frequencies (~vr = ~Ω). Image credit to Rico Visser.
The first complication is that, in a frame that co-rotates with the planetesimal, there will be two extra inertial forces working on the pebbles. The first of these is the centrifugal force
Fcent = mpvp2 r r0 rˆr = mpΩ20(r + x)ˆr, (46)
where vp is the speed with which the frame (thus the planetesimal) is moving, and
once again the φ and z components are negligible. The second of these forces is the Coriolis force Fcor= −2~Ω0× ~v = −2 vyΩ0 −vxΩ0 0 . (47)
The second complication is that there will be a correction to the velocity of the pebbles relative to the frame. This correction will only affect the speed in the y-direction, since that is the only direction the frame is moving in. The correction then is vy, rot= vy, inert− ~Ω0× ~r = vy, inert− r0Ω0. Thus the speeds in the x and z
directions remain the same and can simply be renamed, while both the particle speed in the y direction and the gas velocity must be adjusted as vy,rot = vy,inert− r0Ω0.
This, together with the change in coordinates, means the drag force of equation 44 needs to be adjusted as FD= − vrr + (vˆ φ− vg) ˆφ + vzzˆ ts ⇒ FD,rot= − vxx + (vˆ y− (vk(1 − η) − r0Ω0))ˆy + vzzˆ ts = −vxx + (vˆ y− r 0(Ω − Ω 0) + vhw)ˆy + vzzˆ ts . (48)
With all these corrections, the equation of motion has become
~a(ˆx, ˆy, ˆz) = (r + x)(Ω02− Ω2) −(x2+yGM x2+z2)3/2 − vr ts − 2vyΩ0 − GM y (x2+y2+z2)3/2 − vy−r0(Ω 0−Ω)+vhw ts + 2vxΩ0 − GM z (x2+y2+z2)3/2 − vzts . (49)
This equation can still be reduced further, since Ω can be expressed in terms of Ω0,
thereby reducing the number of parameters, as
Ω = GM? r03/2 = r r0 3/2GM? r3/2 = r r0 3/2 Ω0. (50)
Then, using the approximation r0 ≈ r + x and using the binomial approximation, this can be written as
Ω = r 0 r −3/2 Ω0 = 1 +x r −3/2 Ω0 = 1 −3 2 x r Ω0. (51)
Figure 5: An illustration of the shear field as it impacts the velocities of entities inside our simulation box.
Using this, and a similar derivation for Ω2, the equation of motion becomes
~a(ˆx, ˆy, ˆz) = 3xΩ02−(x2+yGM x2+z2)3/2 −vrts − 2vyΩ0 − GM y (x2+y2+z2)3/2 − vy+32xΩ0+vhw ts + 2vxΩ0 − GM z (x2+y2+z2)3/2 − vz ts , (52)
if we neglect terms of order x2and higher. This then is the final form of the equation of motion, which will be used in the simulations to model the behaviour of pebbles.
3.1.1 Shear
A particular part of the equation of motion worth highlighting is what is called the shear. It is the factor 32xΩ0 in the wind resistance term in the y-direction. This
term is a result of the frame of reference moving with the Keplerian frequency Ω0.
However, the Keplerian frequency is not Ω0 for any value of x except for 0, leading
to the correction term of the shear as derived in the last section. An illustration of the shear term is provided in figure 5.
This means both the gas and the pebble velocities need to be corrected with this term. The gas velocity was already corrected in equation 52, and the initial pebble velocity can be corrected by adding the shear term to the pebble speed in the y direction, which is w from equation 33, and correcting it for the movement of the frame, making it vy, ∞= w τs2 + 3 2xΩ0 (53)
4
Numerical methods
With all the physics of the system now worked out, we can start applying it numer-ically. Quite a number of approximations, assumptions and practical considerations
Definition Value Units
r Orbital distance of planetesimal 3.74 · 1011 m
τs Stokes number 0.01 unitless
¯
m Mean molecular mass 3.99 · 10−27 kg
d Kinetic diameter 2.89 · 10−10 m
αt Turbulent strength 10−4 unitless
ρs Particle density 103 kg m−3
cre Coefficient of restitution 0.01 unitless
Table 2: This table contains the assumptions used for the free parameters r and τs,
and other assumptions based on Visser et al. 2020 [9], Weidenschilling 1977A [10], and Ormel’s contribution to [13]
have to be made to accurately yet efficiently simulate pebble accretion, and these will all be explained in this section.
First, the assumptions need to be outlined. For this thesis, certain values for cer-tain parameters have been chosen, and these have been outlined in table 2. We chose to look at situations with a fixed distance from the star, a fixed Stokes number, a fixed mean molecular gas mass, and a fixed kinetic diameter of the gas. Futhermore, a specific planetary radius was chosen for each measurement.
4.1 Rebound
To start with, the backbone of the simulation is an open source piece of software called Rebound (Rein & Liu, 2012 [14]), along with its extension package Reboundx (Tamayo & Rein, 2020 [15]). Rebound is an N-body integrator, meaning it is software designed to numerically integrate a large number of particles under the effect of some forces causing these particles to move. It is written in C, and has a Python module that can call any of its functions, which I used.
Rebound has a number of different integrators, useful for different purposes. I had Rebound use the IAS15 integrator, described in detail in Rein & Spiegel 2015 [16]. Its name is short for Integrator with Adaptive Step-size control, of the 15th order, which means it is an integrator that itself chooses a step-size that is needed to accurately calculate the forces it is given, be they velocity dependent or not. The error tolerance is determined not by an arbitrary numerical value, but by the global error estimate shown in equation 9 from the aforementioned paper by Rein & Spiegel. This means it overall guarantees both speed and accuracy, which is why it was picked.
The main code of Rebound specialises in gravitational interactions, yet it has options for adding any custom force, and the Reboundx package makes implementing these forces even easier. These functions can be written in both C and Python.
The equation of motion given in equation 52 for the pebbles was implemented as such a custom force, simply in one function, and in C, for efficiency. With this implementation, no inter-pebble interactions are taken into account, as they would take considerable computational capacity to model, while their effects would be negligible (Ormel 2017, [13]). This way, the complexity of this part of the simulation
stays at O(N ), as opposed to nearly O(N2) if interactions between particles were to be taken into account.
Rebound also has useful modules that allow programs, while integrating the move-ment of particles, to check for collisions between these particles, and to resolve these in a number of different ways.
In my simulation I had Rebound use a tree algorithm (of complexity
O(N log(N )) to search for collisions, as opposed to a direct search algorithm, as that has a complexity of O(N2), and thus would take significantly longer. What the algorithm searches for is overlapping particles that will in this time step move closer to one another.
If it finds such a particle pair, the collision is resolved by treating both particles as hard spheres, meaning they will bounce off of one another. They do this with a coefficient of restitution of 0.01, meaning the particles lose a lot of their energy in the collisions.
4.2 Pebble injection
Of course, Rebound is limited in scope and relies on user input to accurately model non-standard cases. The most important of these inputs for our model is how and when pebbles are injected into the simulation.
4.2.1 Pebble flux
In section 2.6 we derived the mass flux in the form of equations 34, for the mass flux through the disk surface, and 40 for the mass flux through any arbitrary surface ring the 3D disk. However, for the simulation it would have been highly impractical to simulate the entire disk, or even a ring in it. My simulation was confined to only a small box that rotates along with the planetesimal in its center, and therefore only the mass flux through this box was needed. Thus an equation analogous to equation 40 was necessary for the box. For this, the same procedure used to arrive at equation 40 was used. By scaling the surface represented by 2πrH to Abox = LxLz for the
surface at which pebbles enter the box, and the speed at which pebbles pass through this surface represented by drdt to vinj (the speed of injection into the box), I obtained
the mass flux through the box. Thus ˙
Mflux,box = LxLzvinjρd. (54)
Now, since ρd is a constant in both equations, ˙Mflux,box could be written in terms
of ˙Mflux, whose numerical value could be taken by filling in equation 34, where we
took the age of the disk to be 106 years. The whole equation then became ˙
Mflux,box =
LxLzvinjM˙flux
2πrHpdrdt
, (55)
where r was the distance between our planetesimal and the star around which it orbits. Then only the speeds vinj and drdt had to be defined. Luckily, that work had
already been done in sections 2.5 and 3.1.1, as vinj is simply vy,∞ from equation 53,
and drdt was u from equation 33, making the final mass flux equation ˙
Mflux,box =
LxLz(τ sw2 +32xΩ0) ˙Mflux
2πrHpu
. (56)
With equation 56 in mind, I needed to define its parameters. The size of the box I based on the concept of an impact parameter (Visser & Ormel, 2016 [17]). The impact parameter denotes the radius of a certain region where, if a pebble is released inside there, it will accrete onto the planetesimal. Lxand Lzwere chosen to
be 5 times the impact parameter in their respective directions, as pebbles released outside that region do not, to a significant degree, exhibit the behaviour that was interesting to this thesis, and can therefore be ignored.
4.2.2 Injection
Of course, the mass flux per unit time was necessary, but in my simulation it was not very efficient to check for something as often as once per second, especially since the simulation usually stretched over multiple millions of them. In the simulation I used a particular time called the crossing time, defined as the typical time it would have taken a pebble to cross our simulation box if no planetesimal perturbed its orbit. By scaling injection to this time, there is an effective maximum number of particles in the simulation at once, as after one crossing time, the amount of particles leaving the box should be equal to the number of particles entering the box. The planetesimal in the box will slightly perturb this relation, but it should still hold roughly. The amount of mass that had to be added per time crossing then was
Mtcross, box= ˙Mflux,boxtcross. (57)
The moments to inject pebbles were then chosen to be exact fractions of tcross,
namely 5001 of it. The fraction 500 was chosen as a compromise between the resolu-tion and speed of the simularesolu-tion, allowing simularesolu-tions of up to 10000 particles over short enough time scales. Thus every time tcross500 passed, M tcross,box
500 kg of pebbles
would need to be added.
Naturally, mass could not be injected directly, only in the form of pebbles, and thus we need to know the number of pebbles M tcross,box
500 kg represents. Before every
simulation I chose a particular Stokes number for the pebbles to investigate, and with that the radius of these pebbles was also predetermined. The mass per pebble was then simply the density times the volume of a pebble (assuming them to be perfect spheres). The amount of pebbles that would then need to be injected into the box per unit time thus was
˙ M flux,box
mpebble . However, this constituted too many particles to
simulate in a timely fashion.
An approximation was therefore required. Luckily, as described by Rein et al. 2009 [18], the situation could still be accurately modelled, provided two criteria are met.
First, the mean free path needed to be constant between the physical and nu-merical systems, meaning that Nnumσnum = N σ, where N is the number of particles
in a system, and σ is the collisional cross section. Second, gravitational time scales needed to remain large in comparison to collisional ones.
The first criterion could easily be satisfied by increasing the numerical radii of pebbles according to snum= r N Nnum s. (58)
An interesting result of Rein was that if Nnum were to be (much) more than 10, the
second criterion was easily satisfied if the first one was.
Thus, I was free to choose a number of particles to simulate, provided that it was larger than 10. The Nnum settled on was 2500, which was as many as could be
simulated in a timely fashion, and still maintain reasonable accuracy. A discussion of preliminary measurements relating to this is given in appendix B.
Using the injection method already discussed, but simply replacing the physical amount to be added per timestep with Nnum, I could ensure there would never be
(many) more than Nnum particles in the simulation. And since we already had the
amount of pebbles N that would have to physically be added as M tcross,box
mp , the ratio N
N num could easily be calculated, and with it snum.
Another method that could be used to calculate the numerical radius of the pebbles is snum2=6.42 · 109 m2 LxLyLz RH3 R 100km 3 · 104 Nnum s 0.01m ˙ M 500 MEMyr−1 h r au i1/4 vhw 32m s−1 αt 10−4, (59)
where R is the radius of the planetesimal, RH is the hill radius of said planetesimal,
defined as r3MM p
?
1/3
where r is the distance at which the planetesimal orbits around the central star, as it is in equation 59. This equation was derived by my daily supervisor Rico Visser for the purposes of this thesis, and will serve as a check for the accuracy of the simulations.
4.3 Numerical values
Naturally, to arrive at quantitative conclusions, constants have to be filled in, and thereby quantitative definitions for parameters have to be derived. Using the con-stants outlined in table 1, the values for the parameters in table 3 have been calcu-lated.
5
Results
5.1 Model testing
With both the model used for the simulation and the theory behind it explained, I can now move on to the results that were obtained.
First, it was important to test whether the physics in our model worked as they should. For this, I turned to a paper by my daily supervisor, Visser et al. 2020 [9].
Definition Value Units equation
T Temperature 107.5 K eq. 1
Σ Gas surface density 3289 kg m−2 eq. 2
H Scale height of gas in disk 1.21 · 1010 m eq. 3
Ω0 Keplerian frequency planetesimal 5.036 · 10−8 s−1 vkr
vhw Headwind velocity 32 m s−1 ηvk
ρg Gas density 1.08 · 10−7 kg m−3 eq: 11
˙
Mflux Mass flux 247 ME Myr−1 eq: 34
snum Numerical pebble radius 2.47 · 105 m eq: 59
u Radial pebble speed 0.64 m s−1 eq: 33
w Transverse pebble speed (shearless) 32 m s−1 eq: 33
Rpa Radius onset pebble accretion 1.83 · 105 m eq: 60
Table 3: This table contains the parameters calculated from the constants and assumptions from table 1 and table 2.
That contained the plots that can be seen on the left side of figure 6, which I set out to replicate, which, as is visible on the right side, succeeded. These experiments show that my model works at the edges of a number of different regimes by testing for a number of different planetary radii. First, with the 200 km test, the lower limit of pebble accretion is tested. Pebble accretion starts at the radius of
Rpa= 160 km · vhw 30 m s−1 ρs 2500 kg m−3 −0.36 r au 0.42τs 0.1 0.28 , (60)
in the limit where τs< 1 (Visser et al. 2016, [17]). Below this radius is what I will
be calling the ballistic regime, where the planetary gravity is not large enough to affect the pebble trajectories much past their ballistic trajectories.
Around this radius pebbles are mostly focused by the planet, not so much ac-creted, whereas for increasing planetary radii more and more pebbles are accreted instead of focused by the planets. This is the regime this thesis will be looking at, as the focused particles have the largest probability of colliding with one another and accrete, where they otherwise would not.
Second, with the 400 km test the start of the main pebble accretion regime, also called the wind regime (Ormel 2017, [13]) is tested. Here pebbles are mainly accreted, not focused, and no other factors except gravity play a dominant role. The reason pebbles are mainly accreted rather than focused with a larger planet radius is that the focusing length over which pebbles would be focused scales differently with the planet radius than the impact parameters (thus the area from which particles would accrete onto the planetesimal) do. Exactly how this scales differently however, falls outside of the scope of this thesis.
Third and last, with the 1000 km test the start of the shear dominated regime is tested. In this regime the shear term discussed in section 3.1.1 starts to dominate the other terms in the relevant velocities, and one even sees the co-rotation line, the line where vhw and the shear are equal and opposite, past which pebbles are
’pushed’ in the opposite direction (from the point of view of the planet, remember that this is a co-rotating frame and the shear is an artifact of that).
It is important to check each of these regimes, as only if all regimes are described correctly, can the model be said to work correctly. Even though only the first regime will be analyzed in this thesis, the model used is valid in all regimes, as can be seen by the results in figure 6.
With this it was clear that the physics worked as they should, and the actual research could begin.
5.2 Initial conditions
For this thesis, 8 simulations were run. These were for four planetary radii, which were 80, 100, 200 km and the radius for the onset of pebble accretion, as calculated from equation 60 and given in table 3. Each planetary radius was run once with collisions, and once without collisions, to assess possible differences in mass accretion and spin accumulation.
Parameter Value all x0 33.25 z0 0 yrelease 1662.23 200 km Lx 36.45 Lz 36.77 Rpa Lx 35.70 Lz 36.20 100 km Lx 21.28 Lz 23.40 80 km Lx 20.08 Lz 24.93
Table 5: Release dimensions of pebbles in terms of planet radius.
Each simulation ran for an in-simulation time of 5tcross so that an equilibrium could be established
in each case. At 2500 evenly spaced points in time during these simulations data was gathered, and new pebbles were injected as discussed in section 4.2.2. Each simulation was configured so that 2500 pebbles would be added each tcross, which makes
12500 particles in total through the system. The pebbles were released within the ellipsoid area formed by their release parameters Lxand Lz,
which are shown in table 5, around x0and z0, which
were experimentally defined as a multiple of their impact parameters. The release parameters were chosen such that if a pebble were released at the maximum extent in the positive direction, it could still collide with a particle released at the other extreme. Beyond these parameters I found the fo-cusing to be too weak to lead to significant results. The exact starting positions of the pebbles in the x-direction were determined by a weighted for-mula provided to me by my daily supervisor Rico Visser, which can be seen in code in appendix A, while in the z-direction it was a simple uniform random fraction of the release pa-rameters. The initial starting speeds are then w + the shear contribution using the
Radius Simulation Equation 59
200 km 534 km 494 km
Rpa 462 km 425 km
100 km 115 km 187 km
80 km 82.5 km 185 km
Table 4: Numerical radii of pebbles in simulation against the values obtained from equation 59.
Figure 6: On the left side: plots of pebbles of Stokes number 0.1 from Visser et al. 2020 [9], where blue lines are pebbles that accrete, and gray ones do not.
On the right side: plots of pebbles of Stokes number 0.1 from my model, where red lines are pebbles that accrete and gray ones do not.
generated starting x-coordinate and u from table 3.
The physical radius of the pebbles was 3.94 cm, which was the size that gave a Stokes number of 0.01. This radius was in the simulation blown up following equation 58. The numerical radii in the simulations and the ones determined from equation 59 are both shown in table 4.
5.3 Ballistic regime
Four of the simulations covered planets just below the size of the onset of pebble accretion. These are the tests at 80 and 100 kilometers. The parameters that are of importance to the conclusion of this thesis are the mass accreted in kg and the spin accumulation without units. The mass accreted is simply the number of numerical particles that accrete, times the mass per particle, times the ratio between physical and numerical particles. The spin accumulation is defined as the spin conferred on the planetesimal by impacting particles, and is given by
S = (vyx − vxy)
normalisation, (61)
where for the normalisation for the spin accumulation we chose the escape angular momentum of the planetesimal
√
2GM Rp
3.5 . A positive value for this means a
pro-grade spin contribution, a negtive value means a retropro-grade spin contribution. This formula was provided by my daily supervisor, Rico Visser.
The results of these simulations are shown qualitatively in figure 7, and quanti-tatively in table 6.
The data for 80 km does suggest a significant difference in spin accumulation, and somewhat less so in mass accretion. The data however does not conform to prior research, as Visser et al. 2020 [9] found that, at least without collisions, spin accumulation should average to zero. This is something the data for 100 km planetary radius does seem to conform to.
5.4 Onset pebble accretion regime
The other four simulations cover the onset of the pebble accretion regime, at the exact radius that happens following equation 60, and 200 km. The same parameters as in the last regime are of importance here.
The results are shown quantitatively in table 6 again, and qualitatively in figure 8. Clearly visible is that there is a larger absolute, but smaller relative difference in mass accreted by the planetesimals with and without collisions in this regime than in the ballistic regime. The spin accumulation is relatively close to 0, which could imply a relatively low contribution to spin evolution by pebbles of this size in this regime.
Figure 7: Parameter plots for the ballistic regime. ’With’ means that the data is from a simulation with collisions, ’without’ means the opposite.
Parameter Collision Collisionless 200 km Macc 1.82 · 1016 1.92 · 1016 S 0.040 -0.019 Rpa Macc 1.19 · 1016 1.15 · 1016 S -0.047 0.083 100 km Macc 1.73 · 1014 1.67 · 1014 S 0.450 0.048 80 km Macc 6.78 · 1013 6.43 · 1013 S 0.638 -0.402
Table 6: Data obtained from simulations. Units for Macc and S are kg and none,
Figure 8: Parameter plots for the accretion regime. ’With’ means that the data is from a simulation with collisions, ’without’ means the opposite.
6
Discussion
The methods outlined in this thesis have the potential to lead to interesting and sig-nificant results, and the results attained for this thesis have been sigsig-nificant. They suggest an interestingly large and significant contribution to spin accumulation at a radius of 100 km, and differ significantly from previous research for 80 km. At the same time, the data suggests little to no contribution from collisions at the onset of the pebble accretion regime, which is certainly a significant if not terribly exciting result.
Yet the data sets do have a number of problems. For example, the ballistic data sets have a problem of sample size, as for 80 km only 79 and 75 pebbles accreted with and without collisions respectively, while for 100 km only 130 and 125 did for the same cases. This does cast some doubt on the significance of these results, as it is doubtful whether with this small number of particles the mean has converged to the ’real’ average to a significant enough degree.
The data sets from the onset of pebble accretion do have a larger sample size, as for Rpa 1024 pebbles accreted with and 986 without collisions, while for 200 km 1269
pebbles accreted with and 1339 without collisions. Thus it is much more likely these will have converged, and thus the results in this regime are much more statistically significant.
Particularly for the ballistic regime then it remains unclear how large random fluctuations might be because of this small sample size. The small sample size might explain the divergence from earlier research found in these simulations, yet more or longer tests would have to be run to be sure. Since really only one data set for each radius has been obtained, averaging over multiple data sets with initial conditions identical except for the exact release locations of the pebbles will certainly help eliminate statistical errors to a significant degree.
A problem not related to the size or number of the data sets gathered is the dis-crepancies clearly visible in table 4, where the numerical radii of the pebbles in the simulation do not particularly match in the ballistic regime. The values in the accretion regime fall within 10 % of one another, which is a relatively small differ-ence that falls within expected expected numerical error margins with the numerical approximations that had to be made.
The discrepancy in the ballistic regime seems to imply that either the theoretical or the numerical model might start breaking down in the ballistic regime. This is something that bears further investigation.
A similar discrepancy between numerical and theoretical values is that in equa-tion 56, ˙Mflux is taken to be the result from equation 34. Yet this might well not
be entirely warranted, as that equation is derived for a case where all material is assumed to be in a 2D plane, as opposed to the 3D cloud that is simulated. It is assumed this can be done as at all times Lz remains much smaller than H and
Hp, meaning the effects of the material thinning out in the z-direction should be
negligible for the purposes of these simulations. Finding an expression analogous to equation 34 for the 3D disk would eliminate the need for this assumption though, and would further increase accuracy.
The last numerical value I want to discuss is the number of particles simulated. The optimal number for this in our simulations is as many as is practically feasible to simulate. As discussed in appendix B, a higher number of numerical particles pushes the simulation closer to the physical situation. Yet it has diminishing re-turns, as increasing Nnum asymptotically converges to the physical system. Here
Nnum= 2500 was chosen as it was relatively close to the highest tested value as seen
in figure 9, and these tests could be completed and analyzed within the remaining time, making it the largest Nnum that was practically feasible.
The last point for discussion is the rather limited scope of this thesis. A considerable amount of time and effort went into the accuracy and efficiency of the simulation, and a number of issues including limited computer memory and long run times alongside a rather too ambitious first few attempts at simulations, limited the amount of time available for actual test simulations.
Further opportunities for research include simulations at extra planetary radii, additional Stokes numbers or various orbital radii. Preliminary research has already suggested significant spin accumulation at a Stokes number of 0.1 with a plane-tary radius of 250 km, yet no actual test runs could be completed because of time constraints.
Finally, running these simulations on something more powerful than a mid-range three-year old laptop such as a research cluster or even a supercomputer if one could get access, might greatly accelerate data collection.
In short, this thesis could have benefited from additional time and resources, while the accuracy of the model in one of the investigated limits could stand to bear further scrutiny. That said, the methodology has been meticulously thought out and iterated upon, and is ready to be used again at a moment’s notice for further data acquisition.
7
Conclusion
The data for the ballistic regime suggests significant effects of collisions, primarily on spin accumulation, and bears further investigation to ensure these are not simply statistical anomalies.
At the same time, the data for the onset of the pebble accretion regime suggests little to no significant effects of collisions on either of the investigated parameters, which is of course a significant result in and of itself.
The data overall seems to suggest that this is an interesting avenue to explore for further research, potentially yielding significant contributions to spin accumulation, not only for the surveyed parameter values, but also for different Stokes numbers and planetary radii than surveyed here.
8
Bibliography
References
[1] Jr. Haisch, Karl E., Elizabeth A. Lada, and Charles J. Lada. Disk Frequencies and Lifetimes in Young Clusters. apjl, 553(2):L153–L156, June 2001.
[2] C. W. Ormel and H. H. Klahr. The effect of gas drag on the growth of pro-toplanets. Analytical expressions for the accretion of small bodies in laminar disks. aap, 520:A43, September 2010.
[3] M. Lambrechts and A. Johansen. Rapid growth of gas-giant cores by pebble accretion. aap, 544:A32, August 2012.
[4] D. Prialnik. An Introduction to the Theory of Stellar Structure and Evolution. Cambridge University Press, 2009.
[5] A. Johansen, E. Jacquet, J. N. Cuzzi, A. Morbidelli, and M. Gounelle. New Paradigms for Asteroid Formation, pages 471–492. University of Arizona Press, 2015.
[6] S. J. Weidenschilling. The Distribution of Mass in the Planetary System and Solar Nebula. apss, 51(1):153–158, September 1977.
[7] C. Hayashi, K. Nakazawa, and Y. Nakagawa. Formation of the solar system. In D. C. Black and M. S. Matthews, editors, Protostars and Planets II, pages 1100–1153, January 1985.
[8] Y. Nakagawa, M. Sekiya, and C. Hayashi. Settling and growth of dust particles in a laminar phase of a low-mass solar nebula. Icarus, 67(3):375–390, September 1986.
[9] R. G. Visser, C. W. Ormel, C. Dominik, and S. Ida. Spinning up planetary bodies by pebble accretion. Icarus, 335, January 2020.
[10] S. J. Weidenschilling. Aerodynamics of solid bodies in the solar nebula. mnras, 180:57–70, July 1977.
[11] F. L. Whipple. On certain aerodynamic processes for asteroids and comets. In Aina Elvius, editor, From Plasma to Planet, page 211, January 1972.
[12] M. Lambrechts and A. Johansen. Forming the cores of giant planets from the radial pebble flux in protoplanetary discs. aap, 572:A107, December 2014. [13] Chris W. Ormel. The Emerging Paradigm of Pebble Accretion, pages 197–228.
Springer International Publishing, Cham, 2017.
[14] H. Rein and S. F. Liu. REBOUND: an open-source multi-purpose N-body code for collisional dynamics. aap, 537:A128, January 2012.
[15] Daniel Tamayo, Hanno Rein, Pengshuai Shi, and David M. Hernand ez. RE-BOUNDx: a library for adding conservative and dissipative forces to otherwise symplectic N-body integrations. mnras, 491(2):2885–2901, January 2020.
[16] Hanno Rein and David S. Spiegel. IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits. mnras, 446(2):1424–1437, January 2015.
[17] Rico G. Visser and Chris W. Ormel. On the growth of pebble-accreting plan-etesimals. aap, 586:A66, February 2016.
[18] H. Rein, G. Lesur, and Z. M. Leinhardt. The validity of the super-particle approximation during planetesimal formation. aap, 511:A69, February 2010.
A
Weighted randomizer
Python code provided by my daily supervisor, Rico Visser.
Release distance in the x-direction as a function of headwind, the left x-bound, right x-bound, Ω0, a random number between 0 and 1, and τs.
def cumdist(vhw, x1, x2, om, x, tau): term1 = vhw / (1 + tau**2)
term2 = (x2 - x1) term3 = 0.75 * om term4 = (x2**2 - x1**2)
N = (term1 * term2) + (term3 * term4) part1 = 2 * N * (3 * om)**-1
part2 = ((3 * om * x * N**-1) + (term1**2 * N**-2)
+ (3 * om * N**-2 * (term1 * x1 + 0.75 * om * x1**2))) part3 = -term1 * N**-1
B
Convergence test for average collisions per numerical
particle
A simulation at planetary radius 200 km with Stokes number 0.01 was run for Nnum
1000, 2500, 5000 and 10000, to check how quickly the average number of collisions per particle over time converged to the physical value. The results can be seen in figure 9, and imply that the difference in this value for a difference in Nnum of a
factor 10 is roughly a factor 2. They also show that an Nnum of 1000 is too low, as
increase in collisions per particle visible in figure 9 shows that at their generation, the pebbles already overlap and thus collide, which is not physical. The initial idea was to use Nnum of 10000 for the simulations, but due to time constraints this had
to be reduced.
Figure 9: Average number of collisions per particle over time to show the convergence of this value to the physical value.
