Dust Evolution in Protoplanetary Disks on the Path to Planet Formation Master’s thesis in Physics MARKUS HJÄLT DEPARTMENT OF SPACE, EARTH AND ENVIRONMENT CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2024 www.chalmers.se www.chalmers.se Master’s thesis 2024 Dust Evolution in Protoplanetary Disks on the Path to Planet Formation MARKUS HJÄLT Department of Space, Earth and Environment Division of Astronomy and Plasma Physics Chalmers University of Technology Gothenburg, Sweden 2024 Dust Evolution in Protoplanetary Disks on the Path to Planet Formation MARKUS HJÄLT © MARKUS HJÄLT, 2024. Supervisor: Dr. Timmy Delage, Department of Space, Earth and Environment Examiner: Prof. Jonathan Tan, Department of Space, Earth and Environment Master’s Thesis 2024 Department of Space, Earth and Environment Division of Astronomy and Plasma Physics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Dust density shown as function of particle sizes and distance from the star. The figure is taken from a simulation with fragmentation velocity 1 ms´1 at time 1 Myr. Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Printed by Chalmers Reproservice Gothenburg, Sweden 2024 iv Dust Evolution in Protoplanetary Disks on the Path to Planet Formation MARKUS HJÄLT Department of Space, Earth and Environment Chalmers University of Technology Abstract Since the discovery of the first exoplanet 51 Pegasi b in 1995 thousands of exoplanets have been found. The Kepler mission, which goal was to look for the first Earth-like exoplanets, found that multiplanetary systems with super-Earth and sub-Neptune size planets seem to be very common. However, the early planet formation theories of formation in the outer protoplanetary disk followed by an inward migration could not explain these systems in which the planets were found in tight non-resonant orbits. Hence the Inside-Out Planet Formation theory was proposed in which the planets are formed in situ at a pressure maximum. This pressure trap, caused by magnetorotational instabilities, halt the radial drift of pebbles and the pebbles can coagulate into protoplanets. Since the grain growth process and the radial drift is of importance for this theory these processes have been studied in detail in this thesis. Simulations have been done of the dust evolution around a solar mass T- Tauri star using the Python-based program DustPy. Three disk models have been studied with three different fragmentation velocities, 1 ms´1, 5 ms´1 and 10 ms´1. The models adopted a gas accretion rate of 9Mgas “ 10´9 Mdyr´1 and a viscosity parameter of α “ 10´4 . It was found that higher fragmentation velocities allowed for larger grain sizes in young disks, but that the representative grain sizes seemed to converge at late times, independent of the fragmentation velocity. At 1 au the dominant grain sizes found were 0.59 cm for 1 ms´1 and 11.06 cm for 10 ms´1 at time 0.1 Myr. At time 5 Myr the grain sizes had almost converged with the dominant grain sizes being 0.59 cm for 1 ms´1 and 0.53 cm for 10 ms´1. A similar tendency could be seen in the outer disk at 10 au where the dominant grain sizes at 0.1 Myr were 0.27 cm and 2.69 cm for the 1 ms´1 and 10 ms´1 fragmentation cases respectively. At time 5 Myr the representative sizes had decreased and converged to 0.012 cm and 0.011 cm for the 1 ms´1 and 10 ms´1 fragmentation cases respectively. When studying the dust flux and mass delivery to a region within 1 au it was found that the dust mass delivered after 1 Myr was 77 M‘ for 1 ms´1 and 99 M‘ for both the 5 ms´1 and 10 ms´1 fragmentation case. The mass delivered by grains larger then 0.06 cm was found to be 94.4%, 99.8% and 99.9% in the different fragmentation cases respectively. The time required to deliver about 3 M‘ to the region within 1 au was also found to be 9 ˆ 104 yrs, 8 ˆ 103 yrs and 5 ˆ 103 yrs respectively for the different fragmentation cases. Keywords: protoplanetary disk, planet formation, Inside-Out Planet Formation, grain growth, Smoluchowski equation, DustPy. v Acknowledgements First of all I would like to thank Prof. Jonathan Tan for inviting me to his research team and giving me the opportunity to carry out this master thesis project. I would also like to thank my supervisor Dr. Timmy Delage whose knowledge, support and guidance has helped me navigate through this project. Furthermore, I would also like to direct a big thanks to my family who has always been there with their endless support and encouragement throughout my studies. Markus Hjält, Gothenburg, June 2024 vii List of Acronyms Below is the list of acronyms that have been used throughout this thesis listed in alphabetical order: DZIB Dead-Zone Inner Boundary GI Gravitational Instability GMC Giant Molecular Cloud IOPF Inside-Out Planet Formation ISM Interstellar Medium MHD Magnetohydrodynamics MRI Magnetorotational Instability PDF Probability Density Function (in thesis not capitalized) STIP System with Tightly packed Inner Planets ix Nomenclature Below is the nomenclature of indices, parameters, and variables that have been used throughout this thesis. Indices k,l,p,q Indices for summation in the discrete Smoluchowski equation Parameters B Magnetic field vgas Gas velocity η Magnetic diffusivity Ω Angular velocity σ Conductivity β Hall parameter pgas Thermal gas pressure T Temperature ρgas Gas mass volume density cs Sound speed Φ Gravitational potential W Shear stress tensor ΩK Keplerian angular frequency ηgas Pressure gradient parameter vK Keplerian velocity Hgas Gas vertical scale height Σgas Gas surface density Trϕ Shear stress tensor in rϕ-plane xi α Disk viscosity parameter ν Kinematic disk viscosity ρs Dust bulk density vdust Dust velocity ρdust Dust mass volume density τs Dust particle stopping time St Stokes number Hdust Dust vertical scale height δz Verticle mixing parameter δr Radial mixing parameter Ddust Dust diffusion Σdust Dust surface density Ftot Total dust flux Fadv Advection flux Fdiff Diffusion flux ndust Dust number density Variables t Time r Radial direction ϕ Azimuthal direction z Vertical direction a Grain radius xii Constants Below is the list of constants that have been used throughout this thesis: c “ 2.99792458 ˆ 108 ms´1 Speed of light e “ 1.60218 ˆ 10´19 C Elementary charge kB “ 1.380649 ˆ 10´23 JK´1 Boltzmann constant G “ 6.67430 ˆ 10´11 Nm2kg´2 Gravitational constant Md “ 1.99 ˆ 1030 kg Solar mass Rd “ 6.96 ˆ 108 m Solar radius xiii Contents List of Acronyms ix Nomenclature xi Constants xiii List of Figures xvii List of Tables xix 1 Introduction 1 1.1 The Quest of Our Origins . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Star Formation and Protoplanetary Disks . . . . . . . . . . . . . . . . 2 1.2.1 Planet Formation Theories in Protoplanetary Disks . . . . . . 3 1.3 The Search for Earth-like Exoplanets . . . . . . . . . . . . . . . . . . 4 1.4 Inside-Out Planet Formation . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Theory 11 2.1 Magnetohydrodynamics and the MRI . . . . . . . . . . . . . . . . . . 11 2.1.1 Ideal MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.2 Non-ideal MHD . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Gas in Protoplanetary Disks . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1 Thermal Gas Structure . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2 Dynamical Gas Structure . . . . . . . . . . . . . . . . . . . . 14 2.2.3 Gas Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Dust in Protoplanetary Disks . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 Dust Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.2 Grain Growth Process . . . . . . . . . . . . . . . . . . . . . . 23 3 Methods 25 3.1 DustPy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1.1 Simulation Setup and Parameters . . . . . . . . . . . . . . . . 26 3.1.1.1 Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.1.2 Star . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.1.3 Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1.1.4 Dust . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 xv Contents 3.1.2 Simulation Boundary Conditions . . . . . . . . . . . . . . . . 29 3.2 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4 Results & Discussion 31 4.1 Grain Size Distributions . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Dust Fluxes and Mass Delivery . . . . . . . . . . . . . . . . . . . . . 39 4.3 Affect of Radial Resolution . . . . . . . . . . . . . . . . . . . . . . . . 41 5 Conclusions 43 5.1 Dominant Grain Sizes . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.2 Dust Mass Flux and Mass Delivery . . . . . . . . . . . . . . . . . . . 43 5.3 Radial Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Bibliography 45 A Appendix 1 I A.1 Proof of Equal Expressions in Gas Evolution . . . . . . . . . . . . . I A.2 Proof of Expression for Radial Gas Velocity . . . . . . . . . . . . . . II xvi List of Figures 1.1 A sketch over the different classes of the star forming process. To the left Class 0 shows the first collapse into the first core. Class I shows the second collapse and how a disk is formed around the protostar. Class II shows the pre-main-sequence star with the disk around it. Most of the gas and dust from the first core has been dispersed. To the right the Class III phase is shown where only the star and a debris disk is left from where rocky and icy planets could grow. . . . . . . . 3 1.2 A sketch over how the MRI arises. Two fluid elements are linked vertically by a magnetic field line. In a) the elements are slightly per- turbed radially shown by the black arrows. In b) the shear leads to an azimuthal displacement shown by the orange arrows. This leads to magnetic tension, shown by the green arrows in c), through which angular momentum is transferred. This causes the radial displace- ment of the fluid elements to increase, shown by the black arrows in d), leading to an instability. Inspiration to this sketch is taken from [1]. 6 1.3 A sketch over the IOPF theory. The arrow shows a pebble drift towards the DZIB (top row). Pebbles get trapped at the pressure maximum (second row). At the pressure maximum pebbles coagulate and grow into protoplanets (third row). This clear the region of dust, allowing for stellar irradiation further out in the disk. The DZIB retreat and causes a new pressure maximum at which pebbles get trapped, eventually forming an additional planet (bottom row). This process continues as long as there is a supply of dust and pebbles. Inspiration to this sketch is taken from [2]. . . . . . . . . . . . . . . . 8 3.1 The initial gas profile where the gas surface density is shown in the top row and the initial gas temperature profile is shown in the bottom row. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1 Probability density functions based on the simulation with the frag- mentation velocity of 1 ms´1. The pdfs are shown at two different radii and at each radii five times are shown between 0.1 Myr and 5 Myr. The peaks represent the most probable grain size at that given time and radius. . . . . . . . . . . . . . . . . . . . . . . . . . . 32 xvii List of Figures 4.2 2D plots showing the density evolution for given grain sizes and given radii at four times between 0.1 Myr and 5 Myr for the fragmentation velocity 1 ms´1. The blue line shows the fragmentation barrier. The green line shows the drift barrier. The white line shows where St=1. . 33 4.3 Probability density functions based on the simulation with the frag- mentation velocity of 10 ms´1. The pdfs are shown at two different radii and at each radii five times are shown between 0.1 Myr and 5 Myr. The peaks represent the most probable grain size at that given time and radius. . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.4 2D plots showing the density evolution for given grain sizes and given radii at four times between 0.1 Myr and 5 Myr for the fragmentation velocity 10 ms´1. The blue line shows the fragmentation barrier. The green line shows the drift barrier. The white line shows where St=1. . 36 4.5 The dust density is shown over time at two different radii, 1 au (to the left) and 10 au (to the right) for the three different fragmentation velocities, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green). Note the different scales on the y-axis between the two radii. . . . . . . . . 37 4.6 The pdfs for grain sizes shown for the three different fragmentation velocities, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green), at 1 au for two different times, 0.1 Myr (to the left) and 5 Myr (to the right). 37 4.7 The pdfs for grain sizes shown for the three different fragmentation ve- locities, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green), at 10 au for two different times, 0.1 Myr (to the left) and 5 Myr (to the right). 38 4.8 Graphs over the mass flux through the 1 au radius over time (to the left) and graphs over the total mass delivered to the region within 1 au (to the right). In both cases three different fragmentation velocities are shown, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green). . . . . 39 4.9 Mass weighted pdfs over the dominant grain sizes for the mass deliv- ered to the region within the 1 au radius at three different times for the three different fragmentation velocities. . . . . . . . . . . . . . . . 40 4.10 Total mass delivered through the 1 au radius over time. Three frag- mentation velocities are shown, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green). The 3 M‘ level is also shown (black dashed line). . . 41 4.11 Probability density functions for grain sizes shown at two different radii, 1 au (to the left) and 10 au (to the right). At each radius two times are shown 0.1 Myr (solid line) and 5 Myr (dashed line). For each radius and time two graphs are shown, one based on simulations with 64 radial grid cells (blue graphs) and one based on simulations with 128 radial grid cells (red graphs). . . . . . . . . . . . . . . . . . 42 xviii List of Tables 3.1 The simulation snapshots and spacing summarized. . . . . . . . . . . 26 3.2 The initial stellar parameters used in the simulations. . . . . . . . . . 26 3.3 The initial gas parameters used in the simulations. . . . . . . . . . . 27 3.4 The initial dust parameters used in the simulations. . . . . . . . . . . 28 4.1 The representative grain sizes for the fragmentation velocities 1 ms´1 and 10 ms´1 at 1 au and 10 au at times 0.1 Myr, 1 Myr and 5 Myr summarized. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 xix List of Tables xx 1 Introduction In this chapter the relevance of the work carried out in this thesis is motivated and put into perspective. I start discussing the quest of our origins and planet formation theories broadly and finally narrowing it down to the objectives for this thesis. 1.1 The Quest of Our Origins For millennia, humans have cast their gaze upward, mesmerized by the night sky filled with stars spread across the plane of the Milky Way galaxy. Our ancestors might for thousands of years have looked up to the sky, being amazed by the lu- minous dots, pondering mysteries far beyond their comprehension. Questions like what those distant luminaries could be, what the big, bright, luminous dot was that appeared during the days spreading heat and light, and likely also questions about life itself and our origins. Answers to these questions have over time been found thanks to the curious nature of humans. The geocentric view of the Universe was broken thanks to obser- vations of bodies in our solar system. As the telescopes got better and the technology advanced new planets orbiting the Sun was also found. The idea of that the planets orbiting our Sun could be the only planets existing and that our solar system would be unique stood as far as until 1995, when the first exoplanet, 51 Pegasi b, was discovered [3]. This newly found planet was nothing like what would have been expected based on the knowledge we had at that time, solely based on our own solar system. This was a Jupiter-sized planet with an orbital period of about 4 days [3], which was nothing like what had been observed earlier due to no analogue in our solar system. Since the discovery of the first exoplanet multiple planets of this kind, known as hot Jupiters, have been found [4]. However, the discovery of 51 Pegasi b broke down the idea that all planetary systems would look somewhat similar to our own and indicated that there might be a huge variety in the architecture of exoplanetary systems. After the discovery of the first confirmed exoplanet the search for additional exoplanets and planetary systems beyond our own increased drastically. Multiple missions and methods were developed in order to find and observe new planetary system candidates. One of these methods is the Transit method used by NASA’s Kepler mission, which objective was to find and observe the first Earth-like exoplanet around a solar-like star [5]. This mission and its outcome is discussed in further detail in Section 1.3. 1 1. Introduction The data collected by all the missions designed to look for exoplanets have made us realise that most of the stars seem to have planetary systems and that our own solar system is far from being alone. This insight lead to questions about planet formation and consequently how the Earth once upon a time was formed. One popular theory for planet formation is the Nebular hypothesis which proposes that the planets form from gaseous disks filled with dust around newly formed stars [6][7]. The Nebular hypothesis suggests that the dust grains in these disks over time will stick to each other and coagulate, eventually forming planetary cores which further could accrete more mass finally forming a planet [7]. Hence these gaseous and dusty disks seem to be of great importance for the formation of planetary systems since they could make up the natal environment for planets. 1.2 Star Formation and Protoplanetary Disks The interstellar medium (ISM) is filled with large clouds made up of interstellar gas and dust. These giant molecular clouds (GMCs) consist dominantly of molecular hydrogen and can have masses up to thousands of solar masses and a mean molecular density of about „103 cm´3. However, the GMCs are not homogeneous, but rather contains multiple regions in which the density is much higher and in which the cloud eventually becomes gravitationally unstable and will start an inward collapse under its own gravity. Such a collapse is the first step in the star forming process and hence these regions, with a density up to „ 105 cm´3, are known as star forming clumps [8][9]. The star forming clumps are at first supported against their weight by thermal pressure, turbulent gas motion and magnetic fields, but at a certain point, when the cores become dense enough to reach the Jeans instability criterion, they become unstable and the collapse will commence [10]. Initially the collapsing cloud is optically thin to the thermal radiation emitted by the dust grains. This causes the cloud to start its collapse isothermally since the cooling rate by thermal emission is much larger than the compressional heating rate. Eventually the compressional heating rate will become greater than the radiative cooling rate halting the isothermal collapse and in the center form an adiabatic core, the so-called first core [11]. During the formation of the first core the cloud will reduce its size considerably. This leads to a decrease in the moment of inertia of the cloud, and since the angular momentum is conserved, the collapsing cloud will increase its angular velocity leading to the gas and dust flattening out into the form of a disk around the newly formed first core. This is known as the Class 0 phase of the star formation and is demonstrated to the left in Figure 1.1. When the central density and the central temperature of the first core reaches ρc „10´7 g cm´3 and Tc „ 2000 K the molecular hydrogen start to dissociate. This is an endothermic reaction which causes the temperature of the gas to decrease again. This leads to a second collapse which is halted when the gas becomes adiabatic and the thermal pressure becomes great enough to counteract the collapse. At the center of this second collapse a second core, in hydrostatic equilibrium, is formed. This second core is called a protostar, onto which material from the first core is falling [12][13]. However, due to the high angular velocity of the gas and dust, part of the gas and dust will not be accreted onto the protostar, but rather end up in a thin 2 1. Introduction Class 0 Class I Class II Class III Figure 1.1: A sketch over the different classes of the star forming process. To the left Class 0 shows the first collapse into the first core. Class I shows the second collapse and how a disk is formed around the protostar. Class II shows the pre- main-sequence star with the disk around it. Most of the gas and dust from the first core has been dispersed. To the right the Class III phase is shown where only the star and a debris disk is left from where rocky and icy planets could grow. magnetized disk around the protostar, a so-called protoplanetary disk. This stage is known as the Class I phase of the star formation, shown as Class I in Figure 1.1. The gas and dust from the first core will now be dispersed by accretion onto the protostar and ejected by bipolar outflows. What is left after the dispersion of the first core is the protostar with an opaque magnetized disk of gas and dust orbiting around it. Throughout this disk the accretion of dust and gas onto the protostar is proceeding continuously [14][15]. The star has now become what is known as a T-Tauri star and it has reached the point at which it can be called a pre-main- sequence star [8] [16]. This phase is shown as the Class II phase of the star formation in Figure 1.1. Eventually all the gas in the protoplanetary disk will have dispersed, leaving a disk containing mostly dust. At this point the pre-main-sequence star does no longer accrete a significant amount of matter, leaving a dusty disk around it, which now is called a debris disk. This debris disk make up the environment in which rocky and icy planets can continue to grow during the „ 100 Myr lifespan of the debris disk [17]. This final stage of the star formation is known as the Class III phase, shown to the right in Figure 1.1. 1.2.1 Planet Formation Theories in Protoplanetary Disks There exist multiple theories for how planets could form in protoplanetary disks. Two of the most popular theories are the Gravitational Instability (GI) and the Core Accretion. The Gravitational Instability is a consequence which arises from the self-gravity of a protoplanetary disk. The self-gravity refers to the diffuse gas gravitational influence on itself, which should not be confused with the external gravitational force exerted on the gas by external objects as, for instance, the protostar. To which degree a disk is self-gravitating can be quantified by Toomre’s Q-parameter where self-gravity becomes increasingly important as Q approaches 1 from above 3 1. Introduction [18]. When a protoplanetary disk becomes gravitationally unstable it can fragment into multiple denser parts, known as protoplanetary clumps. These clumps can then accrete gas resulting in a direct formation of gas giants over a short timescale of several hundreds of years [19]. In this way the GI theory can form planets in the fragmented protoplanetary disk. The Q-parameter has a radial dependence and decreases with increasing radius resulting in the disk being more likely to become gravitational unstable in the outer regions [1]. The Core Accretion theory explains planet formation as a bottom-up process in which initially a seed object forms in the disk. This seed object will over time grow into a solid core eventually creating the embryos for terrestrial planets. At a certain critical mass, some of these cores will have a gravitational field strong enough to initiate a runaway accretion of gas onto the core enabling the formation of gaseous planets [1]. To understand the planet formation in protoplanetary disks through these two theories one need to understand the evolution of these disks which are closely related to the angular momentum and angular momentum transport within the disks. The two most important transport mechanisms of angular momentum are disk winds and turbulence. In the theory of disk winds angular momentum is transported vertically out of the disk through fluid elements tied to open magnetic field lines. These vertical winds can be launched if the magnetic field lines exceeds a critical angle of about 30o to the disk normal. Due to such vertical outflows of gas from the disk surface and the conservation of angular momentum the bulk of the disk will be affected by a magnetic braking torque which can generate gas accretion [1]. Turbulence transports angular momentum through the fact that protoplanetary disks are differentially rotating. This means that gas in an outer region will orbits slower than gas closer to the star. This leads to a shear between the two radially separated gas regions causing the faster orbiting gas elements to transport angular momentum outwards and the gas in the outer regions to move radially inwards. If the two gas regions also are weakly magnetized the shear can result in a magne- torotational instability which drives magnetohydrodynamic turbulence [1]. This is discussed in further detail in Section 1.4. In this thesis the main focus will be on planet formation through the core accre- tion theory with angular momentum transport through turbulence. Hence onwards the GI scenario and the disk winds will not be considered in further detail, although one should bare in mind that these processes likely play and important role in planet formation processes. 1.3 The Search for Earth-like Exoplanets With the launch of the Kepler Space Telescope in March 2009 the goal was to observe the first Earth-sized planets around solar-like stars [5]. To detect these types of planets the Kepler mission used a detection method known as the Transit method. This method detects exoplanet candidates by measuring the flux from the host star. In the event of a transiting object, like a planet, the stellar flux detected by the telescope would decrease indicating a potential orbiting planet [20]. 4 1. Introduction The Kepler mission resulted in a lot of new detected planetary candidates, many of which had sizes of super-Earths, with radii of „ 0.5 ´ 2 RC, to sub-Neptunes with radii of „ 2 ´ 4 RC [21]. When studying the data from the Kepler mission it turned out that many planets were found in multiplanetary systems. The inner planets in these multiplanetary systems were found to orbit quite close to each other and hence these systems have been known as Systems with Tightly packed Inner Planets or STIPs [22][23]. By studying the orbits of these systems it was found that the planets showed little statistics of being in resonant orbits [22]. To explain the formation of the STIPs one early theory proposed was that the planets could have formed in the outer disk by core accretion and then migrated inwards. The migration would occur as consequence of the planet’s interaction with the gas in the protoplanetary disk. However, by making simulations of migrating planets from the outer disk it was proven hard to obtain super-Earth sized planets since the accretion of gas and pebbles often resulted in gas giants or ice giants [24]. Planets that managed to migrate fast enough were by simulations shown to end up in resonant orbits [25], which is the opposite from what was found by observations and the Kepler data. Hence there was a need for another theory of how these STIPs could form and end up in non-resonant orbits. One possible solution to this problem was suggested to be that the planetary embryos could in situ undergo gravitational assembly into final planets, hence avoid the need for planetary migration [26]. This theory seemed to quite successfully reproduce planetary systems similar to the ones observed by the Kepler mission [26]. Based on these results a theory of in situ planet formation known as Inside-Out Planet Formation (IOPF) was developed. This IOPF theory will be at the core of this thesis and is explained in detail in Section 1.4. 1.4 Inside-Out Planet Formation As mentioned in the previous section the Inside-Out Planet Formation or IOPF is a theory for in situ formation of planets that could explain the formation of the STIPs. This theory is based on the fact that in the outer part of the protoplanetary disk dust grains can stick together and coagulate into centimeter to meter sized solids, called pebbles [2]. These pebbles will be affected by a drag force exerted on them by the gas in the disk, so-called gas drag. This gas drag is caused by a radial pressure gradient in the gas that makes the gas in the disk rotate slower than the free orbital velocity, hence exerting drag on freely moving pebbles [27]. This drag will cause the pebbles to start moving radially and the radial velocity is strongly dependent on the size of the particles [27]. If the particles grow to meter-size rocks, the gas drag often either causes these rocks to fragment or causes a high radial velocity that makes these rocks drift into the star. This problem is known as the meter-size barrier since its a persistent problem in planet formation theories to grow pebbles beyond the order of a meter which would inhibit the formation of planets [28]. However, even if the meter-size barrier would inhibit planet formation in most parts of the disk one could argue that the drift caused by the gas drag is crucial for enabling planet formation in the inner disk [2]. In the IOPF theory the pebbles are moving radially inwards until their radial motion is halted at the dead-zone 5 1. Introduction r a b c d Figure 1.2: A sketch over how the MRI arises. Two fluid elements are linked vertically by a magnetic field line. In a) the elements are slightly perturbed radially shown by the black arrows. In b) the shear leads to an azimuthal displacement shown by the orange arrows. This leads to magnetic tension, shown by the green arrows in c), through which angular momentum is transferred. This causes the radial displacement of the fluid elements to increase, shown by the black arrows in d), leading to an instability. Inspiration to this sketch is taken from [1]. inner boundary (DZIB). The DZIB is the boundary at which the temperature is high enough for thermal ionization of alkali metals to commence and give rise to the magnetorotational instability (MRI) [29]. This usually occurs when the temperature exceeds „ 1000 K [30]. The MRI is an instability that occurs due to cylindrical shear flows if a weak magnetic field is present. If a disk contains a magnetic field the magnetic field lines can be thought of as oriented vertically through the disk. These vertical field lines can be slightly perturbed in the radial direction linking fluid elements at different radii. Since the fluid element at the inner radii will rotate faster than the element at the outer radii, due to the shear in the disk, a toroidal field component will be obtained in addition to the initial vertical and radial ones. Due to spring-like properties of magnetic field lines the two fluid elements will be affected by azimuthal forces through the magnetic field linking the both elements. The outer element will feel a force in the direction of its motion and speed up whereas the inner element will feel a force in the direction opposite to its motion and slow down. This will decrease the angular momentum of the inner element which will start moving inwards and increase the angular momentum of the outer element which will start moving outwards. This increases the radial displacement between the two elements causing an instability [1]. In Figure 1.2 a schematic sketch can be seen over how the radial perturbation of two fluid elements linked by a magnetic field line give rise to the MRI. The spring-like behaviour of the magnetic field lines causing the MRI can initiate a magnetohydrodynamic (MHD) turbulence in the disk [31]. This MHD turbulence affects the viscosity of the disk and the turbulence and the viscosity can be re- 6 1. Introduction lated through a dimensionless α-parameter. This relation is known as the α-disk model proposed by Shakura and Sunyaev 1973 and this model is further explained mathematically is Section 2.2.3 [32]. Since the viscosity of the disk is dependent on the MHD turbulence which is a consequence of the MRI, the viscosity will be high close to the DZIB where the MRI is active. From the DZIB there will be a fall-off in viscosity moving radially outwards in the disk due to the decreasing thermal ionization and hence a suppression of the MRI [33]. This fall-off in viscosity will cause a local pressure maximum in the gas pressure in the vicinity of the DZIB. At this pressure maximum pebbles of the order of centimeter in size, which have formed in the outer disk and drifted inwards, are being trapped. The trapped pebbles at the local pressure maximum, due to the trapping also known as a pressure trap, quickly coalesce into a protoplanet. The protoplanet itself will also be trapped at this pressure trap and will continue to grow through pebble accretion. The growth of the protoplanet will continue until it has accreted enough matter to become massive enough to clear a gap of a few Hill radii around its orbit in the disk. Matter on the inside of this cleared gap will drain rapidly through accretion onto the star. The cleared gap in combination with this draining of the inner region will decrease the density of this part of the disk significantly compared to the outer dusty part of the disk. This will allow for direct stellar irradiation to reach the outer rim of the cleared gap at which the MRI can be reactivated causing a new pressure trap at which pebbles can coalesce into a second protoplanet. This second protoplanet can then open a new gap allowing for a second retreat of the DZIB allowing for further planets to form. In this way the IOPF theory proposes that the planets are formed from the inside-out in the planetary system by coagulation of pebbles at the pressure traps caused by the MRI. This process proceeds as long as there is a sufficient supply of pebbles from the outer disk. When the pebble supply is exhausted a system of closely packed inner planets will have been formed [33]. In Figure 1.3 a sketch of how planets could form in this manner is being shown. 7 1. Introduction Figure 1.3: A sketch over the IOPF theory. The arrow shows a pebble drift towards the DZIB (top row). Pebbles get trapped at the pressure maximum (second row). At the pressure maximum pebbles coagulate and grow into protoplanets (third row). This clear the region of dust, allowing for stellar irradiation further out in the disk. The DZIB retreat and causes a new pressure maximum at which pebbles get trapped, eventually forming an additional planet (bottom row). This process continues as long as there is a supply of dust and pebbles. Inspiration to this sketch is taken from [2]. 8 1. Introduction 1.5 Thesis Objectives As has been established in the previous sections the dust evolution and the grain growth process in protoplanetary disks are of high importance when trying to un- derstand how planets can form in these disks. Studies of dust evolution and grain growth in protoplanetary disks have been done previously as presented in the pa- per ”Inside-out Planet Formation. IV. Pebble Evolution and Planet Formation Timescales” [34]. In this paper, hereafter Paper IV, the pebble evolution is studied around a solar mass T-Tauri star in a couple of different cases and models. However, the models used in this paper present some limitations regarding taking fragmen- tation and different fragmentation velocities into account. Furthermore, the dust species used in the paper was limited to pebbles and dust. Inspired by Paper IV the goal of this thesis is to study dust evolution in proto- planetary disks in order to better understand the microphysics and the early stages of the planet formation process through core accretion. This will be done by ex- panding one of the models studied in Paper IV, using an α-parameter value of 10´4 and a gas accretion rate of 10´9 Mdyr´1, by implementing the processes for dust growth. Some questions regarding IOPF that this thesis will try to answer are, ”How much of the delivered dust mass to a given region is made up of grain sizes larger than 0.06 cm?”, since these are the grain sizes thought to be trapped at the DZIB, and, ”What time does it take to deliver 3 M‘ to a given region?”, since this is the mass thought to be required to achieve the gap opening in the disk. 9 1. Introduction 10 2 Theory This chapter will provide the necessary theory required for the study of gas and dust evolution in protoplanetary disks through the core accretion and MRI models. Focus will be put on the math behind the models and theories and the necessary equations will be presented, some of which also will be derived. First MHD and MRI will be discussed in both the ideal and non-ideal MHD cases. Then the gas in protoplanetary disks will be discussed where I start discussing the thermal structure, continuing with the dynamical structure, and at last discuss the gas evolution. Then the focus is turned towards the dust where the dust dynamics and evolution is being discussed. This chapter ends with the theory for the grain growth process being presented, with a special focus on the Smoluchowski equation. I would like to point out that the theory presented about the gas and the dust in protoplanetary disks is largely inspired by the theory presented in the PhD thesis by T. N. Delage [35]. 2.1 Magnetohydrodynamics and the MRI As mentioned in the previous chapter the MRI is an instability that arises due to the shear motion between two ionized gas regions at different radii in the protoplanetary disk. Since the gas close to the star where the temperature is higher than 3000 K can be considered strongly ionized the MRI will take the form described by ideal MHD, whereas in the outer regions where the gas only is weakly thermally ionized the MRI will be described by non-ideal MHD [1]. 2.1.1 Ideal MHD To study the MHD one should start by consider the induction equation which can be written as BB Bt “ ∇ ˆ pvgas ˆ Bq ` η∇2B (2.1) where vgas is the gas velocity, B is the magnetic field and η is the magnetic dif- fusivity. Since the ideal MHD is active in the inner region of the disk where the gas is strongly ionized the gas can be considered a perfect conducting fluid and the magnetic diffusivity can hence be set to zero [36]. This leads to the ideal induction equation which can be written as BB Bt “ ∇ ˆ pvgas ˆ Bq . (2.2) 11 2. Theory One can now through the conservation of mass equation and the momentum equation combined with (2.2) find the condition for MRI to be active in the ideal MHD regime to be dΩ2 dr ă 0 (2.3) where Ω is the angular velocity [1]. 2.1.2 Non-ideal MHD In the non-ideal MHD case the magnetic diffusivity can not simply be set to zero, as in the ideal MHD case, since the gas is only weakly ionized and can hence not be considered a perfectly conducting fluid. The induction equation (2.1) can in this case be written as BB Bt “ ∇ˆpvgas ˆ Bq´∇ˆrηO p∇ ˆ Bq ` ηH p∇ ˆ Bq ˆ B ´ ηAD pp∇ ˆ Bq ˆ Bq ˆ Bs (2.4) and this is known as the non-ideal induction equation [37]. In this equation the first term on the right hand side is the ideal MHD contribution and the second, third and fourth term describes the contribution from Ohmic diffusion, the Hall term and the ambipolar diffusion respectively. In this equation ηO, ηH and ηAD are the diffusion coefficients which are defined as ηO “ c2 4πσO , ηH “ c2σH 4π pσ2 H ` σ2 P q , ηAD “ c2σP 4π pσ2 H ` σ2 P q ´ ηO (2.5) where σO, σH and σP are the Ohmic, Hall and Pedersen conductivities and c is the speed of light. The conductivities can be expressed as σO “ ec B ÿ j njZjβj, σH “ ec B ÿ j njZj 1 ` β2 j , σP “ ec B ÿ j njZjβj 1 ` β2 j (2.6) where Zje is the charge and βj is the Hall parameter defined as βj “ ZjeB mjcγjρn (2.7) where ρn is the density of neutral particles. In the Hall parameter γj is defined as γj “ xσvyj mj ` mn (2.8) where xσvyj is the rate coefficient for momentum transfer, mj is the mass of charged particles and mn is the mean mass of neutral particles [37]. By studying the terms defined above one can find that the Ohmic resistivity is mostly active in the inner disk and originates from collisions between free electrons and neutral particles. The ambipolar diffusion becomes important in the outer part of the disk and is caused by collisions between ions and neutrals. In the midplane of the disk the Hall effect becomes important and arises due to non-collisional drifts [38]. 12 2. Theory This non-ideal MHD affect the MRI since the weaker ionization can result in the magnetic field lines to not be sufficiently coupled to the gas motion. If this coupling is weakened it can allow for the magnetic field to diffuse relative to the gas motion and thus inhibit the MRI. This is primarily caused by the Ohmic resistivity and the ambipolar diffusion which in turn can generate the Hall effect which is a new MHD instability [35]. 2.2 Gas in Protoplanetary Disks With the MHD and MRI established we will in this section focus on the theory for the gas in protoplanetary disks. First the thermal structure of the gas will be described after which the focus will be turned towards the dynamical structure of the gas. When this has been done the equations and theory for the gas evolution will be derived and discussed. 2.2.1 Thermal Gas Structure To find the thermal structure of a protoplanetary disk one can start by consider the ideal gas law pgaspt, r, zqV “ nkBT pt, rq (2.9) where pgaspt, r, zq is the thermal gas pressure, V is the volume, n is the number of gas particles, kB is the Boltzmann constant and T pt, rq is the gas temperature. The t, r and z are parameters indicating time, radial and vertical dependence respectively. The number of gas particles can be expressed as n “ ρgaspt, r, zqV mn (2.10) where ρgas is the gas mass density and mn is the mean molecular mass of neutral gas particles. By inserting this expression for n into (2.9) one get the thermal gas pressure as pgaspt, r, zq “ ρgaspt, r, zqkBT pt, rq mn (2.11) where mn “ µmH and where µ is the mean molecular weight and mH being the atomic mass of hydrogen. The pressure can also be expressed in terms of the adia- batic sound speed cs as pgaspt, r, zq “ ρgaspt, r, zq c2 spt, rq γ (2.12) where cspt, rq “ d γkBT pt, rq mn (2.13) and where γ is the ratio of specific heat and γ “ 1 if the gas is isothermal. 13 2. Theory 2.2.2 Dynamical Gas Structure To study the dynamical structure of the gas in a protoplanetary disk one can start by consider the conservation of momentum equation which can be written as ρgas „ Bvgas Bt ` pvgas ¨ ∇q vgas ȷ “ ´∇ ˆ pgas ` B2 8π ˙ ´ ρgas∇Φ ` 1 4π pB ¨ ∇q B (2.14) where ρgas is the gas mass density, vgas is the gas velocity, pgas is the thermal gas pressure, the term B2{8π is the magnetic pressure, Φ is the gravitational potential and B is the magnetic field [35]. Now, the shear stress tensor, W, can be introduced which takes the form W “ ρgasvgas ¨ vgas ´ B ¨ B 4π “ ¨ ˝ Wrr Wrϕ Wrz Wϕr Wϕϕ Wϕz Wzr Wzϕ Wzz ˛ ‚ (2.15) and with the shear stress tensor (2.14) can be rearranged and written as B Bt pρgasvgasq ` ∇ ¨ W “ ´∇ ˆ pgas ` B2 8π ˙ ´ ρgas∇Φ. (2.16) By studying the shear stress tensor one can see that the tensor is symmetric and that the elements in the tensor, in cylindrical coordinates, would be Wrr “ ρgasδvr,gasδvr,gas ´ BrBr 4π Wrϕ “ Wϕr “ ρgasδvr,gasprΩ ` δvϕ,gasq ´ BrBϕ 4π Wrz “ Wzr “ ρgasδvr,gasδvz,gas ´ BrBz 4π Wϕϕ “ ρgas prΩq 2 ` 2ρgasrΩδvϕ,gas ` ρgasδvϕ,gasδvϕ,gas ´ BϕBϕ 4π Wϕz “ Wzϕ “ ρgasrΩδvz,gas ` ρgasδvϕ,gasδvz,gas ´ BϕBz 4π Wzz “ ρgasδvz,gasδvz,gas ´ BzBz 4π where the δ indicates that the velocity is the fluctuating velocity caused by the MRI driven turbulence. The total gas velocity components in the three directions are vr,gas “ δvr,gas, vϕ,gas “ prΩ ` δvϕ,gasq and vz,gas “ δvz,gas. Worth noticing here is that the azimuathal velocity component vϕ,gas has two terms, rΩ which accounts for the mean gas velocity in the azimuthal direction and δvϕ,gas which accounts for the fluctuations. Now the radial direction of (2.16) can be studied and it can in cylindrical coor- dinates be written as B Bt pρgasvr,gasq` „ 1 r B Br prWrrq ` 1 r BWrϕ Bϕ ` BWrz Bz ´ Wϕϕ r ȷ “ ´ B Br ˆ pgas ` B2 8π ˙ ´ρgas BΦ Br (2.17) 14 2. Theory where the term Wϕϕ is kept since it has a significant radial dependence. By assuming the disk being axisymmetric all derivatives with respect to ϕ becomes zero and (2.17) simplifies to B Bt pρgasvr,gasq ` „ 1 r B Br prWrrq ` BWrz Bz ´ Wϕϕ r ȷ “ ´ B Br ˆ pgas ` B2 8π ˙ ´ ρgas BΦ Br . (2.18) Due to different time scales of radial dynamical equilibrium and accretion, where the disk reaches a radial dynamical equilibrium on much shorter timescales then the timescale associated with the accretion process, (2.18) can be divided into these two parts where the part responsible for the accretion are the terms containing the MRI driven fluctuations [35]. From this one can note that only the term ρgas prΩq 2 will contribute to the dynamical equilibrium in the radial direction. This simplifies (2.18) to ´ρgas prΩq 2 r “ ´ B Br ˆ pgas ` B2 8π ˙ ´ ρgas BΦ Br . (2.19) Now the gravitational potential Φ can be studied in more detail. The gravitational potential is given by Φpr, zq “ ´ GM‹ pr2 ` z2q 1{2 (2.20) where G is the gravitational constant, M‹ is the mass of the central star and r and z are the radial and vertical coordinates respectively. By inserting this expression for the gravitational potential into (2.19) and by studying the first order approximation, at which B2 8π « 0, the following is obtained ´ρgasrΩ2 “ ´ Bpgas Br ´ ρgas B Br ˜ ´ GM‹ pr2 ` z2q 1{2 ¸ . (2.21) By dividing all terms by ´ρgas, evaluating the derivative of the gravitational poten- tial term and by rearranging as well as assuming the protoplanetary disk being thin (z ! r) gives rΩ2 “ 1 ρgas Bpgas Br ` GM‹ r2 . (2.22) This equation can be rewritten as follow rΩ2 “ 1 ρgas Bpgas Br ` rΩ2 K (2.23) where ΩK is the Keplerian angular frequency defined as ΩK “ c GM‹ r3 . (2.24) By introducing the pressure gradient parameter ηgas “ ´ 1 2rΩ2 Kρgas Bpgas Br , (2.25) 15 2. Theory which should not be confused with the magnetic diffusivities, one can use this pa- rameter in (2.23), rearrange the terms and write it as rΩ2 “ rΩ2 K p1 ´ 2ηgasq . (2.26) By taking the square root of both the right hand side and the left hand side of this equation and multiply both sides by ? r the mean gas velocity in the azimuthal direction, without fluctuations, vϕ,gas,mean, is obtained as vϕ,gas,mean “ rΩ “ rΩK p1 ´ 2ηgasq 1{2 (2.27) and by defining the Keplerian velocity vK “ rΩK , the mean azimuthal gas velocity, without fluctuations, becomes vϕ,gas,mean “ vK p1 ´ 2ηgasq 1{2 . (2.28) Now, when the dynamical structure of the gas in the radial direction has been established, lets turn the focus onto the dynamical structure of the gas in the vertical direction. Once again, lets start by considering (2.16), but this time only the terms affecting the vertical direction. This equation can in the vertical direction be written as B Bt pρgasvz,gasq ` „ 1 r B Br prWrzq ` 1 r BWϕz Bϕ ` BWzz Bz ȷ “ ´ B Bz ˆ pgas ` B2 8π ˙ ´ ρgas BΦ Bz (2.29) and by once again assuming axisymmetry the terms containing derivatives with respect to ϕ becomes zero and (2.29) simplifies to B Bt pρgasvz,gasq ` „ 1 r B Br prWrzq ` BWzz Bz ȷ “ ´ B Bz ˆ pgas ` B2 8π ˙ ´ ρgas BΦ Bz . (2.30) One can now note that none of the terms in Wrz or Wzz are contributing to the vertical dynamical equilibrium, only terms responsible for accretion are present, and hence (2.30) in the case of vertical dynamical equilibrium becomes 0 “ ´ B Bz ˆ pgas ` B2 8π ˙ ´ ρgas BΦ Bz . (2.31) Now, a series of assumptions can be done. First, lets assume the first order approx- imation again, B2 8π « 0, which reduces the equation above to Bpgas Bz “ ´ρgas BΦ Bz (2.32) and by now taking the derivative BΦ Bz and once again assume a thin disk (z ! r) gives Bpgas Bz “ ´ρgasΩ2 Kz. (2.33) By using the ideal gas law from (2.12) ρgas can be expressed in terms of pgas and (2.33) can be written as Bpgas Bz “ ´pgas Ω2 K c2 s z. (2.34) 16 2. Theory This equation can easier be written by introducing the gas vertical scale height, defined as Hgas “ cs ΩK , (2.35) and by inserting this into (2.34) it can be written as Bpgas Bz “ ´ pgas H2 gas z. (2.36) This equation can now be rearranged and vertically integrated as ż dpgas pgas “ ż ´ 1 H2 gas zdz ùñ pgaspt, r, zq “ pgaspt, r, z “ 0qexp ˆ ´ z2 2H2 gaspt, rq ˙ (2.37) and since pgas and ρgas are related through the ideal gas law the mass density can also be expressed in a similar way as ρgaspt, r, zq “ ρgaspt, r, z “ 0qexp ˆ ´ z2 2H2 gaspt, rq ˙ . (2.38) Now, the gas surface density can be introduced as Σgaspt, rq “ ż `8 ´8 ρgaspt, r, zqdz “ ż `8 ´8 ρgaspt, r, z “ 0qexp ˆ ´ z2 2H2 gaspt, rq ˙ dz (2.39) and by evaluating this integral one find Σgaspt, rq “ ρgaspt, r, z “ 0q ? 2πHgaspt, rq ùñ ρgaspt, r, z “ 0q “ Σgaspt, rq ? 2πHgaspt, rq . (2.40) By inserting this expression for ρgaspt, r, z “ 0q into (2.38) one can finally express the gas mass density as ρgas pt, r, zq “ Σgaspt, rq ? 2πHgaspt, rq exp ˆ ´ z2 2H2 gaspt, rq ˙ . (2.41) 2.2.3 Gas Evolution Now, when both the radial and the vertical structure of the gas in the protoplanetary disk have been established we are in a position to study the evolution of the gas in the disk. Once again lets consider (2.16), but this time only in the azimuthal direction. One then obtain B Bt pρgasvϕ,gas,meanq` „ 1 r B Br prWrϕq ` 1 r BWϕϕ Bϕ ` BWϕz Bz ȷ “ ´ B Bϕ ˆ pgas ` B2 8π ˙ ´ρgas BΦ Bϕ (2.42) and by assuming axisymmetry, once again making all the derivatives with respect to ϕ zero and by multiplying all the remaining terms by r (2.42) becomes B Bt prρgasvϕ,gas,meanq ` „ 1 r B Br ` r2Wrϕ ˘ ` r BWϕz Bz ȷ “ 0. (2.43) 17 2. Theory By now vertically integrating this equation and remember that vϕ,gas,mean “ rΩ it can be written as B Bt ` r2ΩΣgas ˘ ` „ 1 r B Br ˆ r3ΩΣgasδvr,gas ` r2 ˆ Σgasδvr,gasδvϕ,gas ´ BrBϕ 4π ˙˙ȷ “ 0. (2.44) where the bar indicated that these parameters have been vertically integrated. Here it has also been assumed that the angular momentum only is transported in the radial direction through the disk and the Wϕz-term has hence been discarded [35]. By defining the shear stress tensor in the rϕ-plane as Trϕ “ ρgasδvr,gasδvϕ,gas ´ BrBϕ 4π (2.45) (2.44) can be written as B Bt ` r2ΩΣgas ˘ ` „ 1 r B Br ˆ r3ΩΣgasδvr,gas ` r2 ż `8 ´8 Trϕdz ˙ȷ “ 0. (2.46) Now, lets introduce the continuity equation which in fluid dynamics describes the conservation of mass in a system and in cylindrical coordinates can be expressed as Bρgas Bt ` 1 r B Br prρgasδvr,gasq ` B Bz pρgasδvz,gasq “ 0 (2.47) where once again axisymmetry has been assumed making the derivative with respect to ϕ equal to zero. Vertically integrating this equation and assuming a thin disk (z ! r) gives BΣgas Bt ` 1 r B Br ` rΣgasδvr,gas ˘ “ 0. (2.48) The two equations (2.46) and (2.48) can now be combined. This is done by first expanding the derivatives in (2.46) with the product rule. One then obtain r2ΩBΣgas Bt `Σ B Bt ` r2Ω ˘ ` „ 1 r r2ΩB B ` rΣgasδvr,gas ˘ ` 1 r rΣgasδvr,gas B Br ` r2Ω ˘ ` 1 r B Br ˆ r2 ż `8 ´8 Trϕdz ˙ȷ “ 0 (2.49) and by using (2.48) in the first term in the equation above it will cancel the third term in the equation above giving Σgas B Bt ` r2Ω ˘ ` „ 1 r rΣgasδvr,gas B Br ` r2Ω ˘ ` 1 r B Br ˆ r2 ż `8 ´8 Trϕdz ˙ȷ “ 0. (2.50) It can also be noticed that the term r2Ω does not have any time dependence, mak- ing its derivative with respect to time zero, further simplifying (2.50) which after rearrangement of terms becomes ´rΣgasδvr,gas B Br ` r2Ω ˘ “ B Br ˆ r2 ż 8 ´8 Trϕdz ˙ . (2.51) 18 2. Theory By now assuming that for the protoplanetary disk Ω « ΩK and remembering that ΩK has an r-dependence one can expand the derivative on the left hand side of the equation above giving ´rΣgasδvr,gas 1 2rΩ “ B Br ˆ r2 ż 8 ´8 Trϕdz ˙ (2.52) which after rearrangement of factors and by multiplying both the right hand side and the left hand side of the equation with 2π can be written as 9Mgas “ ´2πrΣgasδvr,gas “ 4π rΩ B Br ˆ r2 ż `8 ´8 Trϕdz ˙ (2.53) and this expression can be identified as the gas accretion rate denoted with 9Mgas. This gas accretion rate can also be written as 9Mgas “ 4π ? r B Br ˆ? r Ω ż `8 ´8 Trϕdz ˙ (2.54) and that this expression is equal to the one presented in (2.53) is shown in Appendix A.1. Now, by putting the expressions for the gas accretion rate presented in (2.53) and (2.54) equal to each other and dividing by ´2π gives rΣgasδvr,gas “ ´2 ? r B Br ˆ? r Ω ż `8 ´8 Trϕdz ˙ . (2.55) This expression can now be inserted into (2.48) giving BΣgas Bt ` 1 r B Br ˆ ´2 ? r B Br ˆ? r Ω ż `8 ´8 Trϕdz ˙˙ “ 0 (2.56) and by rearranging one get BΣgas Bt “ 2 r B Br ˆ ? r B Br ˆ? r Ω ż `8 ´8 Trϕdz ˙˙ (2.57) Now, the shear stress tensor in the rϕ-plane Trϕ can be parameterized as Trϕ “ 3 2αpgas “ 3 2αρgasc 2 s (2.58) where α is a parameter taking values 0 ă α ă 1 [39]. Integrating this relation vertically gives ż `8 ´8 Trϕdz “ 3 2c2 s ż `8 ´8 αρgasdz “ 3 2αc2 sΣgas (2.59) where α indicates that the α-parameter has been vertically integrated and is called the effective disk viscosity parameter. This expression can now substitute the inte- gral in (2.57) and the equation BΣgas Bt “ 2 r B Br ˆ ? r B Br ˆ? r Ω 3 2αc2 sΣgas ˙˙ (2.60) 19 2. Theory is obtained. By now introducing the vertically integrated kinematic disk viscosity ν defined as ν “ αcsHgas “ α c2 s Ω , (2.61) where it has been assumed Ω « ΩK , equation (2.60) can be written as BΣgas Bt “ 3 r B Br ˆ ? r B Br `? rΣgasν ˘ ˙ . (2.62) This is the evolution equation for the gas surface density of a viscously accreting protoplanetary disk. Finally (2.53) and (2.54) can be combined together with (2.61) and (2.58) to obtain an expression for the radial gas velocity which becomes δvr,gas “ ´ 9Mgas 2πrΣgas “ ´ 3 ? rΣgas B Br `? rΣgasν ˘ . (2.63) How this expression is obtained is explicitly shown in Appendix A.2. 2.3 Dust in Protoplanetary Disks With the physics of the gas in the protoplanetary disks now established we change the focus towards the dust. In this section the relevant equations and theory for the dust dynamics and evolution is being presented and discussed. Finally the grain growth process is also being discussed with a special focus on collisional outcomes and the Smoluchowski equation. 2.3.1 Dust Dynamics To describe the dynamics of the dust in protoplanetary disks one can start by as- suming that the dust grains are spherical with radius a and a solid state bulk density ρs. The mass of a dust grain with radius a can then be defined as mpaq “ 4πa3 3 ρs (2.64) and in the same way, if the mass is known, the grain size becomes a “ ˆ 3m 4πρs ˙1{3 . (2.65) In order to now study the dynamics of the dust a drag force density can be introduced that describes the collisions between gas and dust particles. This drag force density can be written as fdrag “ ρdust ˆ vdust ´ vgas τs ˙ (2.66) where ρdust is the dust mass density, vdust is dust particle velocity, vgas is the gas flow velocity and τs is the dust particle stopping time. Now Stokes number can 20 2. Theory be introduced, which is a dimensionless quantity describing the ratio between the dynamical and the stopping timescales. Stokes number can be defined as St “ τsΩK (2.67) and in general a small Stokes number means that the dust particles are strongly coupled with the gas flow whereas a large Stokes number indicates that the dust is highly uncoupled from the gas. The Stokes number can be divided into different regimes, two of which are the Epstein regime and the Stokes I regime. For these two different regimes the vertically integrated Stokes number can be written as Stpt, r, aq “ # π 2 aρs Σgaspt,rq , a ă 9 4λmfp 2π 9 a2ρs Σgaspt,rqλmfp , a ą 9 4λmfp (2.68) where the first regime, a ă 9 4λmfp, is the Epstein regime and the second one, a ą 9 4λmfp, is the Stokes I regime [40]. The λmfp factor is the mean free path in the midplane for the gas particles. Since the dust particles orbit faster than the gas in the azimuthal direction the individual dust particles will feel a headwind resulting in a loss of angular momen- tum. This loss of angular momentum will cause the dust particles to radially start drifting inwards. The radial drift velocity of dust particles of size a can be found to be vr,dustpt, r, aq “ δvr,gaspt, rq 1 ` St2 pt, r, aq ´ 2ηgaspt, r, z “ 0qvKprqStpt, r, aq 1 ` St2 pt, r, aq (2.69) where the first term describes the radial gas flow and the second term describes the radial drift [41]. The azimuthal dust velocity can be written as vϕ,dust “ vKprq ´ ηgaspt, r, z “ 0qvK 1 ` St2 pt, r, aq (2.70) where the first term is the Keplerian velocity and the second term is the azimuthal drift. In the vertical direction it can be found that the vertical drift velocity or the sedimentation velocity of the dust towards the midplane can be expressed as [42] vz,dust “ ´zΩKmin ˆ Stpt, aq, 1 2 ˙ . (2.71) Now, to avoid getting into the details of the vertical dust distribution for each dust grain, z can be replaced by the vertical dust scale height which can be taken as an average height above the midplane for each species a and can be expressed as Hdustpt, r, aq “ Hgaspt, rq d δzpt, rq δzpt, rq ` Stpt, r, aq (2.72) where δz is the vertical mixing parameter that controls the mixing between the gas and the dust in the vertical direction. Since the mixing of gas and dust is caused by the turbulence in the disk the vertical mixing parameter can be assumed to be 21 2. Theory δz “ α, where α is the effective disk viscosity parameter. Furthermore, the turbulent motion of the gas will counteract local concentration gradients of the dust. The gas hence make the dust diffuse and the diffusion of a species with size a can in the radial direction be expressed as Ddustpt, r, aq “ δrpt, rqc2 spt, rq ΩKprq ` 1 ` St2 pt, r, aq ˘ (2.73) where δr is the radial mixing parameter [43]. Since the mixing is a consequence of the turbulence the mixing parameter can be assumed to be δr “ α. By using the definition of the effective kinematic disk viscosity ν, presented in (2.61), the expression for the diffusion can be rewritten as Ddustpt, r, aq “ ν ` 1 ` St2 pt, r, aq ˘ . (2.74) With this foundation the radial transport of the dust can now be studied. The radial transport can be described by a diffusion-advection equation which in cylin- drical coordinates takes the form BΣdustpt, r, aq Bt ` 1 r B Br prFtotpt, r, aqq “ 0 (2.75) where Ftot is the total flux and can be expressed as Ftot “ Fadv ` Fdiff where Fadv and Fdiff are the advection and diffusion fluxes respectively [42]. The advection flux contribution can be expressed as Fadvpt, r, aq “ Σdustpt, r, aqvr,dustpt, r, aq (2.76) and the diffusion flux contribution can be expressed as [42] Fdiff pt, r, aq “ ´Σgaspt, rqDdustpt, r, aq B Br ˆ Σdustpt, r, aq Σgaspt, rq ˙ . (2.77) Combining (2.76) and (2.77) and substituting them into the position of Ftot in (2.75) the final equation for the radial dust transport becomes BΣdust Bt ` 1 r B Br ˆ r ˆ Σdustvr,dust ´ ΣgasDdust B Br ˆ Σdust Σgas ˙˙˙ “ 0. (2.78) What remains to consider in the dust dynamics is now the vertical structure. In a similar way as the gas mass density, shown in (2.41), the dust mass density for each species a can be expressed as ρdustpt, r, z, aq “ Σdustpt, r, aq ? 2πHdustpt, r, aq exp ˆ ´ z2 2H2 dustpt, r, aq ˙ (2.79) and the total dust mass density can be obtained by summing over each species as ρdust,totpt, r, zq “ ÿ a ρdustpt, r, z, aq “ ÿ a Σdustpt, r, aq ? 2πHdustpt, r, aq exp ˆ ´ z2 2H2 dustpt, r, aq ˙ . (2.80) 22 2. Theory Now, with the dust mass density established the number density for each grain size a can be obtained by ndustpt, r, z, aq “ ρdustpt, r, z, aq mpaq (2.81) where mpaq is the mass of the grain species with radius a. The total number density for the dust can finally be expressed by summing over all the dust species a as ndust,totpt, r, zq “ ÿ a ndustpt, r, z, aq “ ÿ a ρdustpt, r, z, aq mpaq . (2.82) 2.3.2 Grain Growth Process The planet formation process is heavily dependent on grain growth, since the build- ing blocks of planets are dust grains that can coagulate into pebbles, boulders and eventually even planetesimals. Due to this fact the grain growth process is of great importance to understand and study when considering the formation of planets. The grain growth process is based on the fact that dust grains of different sizes and relative velocities collide with each other in the protoplanetary disks. The rate of the grain growth is dependent on the outcome of each collision and the outcome is dependent on the grain sizes involved in a certain collision and the relative velocities of the species. As mentioned the grain sizes and their relative velocities can result in a number of different outcomes from a collision. If the relative velocity is sufficiently low the particles might stick together after the collision due to the electrostatic van der Waals forces resulting in one final particle with higher mass[44][45]. This outcome is what is called sticking or coagulation. If the relative velocity between the two species in a collision is higher than a certain velocity the two species will break each other apart into smaller fragments. This threshold velocity is what is known as the fragmentation velocity and this collisional outcome is known as fragmentation [45]. Other possible outcomes from a collision is bouncing where the two grains collide with each other elastically and no change in mass for the two grains is achieved [45]. Furthermore, another collisional outcome is erosion, which occurs when one of the colliding species is much larger than the other colliding species and when their relative velocity is high. The smaller species might shatter in the collision, but might not be large enough to completely fragment the larger species, rather break loose smaller fragments of the larger body. This outcome is called erosion [45]. To simulate the grain growth and taking the different collisional outcomes into consideration one can use the Smoluchowski equation which can be written as Bnpmq Bt “ ż `8 0 ż m1 0 Kpm, m1, m2 qnpm1 qnpm2 qRpm1, m2 qdm2dm1 ´ npmq ż `8 0 npm1 qRpm, m1 qdm1. (2.83) This equation describes the evolution of the number distribution npmq of dust grains with mass m [40]. In this equation the first term on the right hand side describes the supply of newly formed grains of the type npmq. In this term Kpm, m1, m2q is 23 2. Theory the collision Kernal tensor that holds information about the collisional outcome and describes the amount that gets added to npmq from collisions between the species m1 and m2, npm1q and npm2q is the number distributions of the species m1 and m2 respectively and Rpm1, m2q describes the reaction rate between the species m1 and m2. It can be noted that the second integral only is evaluated to m1 and the reason for this is to avoid counting the reactions between m1 and m2 twice, since this is the same reaction as m2 colliding with m1. The second term on the right hand side of this equation describes the amount of particles in the npmq distribution that gets removed as a consequence of collisions with other particles. The continuous distribution npmq can be discretized by creating mass bins as nk “ ż k` 1 2 k´ 1 2 npmqdm (2.84) where k is the index of the mass bin and indicate the bin center [40]. The discretized version of the Smoluchowski equation can then be written as Bnk Bt “ Nm ÿ l“1 l ÿ p“1 KklpRlpnlnp ´ nk Nm ÿ q“1 nqRqk p1 ` δqkq (2.85) where Nm is the number of bins and δqk is the Kronecker delta [40]. 24 3 Methods In this chapter the methods used to simulate the dust evolution is explained, where the focus first is put on the program used, DustPy, and how the simulations were set up in this program. Here the parameters and boundary conditions used are also being discussed. Finally the methods used for the data analysis are briefly being discussed. 3.1 DustPy To simulate the evolution and grain growth of the dust the Python-based program DustPy was used, developed by S. Stammler and T. Birnstiel [40]. DustPy simulates the evolution of gas and dust by solving the relevant gas and dust equations for a specified number of time steps. The gas evolution is in DustPy simulated by solving the evolution equation for the gas surface density (2.62). However, since the focus of this thesis is to study the dust evolution and grain growth, the gas evolution was in DustPy turned off. Hence the background gas disk and temperature was set to be constant over time in the simulations and the gas surface density profile and temperature profile used is described and shown in Section 3.1.1.3. DustPy simulated the dust evolution by solving the radial dust transport equation (2.78) and the grain growth was simulated by solving the Smoluchowski equation. The collisional outcomes of colliding grains used in this thesis were sticking, frag- mentation and erosion, which are the three default outcomes in DustPy. However, the colliding particles do in DustPy not collide with a single relative velocity, but rather the relative velocities are following the Maxwell-Boltzmann distribution. The fragmenting collision rate is then given by integrating over all the possible relative velocities above the fragmentation velocity in the Maxwell-Boltzmann distribution [40]. It is at this point also important to point out that DustPy only has one spatial dimension, the radial dimension, which also is the distance from the star r. This is the reason why most of the equations in Section 2 has been integrated vertically over the distance z between the dust particles and the disk midplane. However, the settling of dust grains towards the midplane is accounted for in DustPy through the dust scale heights in the relative collisional velocities between different dust species [40]. Now, with this general background and knowledge about DustPy we are in a position to describe how the simulations done in this thesis were set up and carried out. This is in detail described in the following Section 3.1.1 where the setup, initial 25 3. Methods profiles and parameters are presented. 3.1.1 Simulation Setup and Parameters As mentioned in the previous section DustPy simulated the evolution of gas and dust in the protoplanetary disks by solving the evolution equation for the gas surface density and the radial dust transport respectively for a specified number of time steps. In this thesis the simulations were run between 0 ´ 5 ˆ 106 yr with a total of 96 snapshots. These snapshots were equally spaced within each decade, up to 105 yr after which the spacing between each snapshot was 105 yr. The time range of the simulations and the snapshot spacing is summarized in Table 3.1. Table 3.1: The simulation snapshots and spacing summarized. Range (yr) Snapshot spacing (yr) 0 ´ 101 100 101 ´ 102 101 102 ´ 103 102 103 ´ 104 103 104 ´ 105 104 105 ´ 5 ˆ 106 105 3.1.1.1 Grid For the simulations carried out in this thesis a radial grid with 128 grid cells between 0.6 au and 30 au has been used. The mass grid used consisted of 162 bins between 6.99 ˆ 10´15 g and 4.48 ˆ 108 g which corresponds to grain sizes, assuming them being spherical and with a bulk density of 1.67 gcm´3, between 0.1 µm and 4 m. The number of mass bins per decade used was set to 7, since this is the minimum number of bins to avoid an over estimation of the growth rate [46]. 3.1.1.2 Star To study the dust evolution in the protoplanetary disk the choice of star is impor- tant. In this thesis a solar-mass T-Tauri star of radius R‹ “ 3.0Rd and effective temperature T‹ “ 4500 K was chosen. These initial parameters are summarized in Table 3.2. Table 3.2: The initial stellar parameters used in the simulations. Parameter Value stellar mass (M‹) 1.0Md stellar radius (R‹) 3.0Rd stellar temperature (T‹) 4500 K 26 3. Methods Figure 3.1: The initial gas profile where the gas surface density is shown in the top row and the initial gas temperature profile is shown in the bottom row. 3.1.1.3 Gas As previously mentioned, for the simulations done in this thesis the gas evolution was turned off, and the gas considered was only a background gas disk with a fixed surface density and temperature profile. These profiles were set to be the same as the profiles used in Paper IV and these gas profiles are shown in Figure 3.1. The profiles used are based on a viscosity parameter of α “ 10´4 and a gas accretion rate of 9Mgas “ 10´9 Mdyr´1. Besides these fixed profiles some additional gas parameters had to be specified in DustPy. These parameters were all also set to the same values as used in Paper IV and these parameters and corresponding values are shown in Table 3.3. Table 3.3: The initial gas parameters used in the simulations. Parameter Value viscosity parameter (α) 10´4 specific heat ratio (γ) 1.4 mean molecular weight (µ) 2.33mH initial gas disk mass 0.03M‹ 27 3. Methods 3.1.1.4 Dust To simulate the dust a couple of initial parameters needed to be specified. The gas to dust ratio was set to 0.01. The initial grain size distribution in DsutPy followed a so-called MRN-distribution up to a maximum initial particle size which in the simulations was set to 10´4 cm. Since the MRN-distribution takes the form npaq9adistExp (3.1) where distExp is the distribution exponent which also is a parameter and was spec- ified to ´3.5 [47]. Another parameter that needed to be specified was the bulk density or monomer density of the dust particles which was set to 1.67 gcm´3. Since one of the objectives of this thesis was to investigate how different fragmen- tation velocities affected the dust evolution three simulations with three different fragmentation velocities were carried out. These fragmentation velocities were spec- ified to be 1 ms´1, 5 ms´1 and 10 ms´1. The erosion mass ratio was set to the default value of 10, meaning that if the mass ratio between two colliding particles was less than this value, both of them fully fragmented. If the mass ratio was above this value, erosion would occur. In an erosive event the mass excavated from the heavier body was specified to the default value of 1 in units of the lighter particle. In the case of a fragment collision the fragments would have a mass distribution of npmqdm9mfragmentDistributiondm (3.2) where the default value ´11{6 of the fragmentDistribution exponent was used taken from [48]. Furthermore, another parameter that could be specified was to initially allow drifting particles or not, and this parameter was set to false. The reason for this was to avoid a visible particle wave of small particles from the outer disk at time zero in the simulation, since these otherwise could be drifting at the initialization due to a high Stokes number in a part of the disk where the gas surface density was low. All of the used dust parameter values are summarized in Table 3.4. Table 3.4: The initial dust parameters used in the simulations. Parameter Value dust-to-gas ratio 0.01 initial maximum grain size 10´4 cm exp of MRN distribution ´3.5 monomer density (ρs) 1.67 gcm´3 fragmentation velocity 1, 5, 10 ms´1 erosion mass ratio 10 excavated mass 1 exp of fragment distribution ´11{6 allow drifting particles False 28 3. Methods 3.1.2 Simulation Boundary Conditions For each simulation carried out in this thesis inner and outer boundary conditions for the gas and dust had to be specified. The inner boundary for the gas was set to a fixed value and this value was chosen to be the inner boundary value for the Σgas profile shown in Figure 3.1, which equals to 1264 gcm´2. The outer gas boundary was also set to a fixed value, but the value chosen was a floor value of 10´100 gcm´2 to prevent the dust-to-gas ratio to become unrealistically high and to prevent inflow of gas through this boundary. For the dust the default boundary conditions were used. The inner boundary was set to a constant gradient while the outer boundary was set to a constant floor value [40]. 3.2 Data Analysis In order to study and analyze the data returned from the simulations different plots were created. To study the evolution of the grain sizes probability density functions (pdf) were created over the grain sizes. The evolution of the grain sizes were then studied by studying the evolution of these pdfs for the different fragmentation ve- locity cases and at different regions of the disk, where the radii 1 au, 2 au, 5 au and 10 au were chosen to focus on. The pdfs were created by weighting each mass bin with the total dust mass as fpaq “ 2πrAprqρdustpt, r, aq ř a 2πrAprqρdustpt, r, aq “ ρdustpt, r, aq ř a ρdustpt, r, aq (3.3) where 2πr is the circumference of a torus with the torus center at a distance r from the star, Aprq is the surface area of an annulus a grid cell spans with the grid center at r and ρdust is the dust mass volume density. So, by making the calculation 2πrAprq the volume of a torus centered at a distance r from the star is obtained. By then multiplying this volume with the dust mass volume density the mass made up of a certain grain size at a given time and radius is obtained. By then in the denominator summing over all the grain sizes, and hence also masses, the mass weighted pdf fpaq was obtained. 2D-plots were also created to study the dust density evolution. The dust density was plotted and color coded on a mesh with the distance from the star on the x-axis and the grain radius on the y-axis. The density plotted was σd which is related to the dust surface density as Σdustpt, r, aq “ ż 8 0 σdpt, r, a, mqd logm. (3.4) This calculation makes the dust surface density distribution independent of the mass grid. However, since the properties that were being studied was grain sizes and masses, σd was used in the 2D-plots. In these plots lines for the analytical estimates for the fragmentation barrier and drift barrier as well as a line indicating St=1 were plotted. This was done in order to study if the fragmentation or drift was the most limiting factor when studying grain sizes and grain growth in the different 29 3. Methods cases. The analytical estimates for the fragmentation and drift barriers were taken from Birnstiel et al. [49]. To study the dust flux through the 1 au radius the dust accretion rate was cal- culated as 9Mdustpt, r, aq “ ´2πrΣdustpt, r, aqvr,dustpt, r, aq. (3.5) The mass delivered through a given radius as 1 au could then be calculated as a function of time t and dust species a as Mdustpt, aq ż t 0 9Mdustpt 1, r “ 1 au, aqdt1. (3.6) The total mass delivered by grains larger than a certain grain size was then obtained by summing up the delivered mass of grain sizes above the desired threshold as Mdustptq “ ÿ a Mdustpt, aq. (3.7) 30 4 Results & Discussion In this chapter the results from the simulations carried out in this thesis is being presented and discussed. First the grain size distributions are considered for the different fragmentation velocity cases at different times and at different parts of the disk. Then the focus is turned onto the dust mass flux and dust mass delivery to the inner parts of the disk. Finally, the results found when investigating the radial resolution’s affect on the grain size distributions are being presented and discussed. 4.1 Grain Size Distributions In order to study the grain size distribution probability density functions over the grain sizes were created as described in Section 3.2. The pdfs for the fragmentation velocity 1 ms´1 is shown in Figure 4.1. For low fragmentation velocities as 1 ms´1 the dust particles tend to easily frag- ment. This results in many smaller particles since the low fragmentation velocity prevents them from sticking and growing into larger pebbles. The smaller grain sizes means less gas drag from the gas disk and hence also a less effective radial drift. This enables for the disk to reach a stable distribution of grain sizes quickly, especially in the inner disk at small radii as 1 au due to evolution on dynamical time scales. It can from Figure 4.1 be seen that the grain sizes have reached a stable distribution after about 0.1 Myr at which the representative grain size is 0.59 cm at 1 au. After this point the position of the peak of the distribution at 1 au stays constant over time and the dominant grain size after 5 Myr remains 0.59 cm. Considering the outer parts of the disk the grain size distribution seem to be less stable over time. The evolution in this part of the disk is slower due to the dynamical time scale and the gas drag allows for the larger particles to drift inwards replenishing the inner regions and depleting the outer regions of larger grain sizes. This causes the representative grain size at an outer radius of 10 au to drift towards smaller grain sizes with time and this is what can be seen in Figure 4.1, where the dominant grain size at 0.1 Myr is 0.27 cm compared to 0.012 cm after 5 Myr. That the grain sizes in general are smaller in this outer region of the disk can also be explained by the lower dust density in this outer part of the disk, which can be seen from Figure 4.2. In Figure 4.2 it can be seen how the dust density, shown by the color coding, continuously decreases in the outer region. In the inner region the dust density seem to be more constant over time until very late times at which the density decreases. The constant density in the inner region might be explained by a replenishing of 31 4. Results & Discussion Figure 4.1: Probability density functions based on the simulation with the frag- mentation velocity of 1 ms´1. The pdfs are shown at two different radii and at each radii five times are shown between 0.1 Myr and 5 Myr. The peaks represent the most probable grain size at that given time and radius. 32 4. Results & Discussion Figure 4.2: 2D plots showing the density evolution for given grain sizes and given radii at four times between 0.1 Myr and 5 Myr for the fragmentation velocity 1 ms´1. The blue line shows the fragmentation barrier. The green line shows the drift barrier. The white line shows where St=1. 33 4. Results & Discussion dust through drift from the outer regions and the decrease in density at late times is likely due to the fact that the dust in most of the disk has almost been depleted at this point. That the low fragmentation velocity of 1 ms´1 is the limiting factor for the grain growth and the grain sizes can also be seen from the blue line in Figure 4.2, which represents the analytical solution to the fragmentation barrier, whereas the green line, representing the drift barrier, does not seem to affect the grain sizes except for in the outer disk at late times. A higher fragmentation velocity of 10 ms´1 will allow for the dust particles to grow into larger grain sizes before fragmenting and hence make up a larger spectrum of grain sizes than particles limited by a lower fragmentation velocity, as 1 ms´1. In the case of a higher fragmentation velocity the larger dust grains will be affected by a higher gas drag causing their radial drift to be more effective compared to smaller grains. Hence these particles initially can grow quite large quickly, after which they will have a quick radial drift inwards due to to the higher gas drag. This can be seen by studying Figure 4.3 in which pdfs for the 10 ms´1 fragmentation velocity case is shown. Considering this figure one can see that, compared to the 1 ms´1 case, the disk hasn’t become quite as stable in the inner region of 1 au where the representative grain size drifts with time from 11.06 cm at 0.1 Myr to 0.53 cm at 5 Myr. This is explained by the fact that the higher fragmentation velocity allows for larger grain sizes to be created making the spectrum of grain sizes larger where the largest particles will be the most affected by the gas drag, having the more efficient radial drift. This will depleted the larger grain sizes quicker then the smaller grain sizes in a given part of the disk, causing the position of the peak at even the 1 au radius to drift towards smaller grain sizes over time. By considering the outer region of the disk with fragmentation velocity 10 ms´1 one can notice that this part also is less stable over time compared to the outer region of the 1 ms´1 case. At 10 au the dominant grain size at 0.1 Myr is 2.69 cm compared to 0.011 cm at time 5 Myr. These smaller grain sizes in the outer region compared to the 1 au-region can, in the same way as for the 1 ms´1 fragmentation case, likely be explained by the lower dust density in this region together with the drift of larger grains from the outer to the inner regions. That this is the case can be seen by studying Figure 4.4. In Figure 4.4 it can be seen how the dust density in most of the disk decreases rapidly compared to the 1 ms´1 fragmentation case shown in Figure 4.2. The reason for the quicker density decrease and depletion of dust in the 10 ms´1 fragmentation case is, as discussed, the more efficient gas drag and radial drift. This can also be seen from the fact that the green line in Figure 4.4, indicating the analytical solution for the drift barrier, is the more limiting barrier for the grain sizes compared to the fragmentation barrier which was more dominant in the 1 ms´1 fragmentation case. 34 4. Results & Discussion Figure 4.3: Probability density functions based on the simulation with the frag- mentation velocity of 10 ms´1. The pdfs are shown at two different radii and at each radii five times are shown between 0.1 Myr and 5 Myr. The peaks represent the most probable grain size at that given time and radius. 35 4. Results & Discussion Figure 4.4: 2D plots showing the density evolution for given grain sizes and given radii at four times between 0.1 Myr and 5 Myr for the fragmentation velocity 10 ms´1. The blue line shows the fragmentation barrier. The green line shows the drift barrier. The white line shows where St=1. To more easily see that the dust density is decreasing earlier for higher fragmen- tation velocities Figure 4.5 can be studied. From this figure it is seen that higher fragmentation velocities and hence larger grains indeed has a quicker density de- crease due to the higher gas drag. Before the density decrease one can observe an increase in the dust density which can be interpreted as a wave front of pebbles drifting through that region. One can note that the density in this wave front is higher for higher fragmentation velocities, explained by the larger amount of big particles which are affected by the gas drag and hence drifting inwards through the disk. The density decrease and the wave front of pebbles at two different radii is shown in Figure 4.5. Worth noticing is the scale on the y-axis where the density is significantly lower at the 10 au-region compared to the 1 au-region. 36 4. Results & Discussion Figure 4.5: The dust density is shown over time at two different radii, 1 au (to the left) and 10 au (to the right) for the three different fragmentation velocities, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green). Note the different scales on the y-axis between the two radii. Figure 4.6: The pdfs for grain sizes shown for the three different fragmentation velocities, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green), at 1 au for two dif- ferent times, 0.1 Myr (to the left) and 5 Myr (to the right). 37 4. Results & Discussion Now, it has been seen how the representative grain sizes evolve at different radii for the three fragmentation velocity cases separately and comparison between them have been discussed. To more easily see the fact that a higher fragmentation velocity allows for larger particle sizes and that the pdf peaks drift towards the same grain size over time in the different cases, given a radius, Figure 4.6 can be studied where the pdf evolution for the three fragmentation velocities at 1 au is shown at an early (0.1 Myr) and at a late (5 Myr) time. To study the difference between the three fragmentation velocities at an outer radius as 10 au Figure 4.7 can be studied. Here a similar pattern as in Figure 4.6 can be noticed, but the grain sizes are in general smaller in this outer region due to a lower dust density as already discussed. Figure 4.7: The pdfs for grain sizes shown for the three different fragmentation velocities, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green), at 10 au for two different times, 0.1 Myr (to the left) and 5 Myr (to the right). To get an overview of the representative grain sizes at different fragmentation velocities, radii and time, the grain sizes discussed are summarized in Table 4.1. Table 4.1: The representative grain sizes for the fragmentation velocities 1 ms´1 and 10 ms´1 at 1 au and 10 au at times 0.1 Myr, 1 Myr and 5 Myr summarized. Radius Time vfrag “ 1 ms´1 vfrag “ 10 ms´1 1 au 0.1 Myr 0.59 cm 11.06 cm 1.0 Myr 0.59 cm 3.00 cm 5.0 Myr 0.59 cm 0.53 cm 10 au 0.1 Myr 0.27 cm 2.69 cm 1.0 Myr 0.25 cm 0.14 cm 5.0 Myr 0.012 cm 0.011 cm 38 4. Results & Discussion 4.2 Dust Fluxes and Mass Delivery The total dust flux through the 1 au radius to the inner region of the disk as well as the total dust mass delivered to that same region is presented in Figure 4.8. The flux and the total mass delivered was investigated for three different fragmentation velocities, 1 ms´1, 5 ms´1 and 10 ms´1, all of which is presented in the same figure. Figure 4.8: Graphs over the mass flux through the 1 au radius over time (to the left) and graphs over the total mass delivered to the region within 1 au (to the right). In both cases three different fragmentation velocities are shown, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green). As can be seen from Figure 4.8 the dust accretion rate is quite constant and similar for all the fragmentation velocities up until about 2 ˆ 102 yrs. Here one can see an increase in the dust flux for all of the fragmentation velocities indicating a wavefront of dust and pebbles through the 1 au radius. Here one should notice that the accretion rate is higher for higher fragmentation velocities which also mean a quicker dust delivery. This can be seen through the dust delivery plot in Figure 4.8 where the dust delivered to the region within 1 au at 0.1 Myr is about 4 M‘ for the fragmentation velocity 1 ms´1, 63 M‘ for the fragmentation velocity 5 ms´1 and 78 M‘ for 10 ms´1. That higher fragmentation velocities result in a higher flux and a quicker mass delivery is likely explained through the fact that a higher fragmentation velocity allows for the dust particles to grow into larger sizes without fragmenting. The larger grain sizes would be affected by a more efficient gas drag causing a quicker radial drift. By studying the dust mass delivered to the region within 1 au one can see that the mass delivered starts to stabilize between all the fragmentation velocities at about 1.5 Myr. The dust mass delivered at this point for the different cases is 95 M‘ for 1 ms´1 and 99 M‘ for both the 5 ms´1 and 10 ms´1 fragmentation case. 39 4. Results & Discussion Then, after 5 Myr, at the end of the simulations, all the fragmentation velocities has stabilized even more and the total mass delivered in all the fragmentation cases was about 99 M‘. At a given time and radius it was investigated which grain sizes that were the most responsible for the mass delivery in the different fragmentation cases. When the mass delivery through the 1 au-radius was investigated the mass weighted pdfs in Figure 4.9 were obtained. In Figure 4.9 the peaks show which grain size is the most dominant for the mass delivered through the 1 au radius after 0.1 Myr, 1 Myr and 5 Myr for the three different fragmentation velocities. By studying Figure 4.9 it was found that the grain sizes responsible for the most delivered dust mass after 1 Myr was 0.65 cm for the 1 ms´1 fragmentation case, 7.16 cm for the 5 ms´1 fragmentation case and 12.33 cm for the 10 ms´1 fragmentation case. The higher fragmentation velocities allow for particles to stick and grow into larger grain sizes. Since the mass goes as mpaq9a3 it’s reasonable that most of the mass is delivered at early times by large grains. This also explain why the curves in Figure 4.9 seem to have stabilized after 0.1 Myr as well. Figure 4.9: Mass weighted pdfs over the dominant grain sizes for the mass delivered to the region within the 1 au radius at three different times for the three different fragmentation velocities. When it for the application of IOPF was investigated how much mass of grain sizes larger than 0.06 cm that was delivered to the 1 au radius it was found that 77 M‘ had been delivered after 1 Myr in the 1 ms´1 fragmentation case, out of which 94.4% where made up of grains larger than 0.06 cm. In the 5 ms´1 case 99 M‘ had been delivered after 1 Myr and 99.8% were made up of grains larger than 0.06 cm. The same total mass had been delivered in the 10 ms´1 fragmentation case, but in this case 99.9% of the mass was delivered by grains larger than 0.06 cm. 40 4. Results & Discussion When investigating how much time it would require to deliver about 3 M‘ to the region within 1 au Figure 4.8 can be studied. An enlargement of the mass delivery plot from Figure 4.8 with the 3 M‘ level marked can be seen in Figure 4.10. The time for delivering 3 M‘ to the region within 1 au can be found to be 9ˆ104 yrs in the 1 ms´1 fragmentation case, 8ˆ103 yrs in the 5 ms´1 fragmentation case and 5ˆ103 yrs in the 10 ms´1 fragmentation case. The quicker delivery for higher fragmentation velocities is explained by the larger grain sizes allowed and hence the more efficient gas drag and radial drift. Figure 4.10: Total mass delivered through the 1 au radius over time. Three frag- mentation velocities are shown, 1 ms´1 (blue), 5 ms´1 (orange) and 10 ms´1 (green). The 3 M‘ level is also shown (black dashed line). 4.3 Affect of Radial Resolution When studying how the resolution of the radial grid would affect the results two cases were compared, simulations ran with 64 radial grid cells and simulations ran with 128 radial grid cells. In both cases the fragmentation velocity 1 ms´1 was used for the sake of comparison. In order to study the resolution’s affect on the results, pdfs of the grain sizes were created for two different regions, 1 au and 10 au, at both of which pdfs for the grain size at 0.1 Myr and 5 Myr were plotted. This was done with data generated by simulations using the two different radial resolutions and the results can be seen in Figure 4.11. 41 4. Results & Discussion Figure 4.11: Probability density functions for grain sizes shown at two different radii, 1 au (to the left) and 10 au (to the right). At each radius two times are shown 0.1 Myr (solid line) and 5 Myr (dashed line). For each radius and time two graphs are shown, one based on simulations with 64 radial grid cells (blue graphs) and one based on simulations with 128 radial grid cells (red graphs). As can be seen from Figure 4.11 the two investigated radial resolutions does not seem to have a major impact on the results, regarding the grain size distributions, since the graphs for the same radius and time overlap each other quite well in the two resolution cases. At the inner region, 1 au, at early times there are a slight deviation between the two resolutions where the 64 grid resolution gives a grain size peak at 0.65 cm compared to 0.59 cm for the 128 grid resolution at 0.1 Myr. This is a deviation of about 9.2%. Studying the same 1 au-region at a later time of 5 Myr the two resolution cases have converged and both cases give a peak position of 0.59 cm. Considering the outer region of 10 au the grain size peaks for the two different radial resolutions match at early times of 0.1 Myr where both peaks indicate a representative grain size of 0.27 cm. The peak positions for the two resolution cases continue to match each other as the time evolves to 5 Myr, at which point both peaks indicate a dominant grain size of 0.012 cm. Worth noticing is also that, even though the peak positions in many of the cases discussed are the same, the peak height differs slightly between the 64 grid and 128 grid cases. The peak heights for the different resolution cases are very similar at early times as 0.1 Myr which can be seen from Figure 4.11. However, going to later times as 5 Myr the peak for the 128 grid resolution is higher then the peak for the 64 grid resolution, indicating a higher probability. This is likely explained by the fact that a higher radial resolution allows for the distributions to be calculated with a higher accuracy at a given radius due to the decreased grid spacing compared to the case with 64 radial grid cells. 42 5 Conclusions In this chapter the conclusions drawn from the results are being presented. This is being done in three sections where the grain size distributions for the different frag- mentation velocities first are considered. Then the conclusions drawn from studying the dust mass flux and mass delivery is being presented. Finally the conclusions drawn from comparing different radial resolutions are stated. 5.1 Dominant Grain Sizes By studying the pdfs and the grain sizes presented in Section 4.1 it can be concluded that the largest grains are obtained in the inner disk at early times. Larger grains are also obtained for higher fragmentation velocities. The dominant grain size at 1 au in an early disk is 0.59 cm at 0.1 Myr for a fragmentation velocity of 1 ms´1 whereas this grain size increase to 11.06 cm for the same radius and time if the fragmentation velocity increase to 10 ms´1. Larger particles are subject to a more efficient radial drift and with time the dust density decreases. Hence the grain size distributions can stabilize at a certain grain size with time, despite different fragmentation velocities. Due to the depletion of dust this stable grain size is smaller than the representative grain sizes at earlier times. At time 5 Myr and at 1 au the representative grain size given by a fragmentation velocity of 1 ms´1 is 0.59 cm whereas the corresponding grain size given by the fragmentation velocity 10 ms´1 is 0.53 cm. Hence the spacing between the pdf peaks for different fragmentation velocities decreases with time towards one mutual grain size. The same pattern occur in the outer part of the disk where the dominant grain sizes at 10 au are 0.27 cm for 1 ms´1 and 2.69 cm for 10 ms´1 at time 0.1 Myr. These grain sizes have at time 5 Myr drifted to 0.012 cm for 1 ms´1 and 0.011 cm for 10 ms´1. The smaller grain sizes in this outer part of the disk, compared to the inner disk, is explained by the lower dust density. 5.2 Dust Mass Flux and Mass Delivery The flux of dust through a given radius, which in this thesis was chosen to be 1 au, have been observed to be higher given a higher fragmentation velocity. This is due to the more efficient gas drag exerted on the larger particles, allowed by the higher fragmentation velocity, resulting in a quicker radial drift. The mass delivery to the inner region hence becomes quicker when higher fragmentation velocities are allowed. For 1 ms´1 the delivered mass to the inner region was 4 M‘ at time 0.1 Myr compared to 78 M‘ for the fragmentation velocity 10 ms´1 at that same 43 5. Conclusions time. However, after a certain time, the delivered mass seem to converge towards the same value, regardless of the fragmentation velocity. This mass is about 99 M‘ and is reached after about 1.5 Myr. In the different fragmentation cases the grain sizes responsible for most of the delivered mass are 0.65 cm, 7.16 cm and 12.33 cm for the 1 ms´1, 5 ms´1 and 10 ms´1 fragmentation case respectively. After 1 Myr the total mass delivered through the 1 au radius were 77 M‘ for 1 ms´1 and 99 M‘ for 5 ms´1 and 10 ms´1. The percentage of the mass delivered by grains larger than 0.06 cm in the three different cases were 94.4%, 99.8% and 99.9% respectively. It is also concluded that the time required for delivering 3 M‘ to the region within 1 au for the different fragmentation cases are 9 ˆ 104 yrs, 8 ˆ 103 yrs and 5 ˆ 103 yrs respectively. 5.3 Radial Resolution Considering the two different radial resolutions investigated in this thesis, 64 radial grid cells and 128 radial grid cells, it has been found that the difference in pdf peak positions are close to non or zero. This can be seen from Figure 4.11, where the graphs overlap each other very well. One can see that a higher radial resolution might increase the accuracy of the grain size distributions shown by the slightly higher peaks in the case where 128 radial grid cells were used. Since the peak positions doesn’t change significantly or at all depending on the resolutions investigated and given the fact that the graphs overlap each other well, it can be concluded that the simulations done have been converged regarding radial resolution. 44 Bibliography [1] P. J. Armitage and W. Kley, From Protoplanetary Disks to Planet Formation. Berlin: Springer, 1 ed., 2019. URL: https://link.springer.com/book/10. 1007/978-3-662-58687-7. [2] S. Chatterjee and J. C. Tan, “Inside-out planet formation,” The Astrophysical Journal, vol. 780, no. 1, 2014. DOI: 10.1088/0004-637X/780/1/53. [3] M. Mayor and D. Queloz, “A jupiter-mass companion to a solar-type star,” Nature, vol. 378, pp. 355–359, 1995. DOI: 10.1038/378355a0. [4] D.-C. Chen, J.-W. Xie, J.-L. Zhou, S. Dong, J.-Y. Yang, W. Zhu, C. Liu, Y. Huang, M.-S. Xiang, H.-F. Wang, Z. Zheng, A.-L. Luo, J.-H. Zhang, and Z. Zhu, “The evolution of hot jupiters revealed by the age distribution of their host stars,” PNAS, vol. 120, no. 45, 2023. DOI: 10.1073/pnas.2304179120. [5] National Aeronautics and Space Administration, ’Kepler: NASA’s First M