Understanding The Accretion Process of High Mass Protostars MARTIN ALMQVIST Space, Earth and Environment CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2025 www.chalmers.se ii www.chalmers.se Master’s thesis 2025 Understanding The Accretion Process of High Mass Protostars MARTIN ALMQVIST Space, Earth and Environment Astronomy and Plasma Physics Chalmers University of Technology Gothenburg, Sweden 2025 Understanding The Accretion Process of High Mass Protostars MARTIN ALMQVIST © MARTIN ALMQVIST, 2025. Supervisor: Brandt Gaches, Space, Earth and Environment Examiner: Jonathan Tan, Space, Earth, and Environment Master’s Thesis 2025 Space, Earth and Environment Astronomy and Plasma Physics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: The simulated accretion disk around a protostar is depicted as a density heatmap. Regions with lower density appear more diffuse, while the densest areas are shown in red-yellow. Typeset in LATEX Gothenburg, Sweden 2025 v Understanding The Accretion Process of High Mass Protostars MARTIN ALMQVIST Department of Space, Earth, and Environment Division of Astronomy and Plasma Physics Chalmers University of Technology Abstract This thesis presents a high-resolution magnetohydrodynamics simulation of a high- mass protostar with accretion rates above 10−4M⊙ yr−1 with radius R∗ = 33R⊙. The thesis will present the results of different models of accretion rate, accretion flow, and the locations on the disk where accretion occurs. The models will feature five protostellar magnetic fields: 50 G, 500 G, 1 kG, 2 kG, and 4 kG. The initial surface density will also be varied for four models: 200 g cm−2, 500 g cm−2, 1 kg cm−2, and 2 kg cm−2. The main finding is that a protostellar magnetic field leads to more substantial but uneven accretion rates. For the 4 kilogauss model, the disk was truncated 88% of the simulation, and 69% of the total cumulative mass came from accretion rates six times the mean accretion rate, which is a sign of magnetospheric accretion. In comparison, the 50 Gauss model had 6% of the total cumulative mass in rates above 6 the mean accretion rate, and the accretion came from a turbulent boundary layer. The accretion mainly occurred in a wide range, but as the protostellar magnetic fields increased in strength, the accretion shifted from an equatorial accretion to happen more at higher latitudes. The lower initial surface density also led to a higher truncation time with the in- creased truncation radius. The two models 200 g cm−2, 500 g cm−2 were seen to have magnetospheric accretion, while the 1 k g cm−2 and the 2 k g cm−2 had a turbu- lent boundary layer with material hitting the protostar at a wide range. The lower models also had a more burst-driven accretion caused by gas buildup at the trun- cation radius, which is disrupted when ram and thermal pressures exceed magnetic pressure. Keywords: High-mass protostars, Magnetohydrodynamics (MHD), Accretion disks, Protostellar magnetic fields, Disk surface density vi Acknowledgements First, I want to express my deepest gratitude to my supervisor, Brandt Gaches. This thesis would not have been possible without your guidance and support. Your willingness to answer my questions at all hours, provide thoughtful feedback, and regularly check in on my progress has meant a lot. I genuinely appreciate your mentorship and everything you’ve taught me. I also want to thank my examiner, Jonathan Tan, for allowing me to participate in this project and overseeing its progress. Your insights and expertise have been incredibly valuable. To my family—thank you for always supporting me throughout my years at Chalmers. It’s been a journey filled with challenges, learning, and fun. I couldn’t have done it without you. Martin Almqvist, Gothenburg, March 2025 List of Acronyms Below is the list of acronyms that have been used throughout this thesis listed in alphabetical order: AU Astronomical Unit CDF Cumulative Distribution Function HPC High-Performance Computing MHD Magnetohydrodynamics PDF Probability Density Function PLUTO A Computational MHD Code x Nomenclature Below is the nomenclature parameters, and the variables that have been used through- out this thesis. Parameters AU Astronomical Unit B∗ Initial Magnetic field strength c Speed of light cs Sound speed G Gravitational constant H Scale height of the disk L Luminosity M⊙ Solar Mass Ptherm Thermal pressure Pkin Kinetic pressure Pmag Magnetic pressure Prad Radiation pressure R⊙ Solar Radius Σc Initial Surface density of the disk T Temperature τ Optical depth vff Free-fall velocity vϕ Azimuthal velocity vr Radial velocity vθ Vertical velocity β Plasma beta (ratio of gas to magnetic pressure) ṁ Mass accretion rate xi ρ Density Ω Angular velocity ωs Fastness parameter ν Kinematic viscosity ξ Ionization fraction xii Contents List of Acronyms ix Nomenclature x List of Figures xv List of Tables xix 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Aim and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 3 2.1 Star formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Accretion Disk . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.2 Alpha Modeling of Disk Viscosity . . . . . . . . . . . . . . . . 6 2.1.3 Turbulent and Rotational Dynamics of Accretion Disks . . . . 7 2.1.4 Rotational and Keplerian Timescales . . . . . . . . . . . . . . 7 2.1.5 Magnetic Field Thresholds and Disk Truncation . . . . . . . . 8 3 Methods 11 3.1 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1.1 Changes from HADES . . . . . . . . . . . . . . . . . . . . . . 12 3.1.2 Heating and cooling . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.3 Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.4 Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.5 Parameter Variations . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.5.1 Stellar Magnetic Field Strength . . . . . . . . . . . . 16 3.1.5.2 Surface Density . . . . . . . . . . . . . . . . . . . . . 16 3.1.6 Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Computational Methods . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3.1 Accretion filling fraction . . . . . . . . . . . . . . . . . . . . . 17 3.3.2 Truncation radius . . . . . . . . . . . . . . . . . . . . . . . . . 18 xiii Contents 4 Results 19 4.1 Accretion Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.1.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 20 4.1.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 21 4.2 Mass Accretion Bursts . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.2.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 22 4.2.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 24 4.3 Location of accretion on the protostellar surface . . . . . . . . . . . . 26 4.3.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 26 4.3.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 30 4.4 Truncation Radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.4.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 32 4.4.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 33 4.5 Disk Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.5.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 34 4.5.1.1 Radial and Azimuthal Flow Variations . . . . . . . . 35 4.5.1.2 Thermal, Magnetic, and Kinetic Pressure . . . . . . 39 4.5.1.3 Magnetic and Plasma . . . . . . . . . . . . . . . . . 39 4.5.1.4 Density profile . . . . . . . . . . . . . . . . . . . . . 41 4.5.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 42 4.5.2.1 Radial and Azimuthal Flow Variations . . . . . . . . 43 4.5.2.2 Thermal, Magnetic, and Kinetic Pressure . . . . . . 48 4.5.2.3 Magnetic and Plasma . . . . . . . . . . . . . . . . . 49 4.5.2.4 Density Profile . . . . . . . . . . . . . . . . . . . . . 50 4.6 Gas-Magnetic Pressure Interactions . . . . . . . . . . . . . . . . . . . 52 4.6.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 52 4.6.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 53 4.7 Fastness Parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.7.1 Magnetic Field Variation (B∗) Series . . . . . . . . . . . . . . 55 4.7.2 Surface Density Variation (Σo) Series . . . . . . . . . . . . . . 56 4.8 Resolution Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5 Conclusion 59 xiv List of Figures 3.1 Example of a filling fraction calculation for the snapshot at day 259 in the B∗ = 4000 G, Σo = 2000 g cm−2 model. The red line represents the sorted mass flux at the protostar’s surface, the green line rep- resents the accretion rate, the black line shows the cumulative mass accretion, and the blue line depicts the cumulative area as a function of the sorted cell index. The vertical dashed lines indicate the indices at which the top 10%, 25%,50%, 75%, and 90% of the total mass is accreted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.1 The figures show the evolution from an initial state of a protostar for four different initial magnetic field strengths of the protostar and a comparison to the 2.4 × 10−4 line which represents the predicted mass accretion for an 8M⊙ protostar [1]. The right image compares changing the initial surface density of the disk, Σo. . . . . . . . . . . 20 4.2 The figure shows the mean and standard deviation of the mass ac- cretion for the two simulation series B∗ Σo. A transient time of 50 days was used to erase the impact of the starting condition. The dots mark the mean accretion for each model, and the lines are the standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.3 Probability density function (PDF) and cumulative distribution func- tion (CDF) of mass accretion rates normalized by the mean accretion rate, ṁ/⟨ṁ⟩mean, for a transient time of 50 days. The left panel (PDF) shows the distribution of accreted mass as a function of the normalized accretion rate for different initial magnetic field strengths of the protostar (B∗). The right panel (CDF) presents the cumulative mass fraction accreted as a function of the normalized accretion rate. The thick lines in the CDF represent the contribution from accretion rates above the mean rate, emphasizing the fraction of mass accreted at higher rates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.4 A box plot for the different initial protostar magnetic field strengths, B∗. The figure shows how much of the accreted mass occurs below and above the mean, and then after 3 × ⟨ṁ⟩ and 6 × ⟨ṁ⟩ for each model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 xv List of Figures 4.5 Probability density function (PDF) and cumulative distribution func- tion (CDF) of mass accretion rates normalized by the mean accretion rate, ṁ/⟨ṁ⟩mean, for a transient time of 50 days. The left panel (PDF) shows the distribution of accreted mass as a function of the normal- ized accretion rate for different surface density values (Σ0). The right panel (CDF) presents the cumulative mass fraction accreted as a func- tion of the normalized accretion rate. . . . . . . . . . . . . . . . . . . 25 4.6 A box plot representing the various initial surface densities, Σo. The figure illustrates the amount of accreted mass that occurs below the mean, and subsequently after 3 × ⟨ṁ⟩ and 6 × ⟨ṁ⟩ for each model. . 26 4.7 Accreting rate versus time and latitude on the stellar surface for four simulations with different initial disk densities, B∗, as annotated in each sub-figure. The color coding represents the magnitude of the accreting momentum flux, with darker shades indicating lower flux and lighter shades indicating higher flux. . . . . . . . . . . . . . . . 27 4.8 Accretion rate and surface area coverage for the B∗-series. The top panel shows the accretion rate for different magnetic field strengths over time. The middle and bottom panels display the fraction of the protostar surface responsible for 90% and 50% of the total accreted mass, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.9 Accretion rate versus time and latitude on the stellar surface for four simulations with different initial disk densities, Σo, as annotated in each sub-figure. The color coding represents the magnitude of the accreting momentum flux, with darker shades indicating lower flux and lighter shades indicating higher flux. . . . . . . . . . . . . . . . . 30 4.10 Accretion rate and surface area coverage for the Σo series. The top panel shows the accretion rate over time for different magnetic field strengths. The middle and bottom panels display the fraction of the protostar surface responsible for 90% and 50% of the total accreted mass, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.11 Figure visualizes the truncation radius as a function of time. The disk is truncated where the black dots are for the B∗ series. The red dotted line represents the mean truncation radius for each model. . . 32 4.12 Figure visualizes the truncation radius as a function of time. The disk is truncated where the black dots are for the Σo series. . . . . . . 33 4.13 Temperature, density, entropy, and pressure profiles for the B∗ series, averaged over a midplane wedge and binned radially into 32 bins. The solid lines represent the profile at t = 259 days of simulation, while the dotted lines correspond to t = 16 days. . . . . . . . . . . . . . . . 34 4.14 The radial, poloidal, and azimuthal velocity components for two snap- shots, day 16 and day 259, for the B∗ series. The Keplerian velocity is subtracted from vϕ. . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.15 The figure shows the variation of poloidal velocity as a function of radius in the disk for different magnetic field strengths B∗. The left figure is the mean polodial velocity, while the right one shows the standard deviation of the same velocity. t = 259 days. . . . . . . . . . 36 xvi List of Figures 4.16 The figure shows the total velocity dispersion and mean velocity pro- files at different times. The left panels represent the total velocity dispersion σ(Vtotal) (solid lines) and the mean velocity ⟨Vtotal⟩ (dot- ted lines). The right sub figures shows the dispersion of the velocity component, vphi − vkep normalized by the sound speed. . . . . . . . . 37 4.17 The figure shows the radial velocity normalized by the free-fall veloc- ity, vr/vff , as a function of radius for different magnetic field strengths B∗ and time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.18 The figure shows the variation of the thermal, magnetic, and kinetic pressure components with radius for two snapshots: day 16 and the final snapshot at day 259 for the B∗ series. . . . . . . . . . . . . . . . 39 4.19 Evolution of the total magnetic field strength and plasma beta in the disk for the B∗ series. The values are computed as the mean over a wedge around the midplane. The left panels show the total magnetic field strength, while the right panels display plasma beta (β) as a function of radius at different snapshots. The black dotted line is where β = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.20 This figure presents a color-coded density map for various B∗ models, taken at a snapshot on day 259. The black streamlines indicate the direction of poloidal velocity, vpol. . . . . . . . . . . . . . . . . . . . . 42 4.21 Temperature, total pressure, density, and entropy profiles for the Σo series, averaged over a midplane wedge and binned radially into 32 bins. The solid lines represent the profiles at t = 259 days, while the dashed lines correspond to t = 16 days. Different colors indicate models with varying initial surface density Σo = 200, 500, 1000, and 2000 g cm−2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.22 The figure illustrates the radial velocity profile for two snapshots, day 16 and the last snapshot day259 for the Σo series. . . . . . . . . . . 44 4.23 Variation of poloidal velocity as a function of radius in the disk for different initial surface densities Σ∗. The left panel shows the mean poloidal velocity, while the right panel displays the standard deviation of the same velocity. The final snapshot of the simulation at t = 259 days was used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.24 Velocity dispersion and mean velocity for the Σo series at different snapshots, showing the evolution of turbulent motion in the disk. . . 46 4.25 The figure shows the radial velocity normalized by the free-fall veloc- ity, vr/vff , as a function of radius for different surface densities Σo and time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.26 Variation of the thermal, magnetic, and kinetic pressure components with radius for two snapshots: day 16 and the final snapshot at day 259 for the Σo series. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.27 Evolution of the total magnetic field strength and plasma beta in the disk for the Σo series. The values are computed as the mean over a wedge around the midplane. The left panels show the total magnetic field strength, while the right panels display plasma beta (β) as a function of radius at different snapshots. . . . . . . . . . . . . . . . . 50 xvii List of Figures 4.28 This figure presents a color-coded density map for various Σo models, taken at a snapshot on day 259. The black streamlines indicate the direction of poloidal velocity, vpol. . . . . . . . . . . . . . . . . . . . . 51 4.29 The figure shows the ratio between the poloidal ram pressure (ρv2 p) plus the thermal pressure over the poloidal magnetic pressure for the disk region (R < 3 AU) for the five different B∗ models. The lines and arrows indicate the direction of the flow. . . . . . . . . . . . . . . 53 4.30 The figure shows the ratio between the poloidal ram pressure (ρv2 p) plus the thermal pressure and the poloidal magnetic pressure for the disk region (R < 3 AU) across four different Σo models. The lines and arrows indicate the direction of the flow. . . . . . . . . . . . . . . 54 4.31 The figures show the fastness parameter, ωs, as a function of time for the five different protostar models of varying magnetic fields. The dotted line is where ωs = 1. . . . . . . . . . . . . . . . . . . . . . . . 55 4.32 The figures illustrate the fastness parameter, ωs, as a function of time for the four different protostar models with varying disk surface density. The dotted line indicates where ωs = 1. The sharp lines represent a numerical error, occurring when the truncation radius is 0, causing ωs to explode. . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.33 The figure shows the mass accretion rate over time for different res- olutions: 128, 256, and 512. The red line represents the theoretical accretion rate of 2.67 × 10−4M⊙/year. . . . . . . . . . . . . . . . . . . 57 xviii List of Tables 4.1 Summary of simulation models with varying protostellar magnetic field strength (B∗) and disk surface density (Σo). . . . . . . . . . . . . 19 4.2 Mass accretion rates for different magnetic field strengths B∗ with fixed Σo = 2000 g cm−2. . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.3 Mass accretion rates for B∗ = 2000 G with varying disk surface den- sities Σo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.4 Fraction of total accreted mass occurring above specific accretion rate thresholds for different B∗ values. . . . . . . . . . . . . . . . . . . . . 23 4.5 Fraction of total accreted mass occurring above specific accretion rate thresholds for different Σo values. . . . . . . . . . . . . . . . . . . . . 25 4.6 Filling fraction (percent) of the accretion flow onto the protostar ac- cording to the top X% of the mass accretion for B∗ (G). . . . . . . . 29 4.7 Filling fraction (percent) of the accretion flow onto the protostar ac- cording to the top X% of the mass accretion for Σo (g cm−2). . . . . . 31 4.8 Truncation Radius and Time Fraction for B∗ Model . . . . . . . . . . 32 4.9 Truncation Radius and Time Fraction for Σo Model . . . . . . . . . 33 xix List of Tables xx 1 Introduction 1.1 Background For as long as humans have existed, we have looked up at the night sky, wondering what those tiny dots of light are. At first, we thought they were just distant lights, part of the heavens, or even divine beings. But over time, as our knowledge grew, we began to understand that these dots were not just lights, they were giant burning spheres of gas, much like our own Sun. With the invention of telescopes, our understanding grew. We realized that space is vast and that our Sun is just one of billions of stars in a galaxy known as the Milky Way. Ultimately, we learned that our galaxy is just one among trillions in the universe. As we studied further, we discovered how stars evolve, how they function, why they shine, and how they ultimately perish—sometimes as spectacular supernovae and other times as quietly fading dwarfs. However, one of the greatest mysteries remains: how are stars born? While we have a good understanding of stellar evolution and death, the star formation process is much more challenging to observe. Stars take millions of years to form, far longer than a human lifetime. We can’t just watch one being born in real-time. Young stars are often hidden behind thick dust clouds, making them difficult to see even with the most powerful telescopes. Because of these challenges, astronomers don’t rely solely on direct observation. Instead, advanced simulations play a crucial role in unraveling the mysteries of star formation. These simulations allow us to recreate the conditions of space, test different theories, and evolve virtual stars over millions of years. By comparing these models with accurate astronomical data, we refine our understanding, confirm or adjust our theories, and push the boundaries of what we know. Each breakthrough in simulation technology brings us closer to uncovering the actual process of how stars are born. By studying the universe, we do more than uncover its secrets, we gain a deeper understanding of our existence. What once seemed like tiny, scattered lights in the sky are now known to be distant suns, each part of a vast and ever-changing cosmos. Yet, for all we have learned, the universe still holds countless mysteries, waiting to be revealed. With every discovery, we take one step closer to understanding not just the stars, but our place among them. 1 1. Introduction 1.2 Aim and objectives This study investigates how the magnetic field strength of a protostar and the sur- face density of its surrounding disk influence the mass accretion rate and the spatial distribution of accretion. Using high-resolution ideal magnetohydrodynamic (MHD) simulations with the PLUTO and HADES codes, this research examines how these properties affect the dynamics of accretion flows, disk truncation, velocity distribu- tions, and pressure profiles in the early stages of protostellar evolution. 2 2 Theory 2.1 Star formation Star formation is an essential process in the universe, it drives the evolution of galaxies and the creation of life itself. Without this process, we would not have had any metals crucial for sustaining life, especially since the high mass stars are the most driven force for the process. It is therefore of high interest to understand how high mass stars are born and the processes that govern the formation, such as the accretion of nearby materials, and the feedback mechanism in the form of jets created as the protostar accretes. The process begins in molecular clouds, often referred to as stellar nurseries, which contain vast amounts of gas and dust. In particular, star formation typically oc- curs within Giant Molecular Clouds (GMCs), massive structures spanning tens to hundreds parsecs, and are typically composed of cold molecular hydrogen. The typ- ical temperature of a GMC is around T ∼ 10 K, and the typical density ranges from 102 cm−3 to 104 cm−3. The mass of these clouds can be on the order of 103 −106 M⊙. Within GMCs, star formation begins when regions of gas collapse under their own gravity, forming dense clumps and filaments. This gravitational collapse is initially governed by the Jeans instability, which sets the conditions under which gravity can overcome thermal pressure, leading to runaway collapse. The Jeans mass defines the minimum mass necessary for a region of gas to become gravitationally unstable The Jeans mass sets the minimum mass required for a region of gas to become gravitationally unstable, MJ = ( 5kBT GµmH )3/2 ρ−1/2 (2.1) where T is the temperature, ρ is the density, µ is the mean molecular weight, mµ is the mass of a hydrogen atom and G is the gravitational constant A gas region that exceeds the Jeans mass becomes gravitationally unstable and collapses. However, thermal pressure is not the only opposing force, turbulence and magnetic fields significantly influence the collapse process, which can delay, prevent, or redirect the collapse. As the gas contracts, it develops dense cores within the GMC. The infalling gas cannot directly fall into the protostar if a collapsing core has enough angular momentum. Instead, it creates a rotational supported accretion disk, which serves as a reservoir, supplying the protostar over time. 3 2. Theory 2.1.1 Accretion Disk The gas surrounding the protostar possesses high angular momentum, which hin- ders direct collapse and facilitates the formation of an accretion disk. If angular momentum were perfectly conserved, the gas would settle into a circular orbit at a characteristic centrifugal radius Rc, Rc = j2 GM∗ (2.2) where j = vr is the specific angular momentum of the gas. However, for material to move inward and accrete onto the protostar, it must lose angular momentum. This loss occurs through viscous torques, which transfer angular momentum outward while allowing the gas to move inward. The balance between gravitational pull, centrifugal forces, and viscosity governs the evolution of the disk. To create a simple model, we can use the thin disk approximation. This allows us to assume that the disk is thin, meaning the vertical height, H, is relatively small compared to the radial component. We also assume that the disk is axisymmetric. By doing so, and assuming that the disk mass is much smaller than the protostar mass, we can neglect the self-gravity within the disk, which allows the angular velocity to align with the Keplerian angular velocity. Ω = √ GM∗ R3 (2.3) The angular velocity is calculated by simply multiplying with the radius, vθ = Ωr. Notably, the rotational velocity depends on the radius for the thin disk. The surface density of a disk, Σ, is determined by the integral, Σ(r) = ∫ z −z ρ(r, z)dz (2.4) Where z is the disk height, which is set to H for a thin disk. Angular momentum transport in an accretion disk occurs due to viscous torques. These torques, denoted as T (r), transfer angular momentum from inner disk regions to outer regions, enabling mass accretion toward the protostar. The transport of Angular momentum outwards in the disk is possible due to the viscosity ν and it is given by: T (r) = 2πνΣr3 dΩ dr . (2.5) The accretion disk can be viewed as radially divided into thin layers, each interacting with its neighboring layer. For gas to move inward, the inner layer must lose angular momentum, which gets transferred outward to the adjacent layer. Simultaneously, the outer layer must transfer mass inward to conserve total mass and maintain a steady-state configuration. The existence of torque arises from the shear flow in a Keplerian accretion disk. Because the disk rotates differentially, the inner layers move at higher angular velocities than the outer ones. This velocity gradient creates a shear between adjacent layers, resulting in the formation of viscous torques. The inner layer exerts a forward drag on the outer layer, while the outer layer pulls back 4 2. Theory on the inner one. This interaction leads to angular momentum being transported outward, while gas slowly spirals inward due to the loss of rotational support. If viscosity were absent, there would be no torque, and therefore no angular momen- tum transport, meaning the gas would remain in stable Keplerian orbits without accreting onto the protostar. Additionally, viscosity converts some of the angular momentum into heat radiated outward from the disk surface. This radiative cooling affects the thermal structure of the disk and influences its emission properties and observed luminosity. The surface density of the disk, is related to mass accretion via the equation: Ṁ = −2πRΣ(r)vr, (2.6) In this context, the radial velocity, denoted as vr, reflects the gas movement towards the protostar, with a negative sign indicating inflow. Similarly, modifying the initial surface density, Σ(R), affects the mass accretion rate, since Σ(R) varies with radius. We can use the result to determine the conservation of mass and angular momentum, R ∂Σ ∂t + ∂ ∂R (RΣvR) = 0. (2.7) This ensures that any change in the surface density Σ is accounted for by a corre- sponding mass transport within the disk. Angular momentum conservation is given by: R ∂ ∂t (ΣR2Ω) + ∂ ∂R (ΣvRR3Ω) = 1 2π ∂T ∂R . (2.8) We can rewrite this into a simple equation using the torque equation, ∂Σ ∂t = 3 R ∂ ∂R [ R1/2 ∂ ∂R ( νΣR1/2 )] . (2.9) This straightforward yet effective equation indicates how the disk will evolve, de- pending on the viscosity, ν. In an accretion disk, the motion of gas is primarily governed by three forces acting in the radial direction: gravitational force, centrifugal force, and pressure gradient force. In a steady-state disk, these forces must balance, ensuring that material follows stable orbits while gradually moving inward due to viscosity. The radial momentum equation for the steady state, assuming Keplerian rotation, is: vr ∂vr ∂r − v2 ϕ r + 1 ρ ∂P ∂r = −GM∗ r2 , (2.10) where ρ is the gas density and (P ) is the gas pressure. The first term represents the change in the radial direction; if vr = 0, this would indicate that the disk is in circular motion. The second term arises as an outward force generated by the centrifugal force due to rotation. The third term is the pressure gradient term of the gas, and the last term accounts for the acceleration caused by the pull of gravity. 5 2. Theory The disk can also be in hydrostatic equilibrium, meaning it is balanced between gravity and pressure forces in the vertical direction of the disk. 1 ρ ∂P ∂z = −GN∗z R3 , (2.11) which for a thin disc has the solution, ρ ≈ ρ0exp ( − z2 2h2 ) (2.12) Where h = cs/ΩK is the scale height and ΩK is the Keplerian rotation rate, ΩK =√ Gm∗/R3. The gas and dust will be influenced by friction, leading to a loss of angular momen- tum in the disk, causing the gas to fall inward. As this occurs, gas and dust release energy due to the loss of gravitational energy. This energy will convert to rotational kinetic energy, and thermal energy will increase and radiate outward from the disk. One can thus study the released energy, or flux, in the form of luminosity generated by the accretion, (Lacc ) to approximate the mass accretion (ṁ), Lacc = GM∗ ṀR∗ . (2.13) 2.1.2 Alpha Modeling of Disk Viscosity In the thin disk model, angular momentum transport was initially assumed to occur through molecular viscosity. However, observations indicate that molecular viscosity alone is too weak to account for the high accretion rates observed in protostellar disks. This suggests that an additional mechanism is necessary for effective angular momentum transport. In accretion disks, turbulence is thought to dominate over molecular viscosity as the primary driver of angular momentum transport. The effective viscosity is often described using the α-prescription introduced by Shakura and Sunyaev (1973) [2]: ν = αcsH. (2.14) Where H is the scale height of the disk and cs is the sound speed of the gas, given by: cs = √√√√kBTg µmp . (2.15) Here, kB is the Boltzmann constant, Tg is the gas temperature, µ is the mean molecular weight of the gas, and mp is the proton mass. The viscosity in the α-disk model is constrained by turbulence, indicating that α reflects the efficiency of angular momentum transport outward by turbulent eddies rather than through molecular interactions. The typical value of α ranges from 10−3 to 10−1, depending on the properties of the disk and the physical processes driving the turbulence. Applying the α-viscosity model, the angular momentum conservation equation 2.9 is modified 6 2. Theory to account for an additional loss of angular momentum. The value of α directly impacts the efficiency of angular momentum transport outward, which regulates the inward motion of gas and the overall accretion rate. The radial velocity for the disk is, vR = − 3 ΣR1/2 ∂ ∂R ( νΣR1/2 ) . (2.16) We can use that equation in 2.6, leading to the following result, Ṁ = 2π 3 R1/2 ∂ ∂R ( νΣR1/2 ) (2.17) where we substitute the viscosity term, Ṁ = −π 6 R1/2 ∂ ∂R ( αcsHΣR1/2 ) (2.18) From here, we can use the scale height H = c2 s√ GM∗/R3 where we apply the Keplerian rotation speed and substitute this into equation 2.18, Ṁ = 6πα c2 s ΩK ∂ ∂R ΣR (2.19) This provides us with a mass accretion expression regarding the α parameter and the surface density Σ. 2.1.3 Turbulent and Rotational Dynamics of Accretion Disks The turbulence in the disk can be measured by comparing azimuthal velocity, vϕ, with the Keplerian orbital velocity, vturbulence = (|vϕ| − |vkep|)/vkep (2.20) For a disk with no turbulence, the difference between the two is zero, while if vturbulence > 1 we have a turbulent disk. 2.1.4 Rotational and Keplerian Timescales Key timescales influence the dynamics of both accretion disks and jets. The star’s rotational timescale is expressed as: trot = 2π δ √ GM∗ R3 ∗ (2.21) where M∗ represents the stellar mass, R∗ denotes the stellar radius, and δ is a dimensionless scaling factor. At a distance of 1 AU, the Keplerian timescale is described as: tkep = 2π√ GM∗ r3 (2.22) In this thesis, we find trot ≈ 15.5 days and tkep ≈ 129 days. To thoroughly analyze the protostar, we aim to simulate at least 2× tkep to accurately capture all behaviors and minimize the impact of initial conditions in the early settings. 7 2. Theory 2.1.5 Magnetic Field Thresholds and Disk Truncation The relationship between the magnetic fields and the ram pressure of the accreting gas explains one theory of how the protostar will accrete mass from its disk. This balance can be expressed as: B2 ∗ 8π ≥ ρv2, (2.23) where B∗ is the magnetic field strength, ρ is the density of the accreting material, and v is the velocity of the infalling gas. The truncation radius, at which this balance is achieved, depends on the protostellar properties and accretion flow dynamics [3] [4]. Supposing the magnetic pressure exceeds the ram pressure of the inflowing material, the inner disk will be truncated, meaning the material will follow the magnetic field lines from the truncation radius into the protostar, as illustrated in the HADES simulation for a low-mass protostar [4]. The magnetic truncation radius can be estimated by substituting parameters for the density and velocity based on mass accretion rates and protostar properties. The equation used in this study is: B∗ ≈ 1 kG × ( fv 1.0 )1/2 ( facc 0.1 )−1/2 ( R∗ 3R⊙ )−5/4 ( Ṁ∗ 10−6M⊙ yr−1 )1/2 ( M∗ M⊙ )1/4 , (2.24) The dimensionless factor fv is the ratio of the radial velocity to the free-fall velocity, fv = vr/vff . The factor facc represents the surface area filling fraction of accretion onto the star and is set to facc = 0.2 [5, ?],[4]. At the same time, fv is shown in early simulations to range between 0.1−0.5. For the high-mass protostar in this work, we use M∗ = 8M⊙, R∗ = 33.4R⊙, and Ṁ∗ = 2.67 × 10−4M⊙ yr−1 [1]. With these values and equation 2.24, one can approximate that the magnetic field strength should be in the interval of 307 − 686 G, depending on the choice of fv for truncation to occur. However, if the mass accretion is ten times higher, the interval would increase to 970 − 2169 G. The HADES simulations [4] for a low-mass protostar show that for B∗ > 1 kG, accretion flows are redirected along magnetic field lines, forming magnetospheric funnels that act as highways for accreting material. At B∗ ≳ 2 kG, accretion predominantly occurs via these funnels, reducing the disk’s inner accretion surface. This is expected for the high-mass protostar, as the accretion rate and magnetic field strength are significantly higher. The balance explained above is achieved at a specific radius called the truncation radius (RT ). At this distance, the magnetic and ram pressure plus the thermal pressure of the accreting material are balanced. The disk will be truncated inside this radius due to the magnetic pressure dominating the ram pressure, which creates the funnel structure. The truncation radius is derived from the fact that magnetic field strength decreases radially as: B(R) ∼ B∗ ( R∗ R )3 . (2.25) By substituting this scaling into the balance equation, the truncation radius RT can 8 2. Theory be expressed as: RT ∼ ( B2 ∗R6 ∗ Ṁ∗ √ 2GM∗ )2/7 , (2.26) where Ṁ∗ is the mass accretion rate and G is the gravitational constant. The value of B∗ will be varied in this thesis to find when the disk gets truncated. The critical magnetic field equation, 2.26, provides a rough estimate of the critical magnetic field strength (Bcrit) needed to truncate the disk. If we use this result, we expect the truncation radius for B = 20 G to be RT = 0.0539 AU from the protostar and roughly RT ≈ 0.44 AU for B∗ = 2000 G. 9 2. Theory 10 3 Methods The simulation in this thesis uses the HADES code that builds on the PLUTO code. PLUTO is a widely used code to solve magnetohydrodynamics (MHD) equations described in the theory section. PLUTO uses advanced numerical schemes such as high-order Runge-Kutta integration and divergence cleaning techniques, which ensure that ∆ · B⃗ = 0 everywhere. It also handles physics like Ohmic dissipation and the Hall effect, which are important for accurately modeling resistive processes in the accretion disk. HADES includes a more detailed treatment of Ohmic resistivity and the Hall ef- fect, which are crucial for modeling the non-ideal effects in accretion disks around protostars [4]. 3.1 Simulation Setup As we use MHD to simulate the behavior of a protostar, the number of cells directly impacts the computational cost and, consequently, the simulation time. There are several ways to divide the area of interest into cells. One method is to divide the area with cells of equal length, meaning a cell near the protostar will be the same size at the outer boundary. A uniform grid would provide consistent resolution across the domain but could be costly regarding computational resources. To address this, one can instead employ radially logarithmically spaced cells that increase in size as one moves further from the protostar. This approach captures all the essential physics occurring near the protostar, where most significant events happen, while still rep- resenting some dynamics as we move outward. Logarithmic grids are standard in accretion physics simulations, representing dynamics across various scales. Both the HADES and the work of [6] are using this setup in their papers. However, one must still use a small enough spacing; otherwise, one will get a lower resolution, impacting the result. This can be seen in [6] where they did a resolution convergence and found that the jet-outflows changed shape in the different resolutions. The simulations grid system is in spherical coordinates, (r, θ, ϕ). To reduce compu- tational cost, we assume axisymmetry around the protostar. This means that there are no cells in the ϕ direction, making the grid a two-dimensional (r, θ), but still solving for all three vector components of velocity and magnetic fields. This setup is commonly called a 2.5D MHD simulation, where the grid is two-dimensional but retains three-dimensional vector components for velocity and magnetic fields. [4] [6] [7] 11 3. Methods 3.1.1 Changes from HADES The same code was used in the HADES simulation as this thesis did, except for two exemptions. The first change involved the tapering factor rfac, which affects both the resistivity and viscosity in the disk. The change was made as early sim- ulation had a problem with numerical instabilities where regions would have had a sharp discontinuity, leading to the code to crash. Instead, a smooth transition was implemented in the two PLUTO C files, res_eta.c and visc_nu.c; // Original code double rfac = tanh(Rp/rmin); // Updated code double rfac = tanh((rcyl - rtrans) / (corot - rtrans)); Where rcyl is the radius in cylindrical coordinates from center of the protostar, rtrans is the transition radius, where resistivity and viscosity begin the transition and it is defined as rtrans = R∗ + rcorot. rcorot is the corotation radius, which is defined as the radius where the disk has the same rotation speed as the protostar, rcorot = GM∗ Ω2 ∗ 1/3. [4] This modification causes the resistivity to decrease gradually instead of suddenly, leading to a much smoother transition. The second modification was adjusting the rotation speed of the gas, vk. Early simulation had a problem with the gas rotating too fast, which led to the centrifugal forces overcoming the gravitational forces. As a result, gas was ejected outward along the midplane. To address this, the epsilon factor, ϵ in the calculation of the vk. vk = sqrt(CONST_G*mstar/(fabs(rsph)*UNIT_LENGTH))/UNIT_VELOCITY; vkep = sqrt(CONST_G * mstar / (fabs(rcyl)*UNIT_LENGTH)); cs = sqrt(CONST_kB *temp/(1.4 *CONST_mp)); eps = cs / vkep; epsFactor = 1.0 - (5.0 * eps*eps / 2.0); vk *= epsFactor; Before the change, vk was assumed to follow a Keplerian profile which neglects the influence of pressure. With the introduction of the epsilon factor and the sound speed, the effect of thermal pressure is included, and the slowing down of the gas is included. 3.1.2 Heating and cooling In this simulation, two cooling methods were used, the continuum and the atomic line cooling, and the X-ray was just used to heat the gas. The heating is calculated via the flux in each cell via, ϵX = nHFXσX , (3.1) where FX is the X-ray flux and σX is the X-ray photo-absorption cross section[4]. This last term describes the reduction in X-ray intensity as the radiation passes 12 3. Methods through a medium with absorbing material. The flux FX is determined using the equation: FX = LX 4πR2 eσX (3.2) where, LX is the X-ray luminosity measured in erg s−1. 3.1.3 Initial conditions The simulations investigate how resolution, magnetic field strength, and disk surface density affect accretion. All other parameters are kept fixed to ensure consistency across simulations. The disk is assumed to be an α-viscous disk with the initial surface density profile following a truncated power-law distribution.[8] [9], Σ(R) = Σ0 ( R R0 )−γ e−R/Rmax , (3.3) where Σ0 is a normalization factor, R0 is a reference radius, γ is the surface density power-law index, and Rmax represents the disk’s outer edge, beyond which the ex- ponential term causes a rapid decrease in surface density [4]. These values are set to R0 = 10 AU and Rmax = 57 AU. Observations of low-mass protostellar disks suggest a γ = 0.4 − 1 range, with a median around 0.9 [10]. However, direct observational constraints for high-mass protostellar disks, such as those around 8M⊙ protostars, are currently lacking. In this theis γ = 1 assuming that the surface density scaling is similar to that observed in lower-mass systems. This is also the value the HADES simulation used[4]. The density of the disk ρ is obtained through the surface density, ρ = Σ(R)√ 2πh exp ( − z2 2h2 ) (3.4) where the height of the disk is obtained via the sound speed in the gas and the Keplerian rotation rate, h = cs/ΩK . The initial temperature profile of the disk is assumed to follow a power-law distri- bution: T (R) = T0 ( R R0 )−0.5 (3.5) where T0 = 1000 K at the reference radius R0. This temperature is consistent with observational measurements of the disk around the 12M⊙ protostar IRAS 20126 + 4104 [11]. The exponent of −0.5 is derived from radiative transfer models for embedded protostellar disks [12], which describe how disk temperature decreases with radius due to stellar irradiation and radiative cooling. The radial velocity profile is set via, vr = − ṁ 2πΣ(R)R, (3.6) 13 3. Methods where ṁ is the target accretion rate, set to ṁ = 10−5. The disk was initialized by a vector potential to magnetize the disk, but due to the axisymmetry, only the θ component is used, Aϕ = 2B1R1 3 − a r−(a−1)/2, (3.7) where B1 = √ 2p1β1, p1 is the pressure at 1AU and plasma beta β1 = ptherm/p1,mag which is the ratio between the thermal and magnetic pressure at radius 1 AU [4]. The magnetic field from the protostar in the form of a magnetic vector potential for a static dipole is also added with the formulae, Aϕ,∗ = B∗R 3 ∗ r2 , (3.8) where B∗ is the protostellar magnetic field strength that will be varied in this thesis. The disk self-gravitation will be neglected in these simulations, but the protostars gravitational profile is given by, Φ(r) = −Gm∗ r , (3.9) where r is the radius from the protostar. The X-ray flux is set differently for two different regions, first, for regions outside the disk, we assume an optically thin X-ray in the region z/h > 4, FX(r) = LX 4πr2 (3.10) where LX is the X-ray luminosity set at LX = 1030 erg s−1. For regions within the disk, up to the surface, the X-ray flux is instead calculated via, FX(r) = LX 3πR2 exp ( −σX ∫ Zsurface Z n(R, z′)dz′ ) (3.11) where n(R, z) is the number density , σX is the X-ray photo-absorption cross section which is set to σX ≈ 10−22 cm−2 for gas what is dominated by hydrogen [13]. This heat generated by the flux will serve as the inner disk’s floor temperature, preventing it from cooling down excessively. The ionization fraction resulting from X-ray irradiation is computed throughout the simulation domain, assuming equilibrium between thermal and X-ray ionization processes and recombination. This ratio is given by, xζ e = √ ζ βn (3.12) where xζ e represents the ionization fraction, ζ is the X-ray ionization rate, β is the recombination coefficient, and n is the local gas number density [4]. The X-ray ionization rate is computed from the flux ζ = FXσX ⟨Eion⟩ , (3.13) 14 3. Methods where Fx is the x-ray flux, σX = 10−22 is the photo-absorption cross-section and ⟨Eion⟩ ≈ 37.0 eV, is the average energy required to ionize hydrogen and eject an electron. Additionally, the disk is ionized by thermal ionization, with the ionization fraction expressed as xζ e ≈ 6.47 × 10−13 ( xK 10−7 )( T 103 )0.75 (2.4 × 1015 cm−3 n )( e−(25188 K/T ) 1.15 × 10−11 ) (3.14) where xK is the amount of available potassium in the disk, set to xK = 10−7 [14] [4]. Finally, by combining these two components, we obtain the total ionization factor xe = xζ e + xζ e. The protostar has a fixed rotation speed of Ω∗ = fbu √ Gm∗ R3 ∗ , (3.15) where fbu denotes the breakup speed fraction, assumed to be 0.5, aligning with the rapid rotation observed in high-mass protostars [15]. For our star, this condition is Ω∗ ≈ 4.69 × 10−6 rad s−1. The initial values are summarized and presented here. The initial magnetic strength of the protostar varies between B∗ = 500, 1000, 2000, and 4000 G, and the initial surface density, Σ0 = 200, 500, 1000, and 2000 g cm−2. The rest of the parameters will be set. The mass of the protostar, m∗ = 8M⊙, the radius, r∗ = 33R⊙, the temperature T0 = 1000 K, the initial density for the envelope ρ = 107 g cm−3, the viscosity parameter, α = 0.1, and the X-ray luminosity LX = 1030 erg s−1. 3.1.4 Boundary condition The inner boundary is located just outside the surface of the protostar and is mod- eled as a perfect conductor. Since a perfect conductor cannot sustain an induced electric field, the protostar does not generate an electric field, nor can one be in- duced within it. Consequently, the electric field at the boundary must also be zero to remain consistent with this assumption. This is ensured by requiring that the electric field produced by the rotating magnetic protostar vanishes at the boundary. Such a condition can only be achieved if the infalling material co-rotates with the protostar [7]. E⃗ = B⃗ × (u⃗ − Ω⃗∗ × R⃗) = 0, (3.16) whereΩ∗ is the angular speed of the protostar, which rotates around θ = 0. We also ensure that the toroidal magnetic field, RBϕ, is conserved at the boundary. To fully uphold the perfect conductor condition, the electric field component in the θ direction must also be zero, Eθ = 0. This prevents numerical artifacts in the magnetic field caused by the rotation of the magnetized protostar while ensuring that the infalling material is absorbed correctly and that no material bounces back. At the outer boundary, fixed at 57.0 AU in all simulations, an "outflow-only" condition is applied. This setup allows material to exit the computational domain freely but prevents re-entry. 15 3. Methods 3.1.5 Parameter Variations This section will cover the three parameters that will be varied in this thesis: initial stellar magnetic field strength (B∗), the initial surface density of the disk (Σo), and the different grid resolutions. 3.1.5.1 Stellar Magnetic Field Strength Several simulations were performed with varying stellar surface magnetic field strengths to investigate the role of magnetic fields in the accretion process and disk trunca- tion. The selected values are 50 G, 500 G, 1000 G, 2000 G and 4000 G to explore a wide range of possible magnetic fields, since the magnetic field of massive proto- stars is not well constrained. Based on theoretical predictions and previous MHD simulations, magnetospheric accretion is expected to start dominate within the field strength range of 1600 − 2400 G, where the magnetic pressure can counteract the ram pressure of the accretion flow, potentially altering the accretion geometry. 3.1.5.2 Surface Density The surface density, Σo, directly affects the accretion disk’s total available mass and the protostar’s corresponding mass accretion rate. The Σo was varied for the four cases: 200, 500, 1000, and 2000 g cm−2 to investigate its influence on the disk structure and accretion efficiency. These values are chosen to approximate a mass accretion rate of Ṁ = 2.67 × 10−4M⊙yr−1. 3.1.6 Grid Since resolution will impact the numerical accuracy, a series of tests are performed to evaluate how different grid spacings impact the results. The simulations are run at three resolutions: 128 × 128, 256 × 256, and 512 × 512 cells. These cells are placed in the interval of [10−2, π − 10−2] radians.The small offset is done to avoid the poles at θ = 0, π. The three grid configurations yield spatial resolutions of approximately ∆r ≈ 5.69 × 10−3, 2.79 × 10−3, 1.38 × 10−3 AU. At the outer boundary, the corresponding resolutions are ∆r ≈ 518×10−3, 264×10−3, 133×10−3 AU. For the rest of the simulation, where other parameters are run, the 512 × 512 cell configuration will be used. 3.2 Computational Methods Each compute node in the DARDEL system consists of two AMD EPYC 64-core processors, totaling 128 physical cores per node [16]. PLUTO partitions the grid out in parallel, so each core is responsible for computing a set of grid cells, ensuring parallel simulation execution. Consequently, for a 128 × 128 grid configuration, 128 cores (i.e., one full node) are utilized, with each core handling approximately an 11×11 patch. For higher-resolution configurations, such as 256×256 and 512×512, the computational requirements scale accordingly to 256 and 512 cores, where each core is responsible for a 16 × 16 and 23 × 23 patch, respectively. 16 3. Methods The simulation evolves until the protostar has completed 259 days of rotation. Since increasing the grid resolution significantly extends the computational runtime, dif- ferent resolutions lead to varying completion times. The 128 × 128 configuration, which has the lowest resolution, completed within a week, whereas the 512 × 512 configuration, which has the highest resolution, required approximately 5 weeks. 3.3 Data analysis The output data from the PLUTO/HADES simulations was stored in VTK format [17]. These files contain 2D arrays of the main physical quantities, such as density, velocity fields, magnetic field components, and temperature. Each simulation gen- erated around 3000 VTK files. To extract raw data from the VTK files, pyPLUTO, a Python tool developed by the creators of PLUTO, was used [17]. The extracted data was then processed using NumPy for numerical operations and visualized with Matplotlib. Together with the software, ParaView [18] was used to generate density maps, velocity streamlines, magnetic field structures, and temperature distributions. When working with large amounts of data, such as processing more than 10 VTK files simultaneously, the analysis was conducted directly on the DARDEL cluster, using parallel processing to speed up data retrieval. The cluster was also used for time-series analysis, examining how different physical parameters evolve, such as accretion rates, disk instabilities, and outflow properties. Additionally, simulation movies were created on the cluster using C code and processed with FFmpeg to help visualize the protostar evolution. 3.3.1 Accretion filling fraction Filling plots were made to understand the fraction of the protostar surface under different conditions, such as where the top 10% of the total mass gets accreted. The cumulative mass fraction is then used to calculate where the 10%, 25%,50%, 75%, and 90% top % is located on the surface. For each snapshot, that is, for each VTK file in the simulation, all cells are sorted along the surface of the protostar with the highest mass flux to the lowest mass flux. The corresponding cumulative sum of the mass accretion and cumulative surface fraction are computed for each index. This is illustrated in figure 3.1 which is done for the last snapshot for the model B∗ = 2000 G, Σo = 2000 gm−2 The filling factor is then obtained for each X%, taking the ratio of ffilling = AX,m Atotal , where the AX,m is the cumulative surface area corresponding on the top X% of the total mass accretion. Atotal is the total surface area of the protostar. 17 3. Methods 0 100 200 300 400 500 Index 0.0 0.2 0.4 0.6 0.8 1.0 Fr ac tio n Mass fraction Area fraction vr m 10 11 10 9 10 7 10 5 10 3 10 1 m (M y r 1 ) | v r (g s 1 ) Figure 3.1: Example of a filling fraction calculation for the snapshot at day 259 in the B∗ = 4000 G, Σo = 2000 g cm−2 model. The red line represents the sorted mass flux at the protostar’s surface, the green line represents the accretion rate, the black line shows the cumulative mass accretion, and the blue line depicts the cumulative area as a function of the sorted cell index. The vertical dashed lines indicate the indices at which the top 10%, 25%,50%, 75%, and 90% of the total mass is accreted. 3.3.2 Truncation radius The truncation radius, RT , is determined for each simulation snapshot by identifying the radial location in the midplane of the disk where the magnetic pressure balances the combined thermal and kinetic pressure of the accreting material, according to equation 3.17, ∆P = Pmag − (Ptherm + Pkin), (3.17) The calculation is performed for each radial cell until the truncation radius is iden- tified, where ∆P = 0. An average within a wedge of eight cells surrounding the disk is taken to reduce noise. This averaging helps mitigate the effects of local fluc- tuations or turbulence that could lead to ∆P = 0, lowering the risk of inaccurate measurements. The truncation radius is also monitored over time for each snapshot. This data will help us understand how the truncation radius changes and obtain a mean truncation value over time. The disk may not always be truncated, and this will also be monitored to determine the fraction of time the disk is truncated across the different models. 18 4 Results This chapter will present the results from the different simulations in this thesis, as described in the method and theory section. The models that were used in this thesis are presented in Table 4.1 Model Name Protostellar B∗ [G] Σo [g/cm2] B50 50 2000 B500 500 2000 B1000 1000 2000 B2000 2000 2000 B4000 4000 2000 Σo200 2000 200 Σo500 2000 500 Σo1000 2000 1000 Σo2000 2000 2000 Table 4.1: Summary of simulation models with varying protostellar magnetic field strength (B∗) and disk surface density (Σo). Both the magnetic field strength and the disk surface density had a direct impact on the mass accretion rate, where the model with B∗ = 2000G, Σo = 2000 had the highest mean accretion rate, while B∗ = 2000, Σo = 200 had the lowest. The B∗ = 4000G, model did have a high mean accretion rate, but the accretion happened in burst, which led to a low accretion between the burst, followed by a peak which led to a mean accretion of 27.39 with a standard deviation of Σo = 138.89 × 10−4 M⊙/year. The disk truncation depends on both magnetic field strength B∗ and disk surface density. For weaker magnetic fields, truncation is rare and occurs only irregularly. However, truncation becomes more frequent and sustained as B∗ increases. This trend is also observed in the simulation series with varying surface densities, where all models experience truncation. The model with the lowest surface dense density had the highest truncation fraction of all models and had the biggest RT radius. 19 4. Results 4.1 Accretion Rate The mass accretion was plotted as a function of time for the different models, to get an overview of how the various models behave. The accretion is the total mass accretion over the surface of the protostar. The change in mass accretion as a time function is visualized in Figure 4.1. Figure 4.1: The figures show the evolution from an initial state of a protostar for four different initial magnetic field strengths of the protostar and a comparison to the 2.4 × 10−4 line which represents the predicted mass accretion for an 8M⊙ protostar [1]. The right image compares changing the initial surface density of the disk, Σo. . 4.1.1 Magnetic Field Variation (B∗) Series The models B∗ = 50, 500, 1000, 2000 G in Figure 4.1 did not experience significant variations in accretion rate, though some fluctuations and a gradual decline over time were observed. In general, the accretion rate increases with B∗ models, but when B∗ is increased to B∗ = 4000 G, the overall accretion rate decreases, and bursts appear around day 90 and 210. The model B∗50 G is the model that is closest to the expected value of 2.4×10−4M⊙/year which a protostar of the size is expected to have [1]. It varies around the target value after day 100 but suddenly drops rapidly at day 210. The mean accretion rate is determined using an average over 209 days, which starts after the transient phase of 50. The transient phase was set after calculating that 50 days roughly corresponds to three rotations of the protostar. The result of this analysis presented in the figure 4.2 and in the table 4.2. It is seen in both Table 4.2 and Figure 4.2, that the B∗ = 50 G model is experiencing the lowest accretion rate, but also, the standard deviation is at the lowest. Its is 20 4. Results Figure 4.2: The figure shows the mean and standard deviation of the mass accre- tion for the two simulation series B∗ Σo. A transient time of 50 days was used to erase the impact of the starting condition. The dots mark the mean accretion for each model, and the lines are the standard deviation. Mean Mass Accretion Rates ⟨ṁ⟩ [M⊙/year] × 10−4 B∗ [G] 50 500 1000 2000 4000 MAIN 4.89 ± 4.19 10.92 ± 8.44 13.66 ± 10.44 39.75 ± 41.22 27.39 ± 138.89 Table 4.2: Mass accretion rates for different magnetic field strengths B∗ with fixed Σo = 2000 g cm−2. however noticeable that compared with the mean rate, the B∗ = 50 G model have higher rate of standard deviation compared with the two model, B∗ = 500 and B∗ = 1000 G models, where the ratio ⟨ṁ⟩/σ(ṁ) drops from 86% for the B∗ = 50 G to 77% and 76%. For the two models, B∗ = 1000 and B∗ = 2000, the standard deviation exceeds the ⟨ṁ⟩, especially for the B∗ = 4000 G, which has an standard deviation of almost three times the B∗ = 2000 G model. The ⟨ṁ⟩ is observed to be increasing as B∗, but this is only true for the models B∗ ≤ 2000. For the B∗ = 4000 model, it is instead observed that the mean accretion has dropped by almost 30% compared with the 2000 G model, 4.1.2 Surface Density Variation (Σo) Series The models Σo = 200, 500, 1000 g cm−2 in figure 4.1 are not stable as time progress. All but Σo = 2000 g cm−2 is experiencing large variation in the accretion and where the Σo = 200 g cm−2 model is overall the model with both the lowest rates between the peaks, but also the peaks itself. The phase between peaks and the peaks increases in magnitude as Σo does. The mean accretion and standard deviation are visualized in Figure 4.2, but also 21 4. Results quantified in Table 4.3, where the same procedure with the use of a transient time as for the B∗ models was done. B∗ = 2000 G, Mean Mass Accretion Rates ⟨ṁ⟩ [M⊙/year] × 10−4 Σo [g cm−2] 200 500 1000 2000 SIGC 3.54 ± 18.01 6.34 ± 32.72 21.34 ± 61.94 39.75 ± 41.22 Table 4.3: Mass accretion rates for B∗ = 2000 G with varying disk surface densities Σo. For the variation of Σo, it is clear that the disk with the highest surface density, Σo = 2000 g cm−2, also had the highest mean mass accretion rate of ṁ = 39.75 × 10−4, M⊙/year. There is a clear trend where the accretion rate decreases as the disk becomes thinner. While the standard deviation, generally decreases with a lower surface density, the Σo = 1000 g cm−2 model experienced the largest variability, with a standard deviation of σ = ±61.94 × 10−4, M⊙/year. This is almost 1.5 times higher than in the Σo = 2000 g cm−2 model and roughly twice as large as in the Σo = 500 g cm−2 model. 4.2 Mass Accretion Bursts This section presents the results on the accretion rates for the different models. For the B∗ series, the accretion rate remained relatively stable for the four lower field strengths. However, as previously demonstrated in earlier chapters, the B∗ = 4000 G model exhibits a trend where accretion occurs primarily in bursts rather than steadily. A similar pattern is observed in the Σo series, where accretion variability increases with higher disk surface densities. 4.2.1 Magnetic Field Variation (B∗) Series The mass accretion is visualized in a histogram where the y axis is the mass accreted, and in the x-axis the mass accretion rate normalized by the mean accretion rate, ṁ/⟨ṁ⟩ is presented. A cumulative distribution function (CDF) plot is also made for easier understanding of which rates most of the mass is accreted to. For the four lower B∗ values, the accreted mass is centered around the mean accretion rate, indicating that these models are relatively stable over time. The accreted mass is lowest for B∗ = 50 G, with increasing accreted mass for higher B∗ values. The B∗ = 2000 G model has two peaks, which could indicate that mass accretion is less stable over time and may experience some bursts. The B∗ = 4000 G model follows an entirely different pattern than the others. Nearly all the accreted mass corresponds to accretion rates above the mean. Specifically, a smaller peak is observed around 15 × ⟨ṁ⟩, followed by more significant peaks at 50 × ⟨ṁ⟩. For the three models with B∗ = 500, 1000, 2000 G, the cumulative mass fraction is similar, with approximately 60% of the total accreted mass originating from accretion rates above the mean. However, the B∗ = 50 G and B∗ = 4000 G models show a significantly higher 22 4. Results Figure 4.3: Probability density function (PDF) and cumulative distribution function (CDF) of mass accretion rates normalized by the mean accretion rate, ṁ/⟨ṁ⟩mean, for a transient time of 50 days. The left panel (PDF) shows the distri- bution of accreted mass as a function of the normalized accretion rate for different initial magnetic field strengths of the protostar (B∗). The right panel (CDF) presents the cumulative mass fraction accreted as a function of the normalized accretion rate. The thick lines in the CDF represent the contribution from accretion rates above the mean rate, emphasizing the fraction of mass accreted at higher rates. fraction of accreted mass at rates above the mean, with values of 80% and 90%, respectively. The fraction of total accreted mass for the different models is also plotted in a box plot for easier visualization of how accretion bursts contribute to the overall accretion process. B∗ (G) ≤ Mean >1.3× Mean >3× Mean >6× Mean 50 0.2219 0.6806 0.0623 0.0195 500 0.3735 0.4841 0.1393 0.0000 1000 0.3844 0.5027 0.0745 0.0078 2000 0.4113 0.5098 0.2769 0.0414 4000 0.0842 0.8823 0.7892 0.6909 Table 4.4: Fraction of total accreted mass occurring above specific accretion rate thresholds for different B∗ values. Figure 4.4 shows the fraction of the total accreted mass occurring at different accre- tion rate thresholds for each magnetic field strength, B∗. Most of the total accreted mass comes from accretion rates above the mean accretion rate for all models. For B∗ = 50 G, approximately 22% of the total accreted mass comes from accretion rates 23 4. Results Figure 4.4: A box plot for the different initial protostar magnetic field strengths, B∗. The figure shows how much of the accreted mass occurs below and above the mean, and then after 3 × ⟨ṁ⟩ and 6 × ⟨ṁ⟩ for each model. at or below the mean, while the remaining mass is primarily distributed around slightly higher accretion rates. For the B∗ = 500, 1000, and 2000 G models, the fraction of mass accreted above the mean is similar, with most of the mass falling within the 1.3 − 3× mean range, and only a small fraction exceeding 6× the mean. What stands out in the figure is that the B∗ = 4000 G model accretes nearly 70% of its total mass from accretion rates above six times the mean accretion rate. In contrast, this fraction remains below 5% for the other models. 4.2.2 Surface Density Variation (Σo) Series Figure 4.5 shows the accretion rate distribution for varied Σo, similar to Figure 4.3 Figure 4.5 shows how the total accreted mass and its distribution shift as the surface density decreases. As Σo decreases, the total accreted mass also decreases, but more importantly, the accreted mass shifts toward higher accretion rates relative to the mean. For Σo = 200, a peak is observed around 30 × ⟨ṁ⟩, while for Σo = 500, this peak is located at 10 × ⟨ṁ⟩. However, unlike Σo = 200, the Σo = 500 model also shows significant accretion at even higher rates beyond its peak. The Σo = 1000 model has one large peak at roughly 20×⟨ṁ⟩ and a minor peak near ⟨ṁ⟩. The Σo = 2000 model exhibits two large peaks, both more centered around the mean accretion rate. The general trend indicates that lower surface density disks accrete more mass at higher rates, whereas thicker disks have their accretion more concentrated around the mean. 24 4. Results Figure 4.5: Probability density function (PDF) and cumulative distribution function (CDF) of mass accretion rates normalized by the mean accretion rate, ṁ/⟨ṁ⟩mean, for a transient time of 50 days. The left panel (PDF) shows the distri- bution of accreted mass as a function of the normalized accretion rate for different surface density values (Σ0). The right panel (CDF) presents the cumulative mass fraction accreted as a function of the normalized accretion rate. The right panel of Figure 4.5 confirms these observations. As the disk transitions from Σo = 2000 to Σo = 200, the total fraction of mass accreted above the mean rate increases from approximately 60% to 95%. For the two intermediate models, Σo = 1000 and Σo = 500, this fraction is around 81% and 95%, respectively. Exact values are presented in Table 4.5. The accretion behavior is further analyzed in Figure 4.6. Σo (g cm−2) ≤ Mean >1.3× Mean >3× Mean >6× Mean 200 0.0478 0.9449 0.9247 0.9123 500 0.0428 0.9515 0.9150 0.8193 1000 0.1843 0.7419 0.5708 0.4516 2000 0.4113 0.5098 0.2769 0.0414 Table 4.5: Fraction of total accreted mass occurring above specific accretion rate thresholds for different Σo values. The bar plot and table confirm that disks with lower surface densities accrete most of their mass at rates above the mean accretion rate. Additionally, the fraction of accreted mass occurring at rates exceeding 6 × ⟨ṁ⟩ is 91%, 82%, and 45% for the Σo = 200, Σo = 500, and Σo = 1000 models, respectively. For the Σo = 2000 model, this fraction drops to about 4%, which is a difference of 41 percentage points compared to the Σo = 1000 model. 25 4. Results Figure 4.6: A box plot representing the various initial surface densities, Σo. The figure illustrates the amount of accreted mass that occurs below the mean, and subsequently after 3 × ⟨ṁ⟩ and 6 × ⟨ṁ⟩ for each model. 4.3 Location of accretion on the protostellar sur- face Where the accretion happens on the protostar will be determined by several different aspects, but one main drive is how strong the magnetic field is compared to gravity. To easier understand what is going on, the accretion momentum flux, ρvr was plotted as a function of time and angle, and can be seen in Figure 4.7. 4.3.1 Magnetic Field Variation (B∗) Series The accretion momentum flux is presented in figure 4.7. The results in Figure 4.7 illustrate that for the lower B∗ models, most of the mass is accreted near the mid- plane of the protostar. However, as B∗ increases, accretion shifts to higher latitudes, with the B∗ = 4000 G model showing accretion concentrated around 45◦. Addition- ally, this model exhibits episodic accretion bursts at t = 90 and t = 210 days. The morphology of accretion also changes with B∗. At B∗ = 50 G, accretion is confined to the midplane as boundary-layer accretion. As B∗ increases to 500 and 1000 G, material flows toward high latitudes and the equatorial region. At B∗ = 2000 G, accretion is increasingly dominated by high-latitude flows, with periodic equatorial accretion events. Finally, in the B∗ = 4000 G model, accretion occurs in well- defined high-latitude streams, while equatorial accretion is minimal and primarily 26 4. Results Figure 4.7: Accreting rate versus time and latitude on the stellar surface for four simulations with different initial disk densities, B∗, as annotated in each sub-figure. The color coding represents the magnitude of the accreting momentum flux, with darker shades indicating lower flux and lighter shades indicating higher flux. burst-driven. To more easily quantify the results in Figures 4.7 and 4.9, two filling factor plots were made, one for each parameter scan, where the amount of accretion over a specific area fraction is visualized. 27 4. Results Figure 4.8: Accretion rate and surface area coverage for the B∗-series. The top panel shows the accretion rate for different magnetic field strengths over time. The middle and bottom panels display the fraction of the protostar surface responsible for 90% and 50% of the total accreted mass, respectively. All models have a high variation in the area fraction as time progresses. The four models,B∗ = 50, 500, 1000 and B∗ = 2000 G, is seen to vary between around (0% − 35%) for the top 90% of the total mass accreted, and around (0% − 10%) for the top 50%. While for the B∗ = 4000, the variance is much larger (2% − 40%) for the top 90% and it reaches almost to 20% for the top 50%. A more detailed analysis is presented in Table 4.6, which presents each fraction’s mean values over time. Note, for the B∗ = 4000 G, model, the area fraction peaks in both 50 and 90% just before the mass accretion peak. At the time of the peaks of mass accretion, the area fraction is below 5% of the surface. The results from figure 4.8 are quantified and presented in the table 4.6. For B∗ = 50 G, the top 10% of the accreted is distributed over 0.67% of the protostar surface , whereas for B∗ = 4000 G, this value decreases to 0.44% of the surface. This trend clearly shows thatB∗ increases, accretion becomes more confined to a smaller area. However, the model B∗ = 500 has the broadest accretion surface area in the top 90% where almost 15% of the surface is used while for B∗ = 50 and B∗ = 4000 G, 28 4. Results Table 4.6: Filling fraction (percent) of the accretion flow onto the protostar ac- cording to the top X% of the mass accretion for B∗ (G). B∗ (G) 10% 25% 50% 75% 90% 50 0.67 1.52 3.39 6.98 12.76 500 0.55 1.23 3.30 7.94 15.01 1000 0.67 1.52 3.39 6.98 12.76 2000 0.49 1.03 2.38 5.16 9.22 4000 0.44 0.88 1.85 3.71 6.66 this value drops to 12.76% and 6.66% of the surface. 29 4. Results 4.3.2 Surface Density Variation (Σo) Series The same analysis for the B-series was done for the Σo-series with the result seen in figure 4.9. The Σo-series shows a similar behavior to the B-series, but as Σo Figure 4.9: Accretion rate versus time and latitude on the stellar surface for four simulations with different initial disk densities, Σo, as annotated in each sub-figure. The color coding represents the magnitude of the accreting momentum flux, with darker shades indicating lower flux and lighter shades indicating higher flux. decreases, the accretion is mostly concentrated in an accretion burst occurring at t ≈ 150 days for Σo = 200. As Σo increases, the accretion exhibits stronger bursts, but the accretion angle remains large compared to the results in figure 4.7. For the Σo series, the same behavior in the peak in area fraction just before the peak in mass accretion rate, is observed in for Σo = 1000 simulation and also in the Σo = 500 but not as distinct. For the lowest density, Σo = 200, the area fraction for both 50 and 90% instead peaks between the peak mass accretion rate, and its most visible in the 50%. As soon as the mass accretion rate increases, the area fraction decreases which indicated that the accretion burst is occurring on small areas of the protostar. 30 4. Results Figure 4.10: Accretion rate and surface area coverage for the Σo series. The top panel shows the accretion rate over time for different magnetic field strengths. The middle and bottom panels display the fraction of the protostar surface responsible for 90% and 50% of the total accreted mass, respectively. The results from figure, 4.10 is quantified and presented in the table 4.7. Table 4.7: Filling fraction (percent) of the accretion flow onto the protostar ac- cording to the top X% of the mass accretion for Σo (g cm−2). Σo (g cm−2) 10% 25% 50% 75% 90% 200 0.55 1.05 2.43 6.19 11.86 500 0.57 1.12 2.42 5.15 9.50 1000 0.45 0.85 1.82 3.92 7.96 2000 0.49 1.03 2.38 5.16 9.22 From Table 4.7, it is seen that the two lowest disk densities have a much broader accretion region across all five percentile levels compared to the denser disk models. This suggests that lower-density disks allow for accretion to happen over a larger surface area of the protostar. However, one noticeable detail is that the Σo = 1000 g cm−2 model has the most localized accretion, with the smallest surface areas 31 4. Results responsible for the accreted mass. Additionally, the Σo = 1000 gcm−2 model is quite similar to the Σo = 500 g cm−2 case in the top 50%, 75%, and 90% ranges, where the values differ only by a few percentage points. These values deviate by less than a percentage point in most cases compared to the B∗ series at all fractions, and the lowest surface density is very close to the values in B∗. 4.4 Truncation Radius The truncation radius for both B∗, Σo series is presented through figures for vi- sualization and quantified using the mean truncation radius. Since the disk is not always truncated in all models, the mean value is calculated only from the instances where truncation actually occurs. The truncation time fraction is also measured. 4.4.1 Magnetic Field Variation (B∗) Series The truncation radius was plotted against radius and can be seen in figure 4.11, Figure 4.11: Figure visualizes the truncation radius as a function of time. The disk is truncated where the black dots are for the B∗ series. The red dotted line represents the mean truncation radius for each model. The figure shows that as B∗ increases, so does the truncation radius. More im- portantly, the frequency of disk truncation also increases, meaning that stronger magnetic fields push the truncation radius outward and cause the disk to be trun- cated more frequently. This behavior is quantified in Table 4.8, which presents the mean truncation radius and the fraction of time the disk remains truncated. Table 4.8: Truncation Radius and Time Fraction for B∗ Model B∗ [G] Mean RT [AU] Truncation Time Fraction % 50 0.1558 0.0393 500 0.1560 0.0237 1000 0.1610 0.1849 2000 0.1641 0.3655 4000 0.2131 0.8790 The table 4.8 confirms the trends observed in figure 4.11. As B∗ increases, the mean truncation radius, ⟨RT ⟩, and the fraction of time the disk remains truncated show a 32 4. Results clear increase. The ⟨RT ⟩ remains almost constant between B∗ = 50 G and B∗ = 500 G, with only a 0.13% increase. Between B∗ = 500 G and B∗ = 1000 G, there is a slightly larger 3.2% increase, followed by a 1.93% increase from B∗ = 1000 G to B∗ = 2000 G. However, in the B∗ = 4000 G model, the truncation radius increases significantly by approximately 30% compared to B∗ = 2000 G. The most significant change occurs in the truncation time fraction. In the B∗ = 50 G model, the disk is truncated only 3.9% of the total time, whereas for B∗ = 2000 G, this value increases to 36%. In the B∗ = 4000 G model, the truncation time fraction rises to almost 88%. 4.4.2 Surface Density Variation (Σo) Series The truncation radius was plotted against radius and can be seen in Figure 4.12, Figure 4.12: Figure visualizes the truncation radius as a function of time. The disk is truncated where the black dots are for the Σo series. In the figure 4.12 it is clear that the lower surface density of the disk is experiencing variation which decreases as the density increases. The same happens with the mean ⟨RT ⟩. Table 4.9 summarizes the finding for the Σo series, Table 4.9: Truncation Radius and Time Fraction for Σo Model Σo [g/cm2] Mean RT [AU] Truncation Time Fraction % 200 0.2434 0.8864 500 0.2051 0.8664 1000 0.1821 0.7241 2000 0.1641 0.3655 The table 4.9 confirms what is observed in figure 4.12. Both the mean truncation radius, ⟨RT ⟩, and the truncation time fraction decrease as the surface density in- creases. From Σo = 200 to Σo = 500 g cm−2, ⟨RT ⟩ decreases by 15.74%. The decline continues with an 11.21% decrease from Σo = 500 to Σo = 1000 g cm−2, followed by a 9.88% reduction from Σo = 1000 to Σo = 2000 g cm−2. The greatest change in truncation time fraction occurs between Σo = 1000 and Σo = 2000 g cm−2, where it drops by nearly 50%. While between Σo = 500 and Σo = 1000 g cm−2 the change is only 16.42%. 33 4. Results 4.5 Disk Dynamics To gain a better understanding of the disk’s physical properties and overall behavior, we several properties such as temperature, density, entropy, and pressure. These properties help characterize the thermodynamic state and structure of the disk. The results of this analysis are presented in this section. 4.5.1 Magnetic Field Variation (B∗) Series The temperature, density, entropy, and pressure profiles were obtained by averaging values within a wedge of 16 cells around the midplane. The radial domain was divided into 32 bins, and the mean values were computed within each bin to create a representative profile. The resulting are presented in figure 4.13. Figure 4.13: Temperature, density, entropy, and pressure profiles for the B∗ series, averaged over a midplane wedge and binned radially into 32 bins. The solid lines represent the profile at t = 259 days of simulation, while the dotted lines correspond to t = 16 days. The figures 4.13 shows that all properties is converging in the outer regions of the disk, regardless of the initial magnetic strength B∗, or what snapshot that was used. Temperature, pressure, and density is experiencing the same trend wheres they are all at the highest values near the protostar, and then as the radial distance increases, so does the quantities decrease. For B∗ = 50, 500, 1000 and 2000 G, the trend is the same with only minor variations in magnitude. The B∗ = 4000 G model stands out in the temperature and density where it gets very low quantities near the protostar followed by large variations for the t = 259 day snapshot, before converging like the rest in the outer regions. Note that the B∗ = 4000 G model also have huge difference in entropy between the two snapshots. 34 4. Results 4.5.1.1 Radial and Azimuthal Flow Variations The velocity of the disk is crucial for analyzing its structure and stability, as it pro- vides insight into the underlying dynamics and turbulence within the system. To better understand these effects, the velocity is decomposed into its three compo- nents, vr, vθ, and vϕ, where the Keplerian velocity is subtracted from the azimuthal component. The velocity profiles for these components are illustrated in figure 4.14, Figure 4.14: The radial, poloidal, and azimuthal velocity components for two snapshots, day 16 and day 259, for the B∗ series. The Keplerian velocity is subtracted from vϕ. From the figure 4.14, it is seen that as B∗ increases in strength, so does all three velocity components. The general trend here is that the velocity for all three com- ponents is decreasing as the radial distance increases. However, it is noticeable that the vϕ component for B∗ = 50, at t = 16 days, first hits a dip in the inner regions and after that increases, while at day t = 259, it instead hits a bump and then con- verges to the same velocity as at t = 16 days. The three velocities are roughly of the same magnitude in the inner region for all models, but the vθ is rapidly decreasing in the outer region and ends up in several of magnitudes lower then the other two components. For the B∗ = 4000 G model, the variations in all three components are strongest, and the B∗ = 50 G model have the least variation. To further analyze the velocity structure within the disk, figure 4.15 illustrates the variation of the poloidal velocity as a function of radius. 35 4. Results Figure 4.15: The figure shows the variation of poloidal velocity as a function of radius in the disk for different magnetic field strengths B∗. The left figure is the mean polodial velocity, while the right one shows the standard deviation of the same velocity. t = 259 days. For all models, the poloidal velocity is generally highest in the inner regions of the disk and decreases with increasing radius. The velocity increases as the magnetic field strength B∗ increases, where the strongest, B∗ = 4000 G, generates a mean poloidal velocity nearly 100 times higher than in the B∗ = 50 G and B∗ = 500 G models. Even though the B∗ = 4000 G model has very high velocities near the protostar, it decreases with radius, and in the outer regions, it is only slightly higher than the B∗ = 2000 G model. For the intermediate values of B∗, the B∗ = 1000 G model has a velocity profile between the lower and higher field models, showing a steady decline without the extreme peaks seen in the B∗ = 2000 G and B∗ = 4000 G models. The B∗ = 2000 G model follows a similar trend to B∗ = 4000 G but at a lower level.The standard deviation of the poloidal velocity follows the same general trend as the mean velocity, with the strongest magnetic field models having the highest variation, and the weakest have the least variation. It is also shown that even if the B∗ = 4000 G has the largest variations in the inner region, the B∗ = 2000 G model has the largest variation in the most outer region. To analyze the turbulent motion of the disk, the total velocity dispersion and mean velocity at different snapshots are plotted , along with the dispersion of the vϕ −vkep component normalized with the speed of sound and can be seen in figure 4.16 It is seen in the figure that the B=2000 G and B∗ = 4000 G models have higher velocity dispersion in the inner regions, but the differences between models become smaller further out in the disk. At day 16, the variations are small, but at later times, the velocity dispersion in the inner regions becomes more different between the models. The B∗ = 4000 G model has larger variations but does not always have the highest values at all radial distance. In the outer regions, both the velocity dispersion and mean velocity become similar for all models and times. In the right panel of Figure 4.16, the differences between the models are larger and increase over time. At its peak (t = 234 days), the B∗ = 4000 G model reaches about 100 × cs, 36 4. Results Figure 4.16: The figure shows the total velocity dispersion and mean velocity profiles at different times. The left panels represent the total velocity dispersion σ(Vtotal) (solid lines) and the mean velocity ⟨Vtotal⟩ (dotted lines). The right sub figures shows the dispersion of the velocity component, vphi −vkep normalized by the sound speed. 37 4. Results while the B∗ = 50 G model at the same radius is around 0.1× cs. Also, even though there are large differences in the inner disk between the models,they all converge in the outer region to the same value. The ratio between the radial and free fall velocity is plotted in figure 4.17 Figure 4.17: The figure shows the radial velocity normalized by the free-fall ve- locity, vr/vff , as a function of radius for different magnetic field strengths B∗ and time. From the figure, an increasing B∗ leads to greater fluctuations in the radial velocity, especially in the inner regions of the disk, while all models tend to stabilize at larger radii. The peak of the fluctuations for the lower magnetic field models, B∗ = 50 G and B∗ = 500 G, is centered further out, whereas for the stronger field models, the 38 4. Results fluctuations peak much closer to the protostar. 4.5.1.2 Thermal, Magnetic, and Kinetic Pressure The pressure in these simulations consists of three components that influence the disk evolution, thermal pressure (Ptherm), magnetic pressure (Pmag), and kinetic pressure (Pkin). These three components are plotted for each B∗ model at two different snapshots, resulting in Figure 4.18. Figure 4.18: The figure shows the variation of the thermal, magnetic, and kinetic pressure components with radius for two snapshots: day 16 and the final snapshot at day 259 for the B∗ series. For the lower magnetic field models, B∗ = 50, 500, and 1000 G, the thermal pressure dominates across the disk, followed by the kinetic pressure, while the magnetic pressure remains the weakest. The magnitude of the thermal pressure remains similar for all models and times, except for the B∗ = 4000 G case, where both the magnitude and shape change significantly in the inner region at day 259. At this time, the magnetic pressure becomes the dominant component near the protostar. In general a clear trend is observed in the pressure profiles, both the magnetic and thermal pressure increase as B∗ increases, while the kinetic pressure remains roughly the same if not slightly weaker. 4.5.1.3 Magnetic and Plasma The magnetic field strength and plasma beta (β) were analyzed as a function of radius and time for the different B∗ models. The results are presented in Figure 4.19. At the protostar surface, the magnetic field strength matches each respective B∗ model, then decreases as the radial distance increases, maintaining nearly the same profile but with a shift in magnitude. By day 16, it is observed that all models 39 4. Results Figure 4.19: Evolution of the total magnetic field strength and plasma beta in the disk for the B∗ series. The values are computed as the mean over a wedge around the midplane. The left panels show the total magnetic field strength, while the right panels display plasma beta (β) as a function of radius at different snapshots. The black dotted line is where β = 1 40 4. Results almost converge in the outer region with very little difference. However, as time progresses, the profiles evolve differently for each model, exhibiting changes in both shape and magnitude. A similar trend is observed in plasma beta (β), where initial similarities increase variation over time. 4.5.1.4 Density profile A density map with overlaid poloidal velocity streamlines is generated to understand the dynamics of gas flow and accretion within the disk. This figure highlights the bulk mass distribution and the direction of motion. The result is presented in Figure 4.20. For the models with lower stellar magnetic field strengths (B∗ = 50, 500, 1000 G), the highest density remains concentrated in the equatorial plane. The poloidal velocity streamlines predominantly point towards the protostar, signifying ongoing accretion, although some deviations are present, likely caused by local turbulence. Additionally, failed winds are evident across all these models, where initially ejected gas is redirected by gravity and then follows the cavity walls, eventually falling back toward the protostar. This effect is particularly clear in the B∗ = 50 G model, where there is no sign of disk truncation. But the disk is still accreting some gas at higher latitudes. For the B∗ = 2000 G model, a ejected wind is observed that pushed more gas away from the higher latitudes, which seems to have an effect, as the mass gets more concentrated towards the equatorial plane if one compares with B∗ = 500 G. The B∗ = 4000 G model exhibits a different density structure. While for R > 1 AU, most of the mass still resides in the equatorial plane, the gas becomes more dispersed at R < 1 AU. Clear signs of disk truncation appear in the inner region, forcing the gas to follow the reconnection magnetic field lines onto the stellar surface rather than accreting through the equatorial plane. This transition in accretion behavior aligns with theoretical expectations for magnetospheric truncation at higher field strengths. 41 4. Results Figure 4.20: This figure presents a color-coded density map for various B∗ models, taken at a snapshot on day 259. The black streamlines indicate the direction of poloidal velocity, vpol. 4.5.2 Surface Density Variation (Σo) Series This subsection presents the Σo series result where temperatures, density, entropy, and pressure profiles were obtained by averaging values within a wedge of 16 cells around the midplane. The radial domain was divided into 32 bins, and the mean values were computed within each bin to create a representative profile. The results are presented in Figure 4.21. 42 4. Results Figure 4.21: Temperature, total pressure, density, and entropy profiles for the Σo series, averaged over a midplane wedge and binned radially into 32 bins. The solid lines represent the profiles at t = 259 days, while the dashed lines correspond to t = 16 days. Different colors indicate models with varying initial surface density Σo = 200, 500, 1000, and 2000 g cm−2 Figure 4.21 shows the temperature, total pressure, density, and entropy profiles for the Σo series. The temperature converges in the outer regions, independent of the initial magnetic strength Σo or snapshot used. Pressure, density, and entropy also show convergence but with differences in magnitude. The Σo = 200 g cm−2 model has the lowest density in the disk, while Σo = 2000 g cm−2 has the highest through both disk and time. The same trend is seen in pressure and temperature. A noticeable dip in temperature appears around 1 AU for Σo = 200 g cm−2, where the Σo = 1000 g cm−2 model instead shows a bump. A similar feature is seen in the pressure profile at the same location. The entropy profile shows a dip at intermediate radii for all models, but the depth varies. The Σo = 200 g cm−2 model has higher entropy in the inner region than the other cases. At larger radii, entropy increases for all models. The profiles at t = 259 days follow the same shape as at t = 16 days but with lower values, showing a gradual change over time. 4.5.2.1 Radial and Azimuthal Flow Variations The analysis of the velocities for the Σo series will be presented in the same way the the B∗ series was. The first figure 4.22, illustrates velocity profiles as function of radius components, vr, vθ, and vϕ, where the Keplerian velocity is subtracted from the azimuthal component. All three velocity components for each model show the strongest variation in the inner region. As the radius increases, both the variation and magnitudes decrease. The most pronounced variations appear in the two lowest Σo models, Σo = 200 and 500 g cm−2, where both vθ and vϕ reach velocities close to 1000 km/s before rapidly 43 4. Results Figure 4.22: The figure illustrates the radial velocity profile for two snapshots, day 16 and the last snapshot day259 for the Σo series. dropping and fluctuating between 0.1 − 10 km/s. For the higher Σo models, these sharp peaks in the innermost regions disappear, and all three velocity components remain within the 0.1 − 10 km/s range for both snapshots. The vθ component behaves differently between the two snapshots. At the earlier time, it decreases rapidly with radius without reaching a high peak, whereas at the later time, it also drops but is shifted further out, showing no clear convergence. vϕ is the only one consistently converging to the same magnitude across all models and times. The figure 4.23 presents the results of the analysis of the variation of poloidal velocity, (vr + vθ), as a function of radius. 44 4. Results Figure 4.23: Variation of poloidal velocity as a function of radius in the disk for different initial surface densities Σ∗. The left panel shows the mean poloidal velocity, while the right panel displays the standard deviation of the same velocity. The final snapshot of the simulation at t = 259 days was used. The overall evolution of all models follows a similar trend, although they differ in magnitude. All four models converge with minimal differences in the outer region for the poloidal velocity’s mean and standard deviation. In the middle region of the disk, the profiles exhibit the same trend, showing decreasing magnitudes as the ra- dius increases. The most notable differences between the models appear in the inner region, where both the Σo = 200 and Σo = 500 g cm−2 models have significantly higher poloidal velocity magnitudes compared to the other two. Specifically, the Σo = 200 model is approximately one order of magnitude higher than the Σo = 1000 and Σo = 2000 g cm−2 models, while the Σo = 500 g cm−2 model reaches a peak that is up to two orders of magnitude higher. Additionally, for the Σo = 200 g cm−2 model, the initial velocity drop occurs at a larger radius compared to the Σo = 500 g cm−2 model. The right-hand panel also shows that the standard deviation peaks at the same radial distance where the mean velocity drops. In general, the standard deviation decreases with increasing radius and as the initial surface density increases. The turbulent motion of the disk is analyzed by studying the total velocity dis- persion and mean velocity at different snapshots. Additionally, the dispersion of the vϕ − vkep component, normalized by the speed of sound, is plotted in figure 4.16. 45 4. Results Figure 4.24: Velocity dispersion and mean velocity for the Σo series at different snapshots, showing the evolution of turbulent motion in the disk. It is seen in the left panels of figure 4.24 that both time and Σo influence the evolution of the total velocity dispersion and mean velocity. A thinner disk generally leads to a higher mean velocity and stronger fluctuations, which become more pronounced over time. This effect is most noticeable in the inner region, where the two lower-density 46 4. Results models exhibit higher magnitudes. Additionally, in the lowest-density model, these higher magnitudes extend further out into the disk. Although all models appear to somewhat converge in the middle region, the Σo = 200 g cm−2 model maintains a slightly higher total velocity dispersion and mean velocity. It is also noticeable that as time progresses, the convergence point for all models shifts outward in the disk. In the right panels of Figure 4.16, the models differences are clearer and increase over time. At its peak (t = 234 days), the Σo = 200 g cm−2 model reaches about 100 × cs, while the other models at the same radius remain around 0.1 × cs. Even though the models have large differences in the inner disk, they eventually converge in the outer region to the same value. Similar to the left panels, this convergence point shifts further outward as time progresses. The ratio between the radial and free-fall velocity is plotted in figure4.25. From the figure, the largest fluctuations occur for the Σo = 200 g cm−2 model, where two peaks are seen at day 184 and 259 in the inner disk. These fluctuations are sig- nificantly larger than those at earlier times. As Σo increases, the overall magnitude of the variations decreases for all models and time. The Σo = 500 model still ex- hibits noticeable fluctuations but is less extreme and more evenly distributed across time. For higher Σo models (1000 and 2000 g cm−2), fluctuations are significantly reduced, and the profiles appear more stable. 47 4. Results Figure 4.25: The figure shows the radial velocity normalized by the free-fall ve- locity, vr/vff , as a function of radius for different surface densities Σo and time. 4.5.2.2 Thermal, Magnetic, and Kinetic Pressure The pressure components are plotted for each Σo model at two different snapshots, resulting in Figure 4.26. The Σo directly impacts the pressures seen in Figure 4.26. As Σo increases, all three pressure components increase in magnitude. Despite these differences, the overall shape of each pressure component remains similar across all models. Both Ptherm and Pkin converge to their respective values at both snapshots in the disk’s outer 48 4. Results Figure 4.26: Variation of the thermal, magnetic, and kinetic pressure components with radius for two snapshots: day 16 and the final snapshot at day 259 for the Σo series. region. The models experience greater variation as time progresses. Still, the most significant difference is seen in Ptherm for Σo = 200 g cm−2, where the differenc