ABSTRACT

Title of Dissertation: MULTISCALE RADIATION-MHD SIMULATIONS
OF COMPACT STAR CLUSTERS

ChongChong He
Doctor of Philosophy, 2023

Dissertation Directed by: Professor Massimo Ricotti
Department of Astronomy

Star formation is a crucial process that lies at the center of many important topics in as-

trophysics: the nature of the first sources of radiation, the formation and evolution of galaxies,

the synthesis of elements, and the formation of planets and life. Recent advances in computing

technology have brought about unprecedented opportunities to deepen our understanding of this

complex process. In this dissertation, I investigate the physics of star formation in galaxies and

its role in shaping the galaxies and the Universe through numerical simulations.

My exploration of star formation begins with a large set of simulations of star cluster

formation from isolated turbulent Giant Molecular Clouds (GMCs) with stellar feedback using

RAMSES-RT, a state-of-the-art radiation-magneto-hydrodynamic (radiation-MHD) code. While

resolving the formation of individual stars, I have pushed the parameters (mass and density)

of the simulated GMCs well beyond the limit explored in the literature. I establish physically

motivated scaling relationships for the timescale and efficiency of star formation regulated by

photoionization feedback. I show that this type of stellar feedback is efficient at dispersing dense



molecular clouds before the onset of supernova explosions. I show that star formation in GMCs

can be understood as a purely stochastic process, where instantaneous star formation follows a

universal mass probability distribution, providing a definitive answer to the open question of the

chronological order of low- and high-mass star formation. In a companion project, I publish the

first study of the escape of ionizing photons from resolved stars in molecular clouds into the in-

tercloud gas. I conclude that the sources of photons responsible for the epoch of reionization,

one of the most important yet poorly understood stages in cosmic evolution, must have been very

compact star clusters, or globular cluster progenitors, forming in dense environments different

from today’s galaxies.

In follow-up work, I use a novel zoom-in adaptive-mesh-refinement method to simulate the

formation and fragmentation of prestellar cores and resolve from GMC scales to circumstellar

disk scales, achieving an unprecedented dynamic range of 18 orders of magnitude in volume

in a set of radiation-MHD simulations. I show that massive stars form from the filamentary

collapse of dense cores and grow to several times the core mass due to accretion from larger

scales via circumstellar disks. This suggests a competitive accretion scenario of high-mass star

formation, a problem that is not well understood. We find that large Keplerian disks can form

in magnetically critical cores, suggesting that magnetic braking fails to prevent the formation

of rotationally-supported disks, even in cores with mass-to-flux ratios close to critical. This is

because the magnetic field is extremely turbulent and incoherent, reducing the effect of magnetic

braking by roughly one order of magnitude compared to the perfectly aligned and coherent case,

which proposes a solution to the “magnetic braking catastrophe.”



MULTISCALE RADIATION-MHD SIMULATIONS
OF COMPACT STAR CLUSTERS

by

ChongChong He

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2023

Advisory Committee:
Professor Massimo Ricotti, Chair/Advisor
Professor Cole Miller
Professor Alberto Bolatto
Professor Carter Hall, Dean’s Representative
Professor Susan Kassin, External Examiner



© Copyright by
ChongChong He

2023



Preface

The research presented in Chapters 2-5 and Appendix A of this dissertation has either

been published or submitted for publication in the Monthly Notices of the Royal Astronomical

Society (MNRAS) or the Astrophysical Journal (ApJ) with minimal modification. Each chapter

corresponds to the following papers:

• Chapter 2 – C.-C. He, M. Ricotti, and S. Geen, “Simulating star clusters across cosmic time

- I. Initial mass function, star formation rates, and efficiencies”, MNRAS 489, 1880–1898

(2019).

• Chapter 3 – C.-C. He, M. Ricotti, and S. Geen, “Simulating star clusters across cosmic time

- II. Escape fraction of ionizing photons from molecular clouds”, MNRAS 492, 4858–4873

(2020).

• Chapter 4 – C.-C. He and M. Ricotti, “Massive prestellar cores in radiation-magneto-

turbulent simulations of molecular clouds”, MNRAS, 10.1093/mnras/stad1289 (2023)

• Chapter 5 – C.-C. He and M. Ricotti, “Magnetic braking fails to work: formation of large

Keplerian disks in magnetically critical giant molecular clouds”, submitted to MNRAS,

(2023).

• Chapter 6 – C.-C. He, “A fast and accurate analytic method of calculating galaxy two-point

correlation functions”, ApJ 921, 59, 59 (2021).

ii

https://ui.adsabs.harvard.edu/abs/2019MNRAS.489.1880H/abstract
https://ui.adsabs.harvard.edu/abs/2020MNRAS.tmp..162H/abstract
https://doi.org/10.1093/mnras/stad1289
https://ui.adsabs.harvard.edu/abs/2021ApJ...921...59H/abstract


The contents of this dissertation form a coherent story. In the first act, I conduct a large

set of simulations to explore the laws of star formation. In the second act, I investigate the

implications of star cluster formation for the source of radiation during the epoch of reionization.

In the third and fourth acts, I delve into smaller scales, focusing on the collapse of prestellar cores

and the formation of circumstellar disks. Finally, as an appendix, I present part of my ongoing

project on stellar dynamics, which is the subsequent evolution of those simulated star clusters.

Throughout, we witness how each of the many characters of star formation drives the behavior

of the natal cloud, which ultimately affects star formation itself.

Originally, I researched star formation to study the source of radiation for cosmic reioniza-

tion but ended up discovering some mysterious and intriguing theories of star formation, which

lured me to dive deeper into the subject. In the end, I attempted to uncover its identity and push

back its frontier, but it appears what I discovered so far is just the tip of the iceberg.

iii



Dedication

To my parents.

iv



Acknowledgments

I am thrilled to present this dissertation as the culmination of my PhD journey in astro-

physics at the University of Maryland. This work would not have been possible without the

guidance and support of many people who deserve my sincere gratitude.

First and foremost, I want to thank my advisor, Professor Massimo Ricotti, for his men-

torship, encouragement, and inspiration. He has taught me how to be a rigorous and creative

researcher, how to communicate my ideas effectively, and how to enjoy the process of discovery.

He has always been generous with his time and resources and has given me the freedom and

confidence to pursue my own interests. I am honored to have worked with him and learned from

him.

I also want to thank my thesis committee members, Cole Miller, Alberto Bolatto, Carter

Hall, and Susan Kassin; my collaborator, Sam Geen; and my research group, Benedikt Diemer,

Jongwon Park, Calvin Osinga, Katya Leidig, Ronan Hix, and Fred Garcia. They have provided

valuable feedback, insights, and suggestions that have improved the quality and clarity of my

work. I am also grateful to my first research advisor, Laurens Keek, for introducing me to the

field of astrophysics as an exchange student at Georgia Tech. He sparked my curiosity and passion

for this subject and encouraged and supported me to pursue a Ph.D. I owe him a debt of gratitude.

I would like to thank all the UMD Astronomy Department Staff – Eric McKenzie, Mary

Ann Phillips, Barbara Hansborough, John Cullinan, Olivia Dent, Natalie Rowe, Susan Lehr,

v



Dorinda Kimbrell, Mark Wolfire, Kevin Rauch, and Peter Teuben – for making the administrative

end of things so easy for us. They have always been helpful, efficient, and cheerful in dealing

with any issues or requests that I had. I also owe a huge debt of gratitude to my wonderful office

mates, Sara Frederick, Robyn Smith, and Pradip Gatkine, for their support during my first years

in the department. I would also like to thank my classmates, Weizhe, Jialu, Laura, Carrie, and Liz

for keeping me company through this journey and helping me adjust to a new country. Finally, I

want to thank all of the graduate students – Rebecca, Ginny, Kyle, Amy, Drew, Ramsey, Milena,

Charlotte, Guangwei, Harrison, Julian, Vicente, Sergio, Tomas, Rye, Alex, Ben, Erica, Teal,

Krista, Dana, Thomas, Mark, Mahmuda, Myra, Zeeve, Vicki, Ell, Joe, Jordan, Arjun, Gokul,

and more – and everyone else in the Department of Astronomy. They have been a source of

friendship, fun, and support throughout my Ph.D. They have also been friendly and supportive

colleagues and created a warm and welcoming atmosphere in the department, which made my

research experience enjoyable and rewarding. I could not have asked for a better place or better

people to spend the past years with.

It’s hard to express how grateful I am for my parents, He Xinyuan and Li Chulan. I would

like to extend my heartfelt thanks to them for their unconditional love and support. They have

always believed in me and sacrificed a lot for me, giving me everything I needed to become who

I am. They have been my biggest fans and cheerleaders. I dedicate this dissertation to them.

This dissertation marks the end of one chapter of my life and the beginning of another. I

am excited about what lies ahead, but I will always cherish the memories and lessons from this

amazing journey.

vi



Table of Contents

Preface ii

Dedication iv

Acknowledgements v

Table of Contents vii

List of Tables x

List of Figures xi

List of Abbreviations xiv

Chapter 1: Introduction 1
1.1 Star Formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.1.1 The star formation rate . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Initial Mass Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.3 Magnetic field and the magnetic braking problem . . . . . . . . . . . . . 7

1.2 Epoch of Reionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Techniques for Simulating Star Formation . . . . . . . . . . . . . . . . . . . . . 13
1.4 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.5 A summary of software and facilities . . . . . . . . . . . . . . . . . . . . . . . . 17

Chapter 2: Simulating Star Clusters Across Cosmic Time: I. Initial Mass Function, Star
Formation Rates and Efficiencies 18

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.1 The IMF and SFE of Molecular Clouds . . . . . . . . . . . . . . . . . . 22

2.2 Numerical Simulations and Methods . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.1 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2.2 Resolution and Sink Formation . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.3 Feedback and Properties of UV source . . . . . . . . . . . . . . . . . . . 32
2.2.4 Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.3.1 Stellar Initial Mass Function . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.2 Star Formation Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.3.3 Star formation law in molecular clouds . . . . . . . . . . . . . . . . . . . 54

vii



2.3.4 Effects of Lowering the Gas Metallicity . . . . . . . . . . . . . . . . . . 58
2.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Chapter 3: Simulating Star Clusters Across Cosmic Time: II. Escape Fraction of Ion-
izing Photons from Molecular Clouds 65

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.2 Numerical Simulations and Methods . . . . . . . . . . . . . . . . . . . . . . . . 70

3.2.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.2.2 Calculation of the Ionizing Escape Fraction . . . . . . . . . . . . . . . . 74

3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.3.1 Sky maps of the Escaping Ionizing Radiation . . . . . . . . . . . . . . . 79
3.3.2 Time Evolution of the Sky-Averaged Escape Fraction . . . . . . . . . . . 81
3.3.3 Time-Averaged Escape Fraction ⟨fMC

esc ⟩ . . . . . . . . . . . . . . . . . . . 84
3.3.4 Escape Fraction of Helium Ionising Photons . . . . . . . . . . . . . . . . 92
3.3.5 Absorption by Dust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
3.4.1 Analytic Modelling and Interpretation of ⟨fMC

esc ⟩ . . . . . . . . . . . . . . 96
3.4.2 Ionising Photons from OB Associations . . . . . . . . . . . . . . . . . . 102

3.5 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

Chapter 4: Massive Prestellar Cores in Radiation-magneto-turbulent Simulations of Molec-
ular Clouds 108

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4.2 Methods and Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.3 Results. I. Turbulent massive disks . . . . . . . . . . . . . . . . . . . . . . . . . 118

4.3.1 Evolution of a 27 solar mass core (Core A-hr) . . . . . . . . . . . . . . . 123
4.3.2 Evolution of a 130 solar mass core . . . . . . . . . . . . . . . . . . . . . 131
4.3.3 The disk thickness is determined by magnetic support and supersonic tur-

bulent motions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.4 Results. II. Low-mass stars form from the fragmentation of massive pre-stellar

cores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.4.1 Fragmentation into low-mass stars before the formation of a steady disk

structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.4.2 Star formation efficiency in cores and multiplicity . . . . . . . . . . . . . 141

4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
4.5.1 Formation of ultra-high-mass stars – a competitive accretion scenario . . 145
4.5.2 UV radiation trapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
4.5.3 Influence of metallicity on disk stability . . . . . . . . . . . . . . . . . . 151

4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

Chapter 5: Magnetic Braking Fails to Work: Formation of Large Keplerian Disks in
Magnetically Critical Giant Molecular Clouds 155

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
5.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163

viii



5.3.1 Morphology of Gas Density and Magnetic Field in Cores . . . . . . . . . 163
5.3.2 B − ρ Relationship . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
5.3.3 Magnetic and Turbulent Support in the Cores and Disks . . . . . . . . . 174
5.3.4 Magnetic Braking Problem . . . . . . . . . . . . . . . . . . . . . . . . . 175
5.3.5 Evidence for Magnetically-Driven Winds . . . . . . . . . . . . . . . . . 180

5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182

Chapter 6: A Fast and Accurate Analytic Method of Calculating Galaxy Two-point
Correlation Functions 184

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
6.2 Analytic Formulas of

⌢

RR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
6.2.1 Rectangles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
6.2.2 Cuboids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192
6.2.3 Circles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
6.2.4 Spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198

6.3 Accounting for Edge Corrections in
⌢

RR . . . . . . . . . . . . . . . . . . . . . . 200

6.4 Computing
⌢

DR Analytically in O(Ng) Time . . . . . . . . . . . . . . . . . . . . 201
6.4.1 Rectangles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
6.4.2 Cuboids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
6.4.3 Circles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
6.4.4 Spheres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205

6.5 Comparisons with the Monte Carlo Method . . . . . . . . . . . . . . . . . . . . 206
6.6 Summary and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209

Chapter 7: Conclusion and Future Work 211
7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214

7.2.1 Resolving the Complete IMF of a Star Cluster . . . . . . . . . . . . . . . 214
7.2.2 Dynamics of Massive Compact Star Clusters and SMBH Growth . . . . . 215
7.2.3 GPU-accelerated Computing Methods for Astrophysics in the Era of Ex-

ascale Computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217

Chapter A: Appendices for Chapter 2 219
A.1 Clump finder criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
A.2 Emission from clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220

Chapter B: Appendices for Chapter 3 222
B.1 Converting column density to escape fraction . . . . . . . . . . . . . . . . . . . 222

Chapter C: Appendices for Chapter 5 223
C.1 Mass-to-flux ratio of a non-singular isothermal sphere . . . . . . . . . . . . . . . 223
C.2 Resolution study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224

Chapter D: Appendices for Chapter 6 225

D.1 Algorithms to Compute
⌢

DR Analytically in O(Ng) Time . . . . . . . . . . . . . 225

ix



List of Tables

2.1 Initial conditions of our 16 simulations. . . . . . . . . . . . . . . . . . . . . . . 29
2.2 A collection of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

3.1 A table of parameters in all simulations. . . . . . . . . . . . . . . . . . . . . . . 71
3.2 A summary of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.3 Escape fraction at the Lyman edge with and without dust extinction. . . . . . . . 95

4.1 Summary of the properties of the simulated cores, protostellar disks, and forming
stars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

5.1 List of zoom-in simulations presented in this chapter . . . . . . . . . . . . . . . 159

6.1 Percentage errors ϵ(r̂) of the calculated
⌢

RR (r̂) when certain edge corrections are
not included . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200

x



List of Figures

1.1 The IMF measured in a radiation-hydrodynamic simulation of the collapse of a
molecular cloud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1 Simulation parameters in this work . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2 Line-of-sight projections of the gas density for three simulations . . . . . . . . . 36
2.3 Same as Fig. 2.2 but showing the density-weighted projection of the temperature. 37
2.4 Temporal evolution of the mass functions of the stars . . . . . . . . . . . . . . . 38
2.5 Comparing IMFs from simulations with different resolutions . . . . . . . . . . . 41
2.6 Numerical experiments of fragmenting each sink particle into stars . . . . . . . . 42
2.7 The IMF of the simulated star clusters along with the best fit power-law . . . . . 44
2.8 Maximum stellar mass in a cluster v.s. the mass of the star cluster . . . . . . . . . 46
2.9 Dimensionless star formation efficiency as a function of the dimensionless time . 48
2.10 Stellar mass and total SFE of the clusters as a function of the initial mass of the

gas cloud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.11 Dimensionless star formation efficiency and dimensionless star formation rate per

free-fall time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.12 Maximum dimensionless star formation rate per free-fall time and the ratio of

star-formation time to sound-crossing time . . . . . . . . . . . . . . . . . . . . . 55
2.13 Same as Fig. 2.7 but for the L-C cloud with various metallicities. . . . . . . . . . 58
2.14 Same as Fig. 2.9 but for the L-C cloud with different gas metallicities, as shown

in the legend. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.15 Comparing gas temperature from simulations with different metallicities . . . . . 60

3.1 Time sequence plot of line-of-sight projections of density-weighted gas density . 73
3.2 Angular distribution of escape fraction of ionizing photons . . . . . . . . . . . . 77
3.3 Time evolution of the LyC emission rate (Q), escaping rate (Qesc), and escaping

fraction (fesc ≡ Qesc/Q) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.4 Time evolution of the SFE, hydrogen ionizing-photon emission rate, and escaping

rate at various metallicities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.5 The total escape fraction of ionizing photons . . . . . . . . . . . . . . . . . . . . 85
3.6 ⟨fMC

esc ⟩ v.s. SFE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.7 Hydrogen-ionizing photons per unit star cluster mass . . . . . . . . . . . . . . . 87
3.8 Mean escaped ionizing-photon emission rate and duration . . . . . . . . . . . . . 89
3.9 Fractional cumulative radiation emission and escaping . . . . . . . . . . . . . . . 90
3.10 Escape fraction of photons as a function of hν . . . . . . . . . . . . . . . . . . . 94
3.11 Ratio of tMS(Mmax) to the measured tuv . . . . . . . . . . . . . . . . . . . . . . 99
3.12 Comparing model ⟨fMC

esc ⟩ with ⟨fMC
esc ⟩ from simulations . . . . . . . . . . . . . . . 99

xi



3.13 Conversion from theR parameter to ⟨fMC
esc ⟩, following Eq. (3.15). . . . . . . . . . 100

4.1 Outline of the simulation method . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.2 Gallery of the collapsing prestellar cores A-hr, B, C, and D . . . . . . . . . . . . 117
4.3 Schematic of the two modes of prestellar core fragmentation and disk formation . 119
4.4 Density slices of the ∼ 27 M⊙ core (Core A-hr) at various ages . . . . . . . . . . 124
4.5 Radial and vertical profiles of the 27 M⊙ core A-hr . . . . . . . . . . . . . . . . 126
4.6 Same as Figure 4.4 but for the very massive core (Core B) . . . . . . . . . . . . . 130
4.7 Same as Figure 4.5 but for the ∼ 130 M⊙ core B . . . . . . . . . . . . . . . . . . 132
4.8 Thermal, magnetic, and turbulent properties of the cores and disks . . . . . . . . 135
4.9 Snapshots of projected density of core A-hr . . . . . . . . . . . . . . . . . . . . 139
4.10 The mass distribution of the stars forming in each core in our ‘zoom-in’ simula-

tions compared to the stars that form from the same core in the baseline run . . . 140
4.11 The growth of Core B and Core A-hr, demonstrating a “turbulent core” scenario

and a turbulent core scenario . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
4.12 The trapping and escaping of H II regions . . . . . . . . . . . . . . . . . . . . . 149
4.13 The escaping of UV radiation from a dense filament due to dynamical motion . . 149

5.1 Visualization of the cores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
5.2 Time evolution of the four cores . . . . . . . . . . . . . . . . . . . . . . . . . . 164
5.3 Time evolution of the magnetic field morphology . . . . . . . . . . . . . . . . . 165
5.4 3D rendering of the magnetic field lines . . . . . . . . . . . . . . . . . . . . . . 165
5.5 The global B − ρ relationship at GMC scale . . . . . . . . . . . . . . . . . . . . 167
5.6 Enhancement of the magnetic intensity at high density . . . . . . . . . . . . . . . 168
5.7 The magnetic, turbulent, and thermal support of the cores and disks . . . . . . . . 173
5.8 Explanation of the weak overall magnetic braking . . . . . . . . . . . . . . . . . 176
5.9 The extremely turbulent and incoherent magnetic field on the circumstellar disks

as a solution to the magnetic braking problem . . . . . . . . . . . . . . . . . . . 179
5.10 Demonstration of magnetically driven outflow from the disk . . . . . . . . . . . 181

6.1 Diagrams showing edge corrections in the calculation of random-random and
data-random pair counts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189

6.2 Comparing the two-point correlation functions calculated from the analytic method
from this work with that from brute-force Monte Carlo method . . . . . . . . . . 207

6.3 Comparing the speed of the analytic approach from this work to that of the tradi-
tional Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208

7.1 Demonstration of the set of proposed simulations in the context of my past work
and the literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215

A.1 Explanation of the sink formation criteria . . . . . . . . . . . . . . . . . . . . . 219
A.2 Ionising photon emission rate as a function of stellar mass . . . . . . . . . . . . . 220
A.3 He0 and He+ ionising photon emission rate as a function of the star cluster mass . 221
B.1 Color plots of f(τ0, T ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222

xii



C.1 Comparing two definitions of the relative importance of the gravitational and
magnetic forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224

C.2 A comparison of disk structure between two simulations with different resolutions 225

D.1 An algorithm for precise calculations of
⌢

DR (r) in O(Ng) time, applying to
rectangular regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226

D.2 An algorithm for precise calculations of
⌢

DR (r) in O(Ng) time, applying to
cuboidal regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 227

D.3 An algorithm for precise calculations of
⌢

DR (r) in O(Ng) time, applying to
circular regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228

D.4 An algorithm for precise calculations of
⌢

DR (r) in O(Ng) time, applying to
spherical regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228

xiii



List of Abbreviations

AGN Active galactic nucleus
ALMA Atacama Large Millimeter/submillimeter Array
AMR Adaptive mesh refinement
AU Astronomical unit
ApJ Astrophysical Journal
BH Black hole
CA Competitive accretion
CMF Core mass function
CO Carbon monoxide
CPU Central processing unit
DD Data-data
DR Data-random
EAGLE Evolution and Assembly of GaLaxies and their Environments
EoR Epoch of reionization
GC Globular cluster
GMC Giant molecular cloud
GPU Graphical processing unit
H Hydrogen
H2 Molecular hydrogen
HM high-mass
HPC High Performance Computing
HST Hubble Space Telescope
He Helium
He II Double-ionized helium
IGM Intergalactic medium
IM Intermediate-mass
IMBH Intermediate-mass black hole
IMF Initial Mass Function
IR Infrared
ISM Interstellar medium
JWST James Webb Space Telescope
KH Kelvin-Helmholtz
LyC Lyman-continuum
MC Molecular cloud

xiv



MHD Magneto-hydrodynamics
MNRAS Monthly Notices of the Royal Astronomical Society
NSC Nuclear star cluster
PDF Probability density function
QSO Quasi-stellar object
RR Random-random
RT Radiation transfer
SFE Star Formation Efficiency
SFR Star formation rate
SMBH Supermassive black hole
SN Supernova
SPH Smoothed-particle hydrodynamics
TC Turbulent Core
TPCFs Two-point correlation functions
TSFE Total star formation efficiency
UV Ultraviolet
YSO Young stellar objects
H I Neutral hydrogen
H II Singly-ionized hydrogen

Simulation names used in Chapter 2:

XXS Extra-extra-small
XS Extra-small
S Small
M Medium
L Large
XL Extra-large
F Fiducial
C Compact
VC Very compact

xv



Chapter 1: Introduction

Stars are the fundamental building blocks of the universe. The problem of how stars form

lies at the center of many important topics in astrophysics: the nature of the first sources of

radiation, the formation and evolution of galaxies, the synthesis of elements, and the formation

of planets and life.

Stars are not born in isolation; they are mostly found in groups called star clusters. These

clusters form from molecular clouds, which serve as nurseries for stellar objects. They are sites

of intense stellar feedback processes such as UV feedback, jets, winds, and supernovae.

It is difficult to formulate a general theory for star formation due to the wide range of

physical processes involved. However, recent advances in computing technology have brought

about unprecedented levels of understanding of this complex process. My research has focused

on understanding the physics of star formation in galaxies and its role in shaping the galaxies

and the universe, using multiscale radiation magneto-hydrodynamics (MHD) simulations of star

formation from molecular clouds.

This chapter begins by reviewing two major problems in star formation: the star formation

rate and the initial mass function, which provide background for Chapter 2 and Chapter 4. I

then review the role of magnetic field in star formation, raising an open question that motivates

the work presented in Chapter 5. Next, I provide a brief introduction to the topic of cosmic

1



reionization, which serves as a background for Chapter 3. Then, I present a summary of the

techniques used to simulate star formation. Finally, this chapter closes with the thesis outline and

a list of software and facilities used in this thesis.

The essence of this thesis and what we hope to convey is that we desperately need more

physically-grounded models of star formation to advance our understanding of this fundamental

process in the universe. We are in the midst of an exciting period of rapid advancement in

computational technology, and future work with the help of large-scale simulations and advanced

observations holds the promise of uncovering the full mystery of star formation.

1.1 Star Formation

1.1.1 The star formation rate

One parameter to characterize how fast stars form in a region is the dimensionless star

formation rate (SFR) per free-fall time,

ϵff =
ṁ∗
mgas

tff , (1.1)

where tff is the free-fall time of the gas at average density and ṁ∗ is the SFR. A gas cloud that

is not supported against collapse will naturally produce stars at a rate corresponding to ϵff ≈ 1,

since the free-fall time is the natural evolutionary timescale for a self-gravitating system with

nothing inhibiting star formation. A direct way to measure ϵff is to estimate the SFR by counting

young stellar objects (YSOs) and estimating their masses and the duration of the YSO phase

(e.g. Krumholz et al., 2012). A second approach is to match catalogs of star-forming regions

2



identified by tracers such as IR or free-free emission with catalogs of molecular clouds identified

by CO or dust emission, and then using the mass, free-fall time, and SFR of the matched clouds

and star-forming regions to estimate ϵff (e.g. Lee et al., 2016). A third approach, available for

extragalactic systems with extensive molecular gas and star-formation tracer maps, is to pixelate

the entire galaxy and estimate masses, densities, and free-fall times in each pixel (e.g. Krumholz

et al., 2012).

These methods, applicable to scales from the Milky Way and its satellites to extragalactic

systems, provide reasonable consistency in results. In a recent review article, Krumholz et al.

(2019) conclude that the preponderance of the current observational evidence favors ϵff ≈ 0.01

for regions ≳ 1 pc in size with small dispersion and systematic uncertainty.

The observed value of ϵff is surprisingly low. A number of authors have proposed theoret-

ical models aiming at explaining this low ϵff . Earlier models propose magnetic regulation as the

reason. Recent measurements of magnetic fields have shown that they are too weak to support

the gas (Crutcher, 2012). The unbound cloud models propose that ϵff is low because most of

the materials observationally defined as molecular clouds are not actually self-gravitating. How-

ever, this hypothesis does not explain the YSO counting studies where the gas is almost certainly

bound (Kauffmann et al., 2013). Consequently, attention has been focused on two possibilities:

turbulence support and feedback regulation.

There is considerable evidence that GMCs obey the relations discovered by Larson (1981):

the 1D velocity dispersion, σ, is supersonic and varies with the size as σ ∝ Lp, where p ≈ 0.5 in

the Milky Way. At a randomly chosen subregion of a cloud, the binding energy per unit mass of

a region of mass M scales as GM/R ∝ R2, assuming a mean density, while its kinetic energy

per unit mass scales as σ2 ∝ R, therefore the virial ratio obeys αvir ∝ 1/R. A virial parameter,

3



or virial ratio, is defined as the ratio of the total kinetic energy to the gravitational energy of a

dynamical system, αvir ≡ 2K/|W|, where the numerical coefficient is chosen so that αvir = 1

indicates virial equilibrium (W = −2K), and αvir > 1 (αvir < 1) means unbound (bound)1. Most

regions smaller than the size of the cloud are unbound (Krumholz & McKee, 2005). Numerous

recent works have focused on modeling the density PDF and the evolution of the self-gravitating

parts of a cloud with increasing accuracy (e.g. Hennebelle & Chabrier, 2011). The consensus of

these models is that turbulence does substantially reduce ϵff but not to ϵff ≈ 0.01 as required by

observations.

The density PDF of a self-gravitating cloud deviates and develops a prominent power-law

tail on its high-density end that causes ϵff to rise with time. This density build-up is likely to

be counterbalanced by feedback processes that break up high-density clumps. Recently, Grudić

et al. (2018) propose that cloud disruption by stellar feedback is so fast and efficient that a typical

cloud never has time to develop power-law tails substantial enough to bring ϵff to large values.

However, the nature and existence of turbulence in molecular clouds are still not well under-

stood. There is considerable debate in the literature between the turbulent theory and the theory

of coherent gravitational collapse (Heitsch et al., 2009). Although there is significant evidence

from analytic theories (Klessen & Hennebelle, 2010), numerical simulations (e.g. Robertson &

Goldreich, 2012), and kinematic measurements with Gaia that accretion flows inevitably drive

turbulence, the question of the low value of ϵff is still unsettled.

1Note that there is another definition commonly used in the literature where αvir ≡ K/|W| and αvir = 0.5
indicates virial equilibrium.

4



1.1.2 Initial Mass Function

Zooming in even further from star-forming molecular clouds, we reach the scale of indi-

vidual stars. One of the properties of stars that is important for determining their observable

characteristics and evolutionary path is their mass. The distribution of stellar masses at birth is

known as the initial mass function (IMF). Almost all inferences of light or physical properties for

unresolved stellar populations, as well as most models of galaxy formation, rely on an assumed

form of the IMF.

The first attempts to measure the IMF were carried out by Salpeter (1955), using stars in

the Solar neighborhood. Observations using the field star, young cluster, and globular cluster

methods all appear to produce roughly consistent results. Current evidence suggests that the IMF

appears to be quite universal in many locations throughout the Milky Way, with the possible

exception of star clusters formed very near the Galactic Center. The IMF has a distinct peak

in the mass range 0.1 − 1M⊙ and falls off as a powerlaw dN/dm ∝ m−α at higher masses,

with α ≈ 2.35, the value originally determined by Salpeter (1955). A widely-used functional

representation of the IMF is a two-piece powerlaw given by Kroupa (2002):

Φ(m) =
dN

d logm
∝





m−(α1−1) (m0 < m < m1)

km−(α2−1) (m > m1)

(1.2)

with m0 = 0.08 M⊙, m1 = 0.5 M⊙, α1 = 1.3, α2 = 2.3, and k = 0.5 (to guarantee continuity

across the powerlaw break). A large uncertainty is in the brown dwarf regime below 0.08 M⊙,

where there is a clear fall-off from the peak, but its exact functional form is poorly determined.

5



There may also be an upper cutoff between 100 and 150 M⊙ (Figer, 2005), although whether

this is an aspect of star formation or a result of a sharp increase in instability and mass loss is

undetermined (Tan et al., 2014).

Various theoretical models have been constructed to explain the universality of the IMF

based on gravoturbulent fragmentation of the host cloud (Padoan et al., 1997; Padoan & Nord-

lund, 2002; Mac Low & Klessen, 2004; Hennebelle & Chabrier, 2008; Hopkins, 2012a). Early

pioneering numerical simulations with pure hydrodynamics and gravity based on this model do

sometimes return a mass distribution that looks much like the empirical IMF, as illustrated in

Figure 1.1. However, this is largely depending on the choice of initial conditions. Simulations

that start with more realistic initial conditions, namely with turbulent density structures or with

supervirial/subvirial initial state, do not appear to collapse as predicted by the empirical IMF

(e.g., Clark et al., 2008). Besides, they were often limited in terms of statistics, resolution, or

cloud size.

More recent work, with increasing computing power, provided more reliable statistics and

IMF distributions (e.g., Bonnell et al., 2011; Girichidis et al., 2011; Krumholz et al., 2011; Bate,

2012; Ballesteros-Paredes et al., 2015). Furthermore, radiative stellar feedback has been invoked

to explain the precise shape of the IMF, using both simulations (e.g., Bate, 2009a; Krumholz

et al., 2011; Gavagnin et al., 2017) and analytic models (e.g., Guszejnov & Hopkins, 2016). The

core idea is that radiative feedback shuts off fragmentation at a characteristic mass scale that

sets the peak of the IMF. Simulations including radiation seem to support the idea that stellar

feedback can pick out a characteristic peak mass. However, the limitation of this hypothesis

is that it has little to say about the powerlaw tail of the IMF. Indeed, by sacrificing resolution,

Gavagnin et al. (2017) captured a large number of massive stars, which emit large quantities of

6



Figure 1.1: The IMF measured in a radiation-hydrodynamic simulation of the collapse of a
500 −M⊙ cloud with uniform initial density. The single-hashed region gives all objects, while
the double-hashed region gives those that have stopped accreting. Credit: Bate (2009b).

ionizing radiation, and argued that this alters the high mass end of the IMF. So far no radiation-

MHD simulation has managed to follow the evolution of a cloud with realistic size and resolve

the formation of individual high- and low-mass stars that match the empirical IMF.

1.1.3 Magnetic field and the magnetic braking problem

Molecular clouds are observed to be permeated by magnetic fields (Crutcher, 1999; Lee

et al., 2017). Observers use the Zeeman effect to measure magnetic fields in GMCs and find

strengths ranging from tens to thousands of µG, with higher-density gas generally showing

stronger fields. For a low-density envelope of a GMC with n ∼ 100 cm−3 (ρ ∼ 10−22g cm−3),

the typical dynamic velocity v ∼ a few km s−1, giving a kinetic energy density

eK =
1

2
ρv2 ∼ 10−11 erg cm−3. (1.3)

7



With a magnetic field strength of 20 µG, typical of molecular clouds on large scales, the energy

density is

eB =
B2

8π
∼ 10−11 erg cm−3, (1.4)

comparable to the kinetic energy density. Therefore, the magnetic field is dynamically significant

in the flow of gas.

The strong magnetic field in high-density prestellar cores can in principle strongly affect the

evolution of angular momentum during the core collapse. The twisting of the magnetic field lines

produced by disk rotation in the flux-freezing regime of ideal magnetohydrodynamics (MHD),

can apply a force counter to the rotation velocity, also known as magnetic braking, effectively

slowing down rotation and increasing radial gas infall. In idealized numerical MHD simulations,

the timescale of the braking can become so short that protostellar disks fail to form or are much

smaller than the observed sizes, a phenomenon known as the “magnetic braking catastrophe” in

the theoretical literature of disk formation (e.g. Allen et al., 2003; Galli et al., 2006; Hennebelle

& Fromang, 2008; Li et al., 2014). What is the minimum required magnetic field strength so that

magnetic braking significantly affects disk formation?

To answer this question, let us first examine how magnetic fields affect the collapse of a

molecular cloud. Magnetic fields support charged gas against gravitational collapse. A common

characterization of the relative importance of the gravitational and magnetic forces in a molecular

cloud or core is the normalized mass-to-flux ratio,

µ ≡ M/ΦB

MΦ/ΦB

=
M

MΦ

, (1.5)

8



where M is the total mass contained within a spherical region of radius R, ΦB = πR2B is the

magnetic flux threading the surface of the sphere assuming a uniform magnetic field strength B,

and

MΦ = cΦ
ΦB√
G

(1.6)

is the magnetic critical mass, the mass at which the magnetic and gravitational forces balance

each other. The constant cΦ is a dimensionless coefficient that depends on the assumed geometry

of the system. For a spherical cloud of uniform density, cΦ =
√
10/(6π) = 0.168. Note that

the definition of the normalized mass-to-flux ratio µ is simplified here because the critical value

depends on the geometry of the gas and the magnetic fields.

In a sub-critical cloud (defined as a cloud with µ < 1), the magnetic field should prevent

the collapse of the cloud core altogether. Analytical predictions (Joos et al., 2012) suggest that

there are no centrifugally-supported disks in models with µ ≤ 10, although there are disk-like

over-densities of gas in which the magnetic fields, rather than the centripetal force, support the

gas against collapse in the radial direction. Observations suggest typical values of µ ≈ 2 − 10

in molecular cloud cores (e.g. Crutcher, 1999; Bourke et al., 2001), and this value could be even

smaller after correcting for projection effects (Li et al., 2013). Theoretically, disk formation

should be completely suppressed in the strict ideal MHD limit for the level of core magnetization

deduced from observation – analytic study and numerical simulations have shown that the angular

momentum of the idealized collapsing core is nearly completely removed by magnetic braking

close to the central object (e.g., Mestel & Spitzer, 1956; Mellon & Li, 2008).

However, these results seem to be in contrast to recent high-resolution observations reveal-

ing the existence of Keplerian disks around Class 0 protostellar objects (e.g. Tobin et al., 2012;

9



Codella et al., 2014; Lee et al., 2017; Johnston et al., 2020). How the catastrophe is averted for

disk formation is still a question under debate. Recent studies have revealed that the catastrophe

can be avoided if the magnetic field and the rotation axis are not aligned (Joos et al., 2012; Gray

et al., 2018), the magnetic field is diffused by turbulent velocity field (Santos-Lima et al., 2012;

Joos et al., 2013), or the magnetic field is less coherent (Seifried et al., 2013). However, these

results largely depend on the choice of artificial initial conditions and fall short on explaining the

existence of large (> 1000 AU) disks revealed by recent radio/mm and optical/IR observations

(van Kempen et al., 2012; Takahashi et al., 2012; Johnston et al., 2015, 2020). This motivates the

work presented in Chapter 5.

1.2 Epoch of Reionization

The epoch of reionization (EoR) is a period in the early universe when the first stars and

galaxies formed and began to emit ultraviolet radiation. This radiation ionized the neutral hydro-

gen atoms that filled the universe, marking the transition from the “Dark Ages” to the “Cosmic

Dawn”. Cosmic reionization involves the coupling of galaxy formation with the physics of grav-

ity and radiation transfer to produce a global phase transition of ionization, making for a complex

problem.

Ionization occurs when photons with energies E ≥ IH = 13.6 eV, a.k.a. ionizing photons,

or Lyman-continuum (LyC) photons, interact with neutral hydrogen atoms. The radiation from

luminous objects will ionize the surrounding IGM once it escapes from its host galaxies. This

ionized IGM is known as cosmological H II regions. Recombination occurs when the Coulomb

force attracts protons and electrons, which is efficient at temperatures T ≲ 104 K. In regions with

10



ionizing radiation, the gas approaches ionization equilibrium when the recombination rate equals

the ionization rate. The recombination rate is proportional to the product of the number density

of protons and electrons, nenp.

The H II regions then evolve. The ionizing photons in a radial direction penetrate until all

of them are absorbed by newly recombined atoms inside the H II region and there is no more

flux left at the ionization front. The radius of the H II region at this initial stage is determined

by the Stromgren radius Rs, by balancing the total number of ionizations and recombinations.

At this stage, the H II region is heated to over 10, 000 K and is embraced by the cold ambient

medium. The gas expands because it has higher pressure than its surroundings, and it transitions

to the second intermediate stage after the time it takes a sound wave to cross the H II region. As

the shock wave sweeps up most of the gas in its path and accumulates mass, the H II region ex-

pands and diffuses. Eventually, the H II region comes into pressure equilibrium with the ambient

medium, and the ionization front stalls out at the final Stromgren radius. Any ionizing radiation

that escapes from the galaxy creates a cosmological H II region, which is the building block of

cosmic reionization.

What sources are responsible for producing the required ionizing radiation in the EoR? Al-

though quasi-stellar objects (QSOs, or quasars) are some of the brightest objects in the Universe,

the latest studies have shown that they only contribute 1− 5% of the UV photon budget at z = 6

due to their low number densities (e.g. Willott et al., 2010). X-rays from compact binaries can

penetrate much deeper into the IGM and create large partially ionized regions to a distance of 100

kpc with ionization fraction between 1% and 2% (Xu et al., 2014). However, observations (e.g.

the Gunn–Peterson trough or Gunn-Peterson test) have shown that the universe is fully ionized

(with a neutral fraction below 10−4) at z ∼ 6. Therefore, stellar radiation from galaxies must be

11



the dominant source of radiation responsible for cosmic reionization. Two important questions

that follow are (1) How abundant are galaxies as a function of luminosity and redshift? (2) What

fraction of ionizing photons escaped from theses galaxies into the IGM?

Recent observations have provided valuable constraints on the nature of the first galaxies

and their role in the EoR. The Hubble Space Telescope Ultra Deep Field (Ellis et al., 2013) and

Frontier Fields (Coe et al., 2015) projects probe galaxies with stellar masses as small as 107M⊙

at z ≳ 6 and galaxies as early as z ≳ 11. There could be an unseen population of even fainter

and more abundant galaxies that may eventually be detected by JWST.

In the calculations of reionization, the key quantity is the ionizing emissivity, ργ (in units of

erg s−1Hz−1Mpc−3). The intrinsic ionizing emissivity can be inferred by integrating the product

of the luminosity function ϕ(L) and star formation rate (SFR), given a relation between total

luminosity L and SFR. To get ργ , we need another parameter, namely the fraction of ionizing

photons that escape into IGM, fesc. It is arguably the most uncertain parameter in the models of

reionization.

By measuring the hydrogen ionizing photon budget and constraining the best-fit Schechter

parameters for the galactic luminosity function, Ouchi et al. (2009) suggests that galaxies at

z = 7 need large fesc (≳ 0.2) to keep the universe ionized. Khaire et al. (2016) use cosmological

radiation simulations and find that the updated QSO emissivity and star formation history have

similar implications on the fesc at z > 5.5. This value is too large with respect to what is observed

in local galaxies.

A number of attempts have been made to calculate the escape fraction of hydrogen LyC

photons from galaxies using analytic models and simulations of galaxy formation (Ricotti &

Shull, 2000; Gnedin et al., 2008; Wise et al., 2014; Ma et al., 2015; Xu et al., 2016). However,

12



because of the complexity of the problem and the uncertainty about the properties of the sources

of reionization, the results are inconclusive. Furthermore, any realistic theoretical estimate of

fesc must take into account the escape fraction of LyC radiation from the molecular clouds in

which the stars are born, ⟨fMC
esc ⟩, a sub-grid parameter in galaxy-scale and in cosmological-scale

simulations.

Typically ⟨fMC
esc ⟩ is set to unity in these simulations, which could dramatically overpredict

the galactic fesc (Ma et al., 2015). Observational constraints of ⟨fMC
esc ⟩ are very limited and suffer

from significant uncertainty (Doran et al., 2013). Therefore, a calculation of ⟨fMC
esc ⟩ from simula-

tions of star formation from molecular clouds is essential to the understanding of the EoR, which

is the motivation of the work presented in Chapter 3.

1.3 Techniques for Simulating Star Formation

The analytic theory is essential for understanding the big picture for star formation, but

only numerical models can capture the nonlinear behaviors of the turbulent collapse of a GMC

and take into account the effects of magnetic fields. On this scale, fluid dynamics and gas cooling

play a dominant role in structure evolution. The complexity of hydrodynamic processes such as

shock heating, atomic radiation cooling, and star formation requires an accurate treatment of the

gas and radiation components.

For the sake of completeness, we give a short discussion of the MHD equations. The

13



standard form of the ideal MHD equations is

∂ρ

∂t
+∇ · (ρv) = 0, (1.7)

ρ

[
∂v

∂t
+ (v ·∇)v

]
= −∇P +

(∇×B)×B

4π
, (1.8)

ρ

[
∂e

∂t
+ (v ·∇)e

]
= −P (∇ · v)− ρL, (1.9)

∂B

∂t
= ∇× (v ×B), (1.10)

where ρ is the mass density, v is the fluid velocity, e is the specific total energy, and L represents

the net loss function that describes the radiative heating and cooling of the gas. To close the

system of equations, we must complement it with an equation of state. For the ISM, a perfect gas

is a good assumption, with P = (γ − 1)ρϵ, where γ is the adiabatic index and ϵ is the specific

energy excluding kinetic energy of the gas.

To get a better physical understanding of the MHD equations and the role of magnetic

fields, we write the Lorentz force as

fL =
(∇×B)×B

4π
= −∇

(
B2

8π

)
+

(B ·∇)B

4π
. (1.11)

The first term is the magnetic pressure gradient and the second term is called the magnetic tension.

In computational astrophysics, two main techniques are employed to solve hydrodynamical

equations: smoothed-particle hydrodynamics (SPH) and grid-based methods. The fundamental

difference is how the fluid is discretized. An SPH solver divides fluid into mass elements (par-

ticles), whose evolution is then followed, in line with the Lagrangian approach. Grid methods,

however, subdivide the computational domain into volume elements (cells) that are fixed in space,

14



following the Eulerian formulation of hydrodynamics.

The natural outcome of an SPH code is it automatically adapts the numerical resolution

in dense regions, with more particles concentrated in those areas, which makes it a valid tool in

situations with large density contrasts such as self-gravitating fluids. Grid-based methods can

increase their spatial resolution using the adaptive mesh refinement (AMR) technique, introduc-

ing new grid elements from the partition of a cell into smaller cells triggered by some specified

refinement criteria.

RAMSES (Teyssier, 2002) is a massively parallel grid-based adaptive mesh refinement code,

originally developed to study the co-evolution of dark matter and gas in a cosmological context.

RAMSES is now a complete tool for simulations of self-gravitating fluid dynamics, equipped

with a magneto-hydrodynamic (MHD) solver (Fromang et al., 2006) and several modules which

implement AGN feedback, star formation recipes (Bleuler & Teyssier, 2014), interstellar medium

cooling functions, and stellar feedback (Rosdahl et al., 2013; Rosdahl & Teyssier, 2015). More

physics including stellar evolution, photoionized metal chemistry, and dust dynamics are recently

implemented (Geen et al., 2021; Katz, 2022; Moseley et al., 2023).

The AMR structure implemented in RAMSES is tree-based. Every cell can be refined into

2dim child cells called octs. Each oct at refinement level l points to all the child cells in the next

refinement level (l+ 1), its parent cell in level l− 1, and all the sibling cells of the parent cell. A

cell can be refined according to several criteria of the user’s choice. The refinement criteria used

in this thesis work is the Jeans length criteria, stating that the local Jeans length must be resolved

by at least Ns cells, where Ns = 10.

The formation of stars is modeled with sink particles, following the implementation of

Bleuler & Teyssier (2014). When the density is higher than a certain value, the Jeans length

15



criterion (with a smaller Ns) is violated at the highest refinement level and the high-density peak

is replaced with a sink particle. Sink particles are allowed to merge and accrete gas according

to different merging or accretion schemes (threshold, Bondi, flux). The dynamical interaction

between sinks can be computed using direct summation and the force between a sink and gas can

be computed using the same technique used for other particles (particle-mesh method).

Rosdahl et al. (2013) have integrated radiation hydrodynamics into RAMSES by enclosing

the radiation transfer (RT) equations using the M1 approximation of the Eddington tensor and

then solving the differential equations using the first-order Godunov scheme. This moment-based

approach has the advantage of being independent of the number of radiative sources and is much

more computationally efficient compared to ray-tracing techniques.

Combining all these techniques presented above, RAMSES-RT has become a highly sophis-

ticated and powerful tool in the field of computational astrophysics for cosmological simulations

as well as small-scale star formation simulations.

1.4 Thesis Outline

This thesis investigates star formation processes from GMC to protostellar disk scales. In

Chapter 2, I explore the physics and laws of star formation from molecular clouds based on a set

of radiation-MHD simulations. Next, I postprocess the simulation data and calculate the escape

of LyC photons from GMCs into the intercloud medium that takes into account contributions

from spatially resolved stars. This part is reported in Chapter 3. In Chapter 4, I investigate the

fragmentation of massive prestellar cores and the formation of circumstellar disks using a series

of extremely high-resolution simulations of prestellar cores. Then, in Chapter 5, I explore the

16



magnetic braking problem in disk formation and propose a solution to the “magnetic braking

catastrophe”. As part of an ongoing project on stellar dynamics, I study star clustering by quan-

tifying the clumpiness of the stars using the two-point correlation function ξ(r) in Chapter 6.

Finally, in Chapter 7, I summarize the main conclusions and introduce future research directions

building on the work presented in this thesis.

1.5 A summary of software and facilities

Computational resources used in this dissertation:

1. The Deepthought2 HPC cluster of the University of Maryland supercomputing resources2

2. The Yorp Cluster of the Department of Astronomy, University of Maryland3

3. Student workstations of the Department of Astronomy, University of Maryland4

Software/codes used:

1. RAMSES-RT (Teyssier, 2002; Bleuler & Teyssier, 2014; Rosdahl et al., 2013)

2. NUMPY (van der Walt et al., 2011), MATPLOTLIB (Hunter, 2007), and YT (Turk et al.,

2011)

3. RAMTOOLS5

4. HEALPIX6 (Zonca et al., 2019), HEALPY (Górski et al., 2005)

2http://hpcc.umd.edu
3https://www.astro.umd.edu/twiki/bin/view/AstroUMD/YorpCluster
4https://www.astro.umd.edu/twiki/bin/view/AstroUMD/StudentWorkstations
5https://chongchonghe.github.io/ramtools-pages/
6http://healpix.sf.net

17

http://hpcc.umd.edu
https://www.astro.umd.edu/twiki/bin/view/AstroUMD/YorpCluster
https://www.astro.umd.edu/twiki/bin/view/AstroUMD/StudentWorkstations
https://chongchonghe.github.io/ramtools-pages/
http://healpix.sf.net


Chapter 2: Simulating Star Clusters Across Cosmic Time: I. Initial Mass Func-

tion, Star Formation Rates and Efficiencies

In this chapter, we present radiation-magneto-hydrodynamic simulations of star forma-

tion in self-gravitating, turbulent molecular clouds, modeling the formation of individual mas-

sive stars, including their UV radiation feedback. The set of simulations have cloud masses

between mgas = 103 M⊙ to 3 × 105 M⊙ and gas densities typical of clouds in the local universe

(ngas ∼ 1.8 × 102 cm−3) and 10× and 100× denser, expected to exist in high-redshift galax-

ies. The main results are: i) The observed Salpeter power-law slope and normalisation of the

stellar initial mass function at the high-mass end can be reproduced if we assume that each star-

forming gas clump (sink particle) fragments into stars producing on average a maximum stellar

mass about 40% of the mass of the sink particle, while the remaining 60% is distributed into

smaller mass stars. Assuming that the sinks fragment according to a power-law mass function

flatter than Salpeter, with log-slope 0.8, satisfies this empirical prescription. ii) The star forma-

tion law that best describes our set of simulations is dρ∗/dt ∝ ρ1.5gas if ngas < ncri ≈ 103 cm−3,

and dρ∗/dt ∝ ρ2.5gas otherwise. The duration of the star formation episode is roughly 6 cloud

sound crossing times (with cs = 10 km/s). iii) The total star formation efficiency in the cloud is

f∗ = 2%(mgas/10
4 M⊙)0.4(1 + ngas/ncri)

0.91, for gas at solar metallicity, while for metallicity

Z < 0.1 Z⊙, based on our limited sample, f∗ is reduced by a factor of∼ 5. iv) The most compact

18



and massive clouds appear to form globular cluster progenitors, in the sense that star clusters

remain gravitationally bound after the gas has been expelled.

2.1 Introduction

Star formation in galaxies is a complex and only partially understood astrophysical phe-

nomenon. It is difficult to formulate a general theory in part because of the wide range of scales

and of physical processes involved. From an observational point of view, quantifying star forma-

tion efficiency (SFE) in nearby molecular clouds has been the focus of much recent research (e.g.,

Lada et al., 2010; Heiderman et al., 2010; Gutermuth et al., 2011). A power-law relationship be-

tween the gas surface density of galaxies and their star formation rate (SFR) was first proposed by

Schmidt (1959) and later tested by large, multi-galaxy data (Kennicutt, 1998). This relationship

has been widely used in cosmological simulations of galaxy formation. However, on sub-galactic

scales the dispersion of star formation rates for a given gas surface density of H I is large, and

other parameters such as gas metallicity (Bolatto et al., 2011; Krumholz, 2013) and stellar surface

density (Leroy et al., 2008) appear to become important. As the resolution of surveys improved,

numerous studies have shown that star formation on kpc scales is more strongly correlated with

H2 surface density (e.g. Krumholz, 2014) rather than atomic gas. Therefore, modern cosmolog-

ical simulations of galaxy formation aim at reproducing the molecular phase of the interstellar

medium (ISM) and adopt an empirical sub-grid recipe for star formation within partially resolved

molecular clouds of the form ρ̇∗ ∝ ρnH2
, where typically n = 1 or 1.5. The molecular phase of

the ISM is treated in simulations using different prescriptions: Robertson & Kravtsov (2008) pre-

computed a grid of models from a photo-chemistry code, Gnedin et al. (2009) directly solved the

19



formation and dissociation equations for H2, but with an increased formation rate to model unre-

solved clumping, and Kuhlen et al. (2012) used an analytic model to estimate the equilibrium H2

abundance. The sub-grid recipe (with grid maximum resolution typically between few parsecs to

few kpc) is calibrated to reproduce observational data in galaxies at z = 0.

However, the conditions in the ISM of high-redshift galaxies are likely different to those

found in the present day. Krumholz et al. (2012) argue that the SFR is in fact correlated to

the local free-fall time set by the gas density, not to the column density. Simulations show that

densities and pressures of star-forming regions in high-redshift galaxies are much higher than

in today’s ISM (e.g., Ricotti, 2002; Wise et al., 2014; Ricotti, 2016). Using Adaptive Mesh

Refinement (AMR) simulations of the first stars and galaxies with parsec-resolution, Ricotti et al.

(2016) found that compact molecular clouds in primordial galaxies can either form gravitationally

bound star clusters that resemble the progenitors of today’s globular clusters (GC)1, or the clusters

may disperse and fill up a large fraction of the dark matter halo of primordial dwarf galaxies. In

this second case the stars would appear as spheroids 20-200 pc in radius, dark matter dominated

and with very low surface brightness. These objects would be identified today as “ultra-faint”

dwarf galaxies observed in the Local Group (e.g., Willman et al., 2005; Zucker et al., 2006a,b;

Belokurov et al., 2007; Walsh et al., 2007; Majewski et al., 2007; Martin et al., 2009). Star

formation in compact star clusters appears to be especially important, even perhaps the dominant

mode of star formation at high-redshift. Thus, in order to make progress in understanding the

formation of the first dwarf galaxies and the sources of reionisation, it is important to focus on

understanding the small-scale physics of this process, which is poorly resolved in cosmological

1The compact bound stellar objects found in the simulations are actually not only globular cluster progenitors,
but also ultra-compact dwarfs and dwarf-globular transition objects, depending on whether the stellar cluster forms
at the center of the halo, in the disk’s spiral arms, or even in satellite minihalos.

20



simulations.

Most numerical work on star formation in molecular clouds focuses on star formation in

the local universe, aiming at explaining observed young star forming regions. In this chapter

we analyse the results of a large grid of simulations of realistic molecular clouds with initial

conditions chosen to reproduce not only local molecular clouds but also clouds that form in

higher density and pressure environments, typical of star formation in high redshift galaxies. We

vary the masses of the clouds, their compactness (central density), and in few cases explore the

effect of changing the gas metallicity and therefore the gas cooling function.

The motivation for this chapter is twofold. The first goal is to deepen our understanding

of the physics of star formation in high-pressure environments to justify and inform the sub-grid

star formation recipe used in cosmological simulations. A closely related important question in

Near Field Cosmology is: how does the formation of self-gravitating bound star-clusters relate to

the star formation efficiency, compactness, mass and gas metallicity of molecular clouds found in

cosmological simulations? We will only touch on this questions in the present chapter, but more

detailed work will be presented in a followup work.

The second goal is to estimate the escape fraction of H I ionising radiation from molecular

clouds as a function of cloud compactness and mass. This is the first necessary step for a realistic

estimate of the escape fraction from galaxies. Ricotti (2002) have shown that, if a non-negligible

fraction of today’s GCs formed at z > 6 with ⟨fMC
esc ⟩ ∼ 1, their progenitors would be a dominant

source of ionising radiation during reionisation. Katz & Ricotti (2014) presented arguments

in support of significant fraction of today’s old GCs forming before the epoch of reionisation.

However, although it is naively expected, it has not been shown with numerical simulations that

⟨fMC
esc ⟩ from GC progenitors forming in compact molecular clouds is higher than ⟨fMC

esc ⟩ in more

21



diffuse clouds. The answer to this question and the contribution of compact star clusters to

reionisation will be presented in the next chapter.

This chapter is organized as follows. In Section 2.1.1 we present a brief review of the

current status of numerical simulations of star cluster formation. In Section 2.2 we provide an

overview of our numerical methods, including details on the initial conditions of our simulations

and the recipes for formation of sink particles and feedback. In Section 2.3 we present some

results from the analysis of our large set of simulations with emphasis on the stellar initial mass

function (IMF) and the star formation rate (SFR) and efficiency (SFE). A summary and conclu-

sions are presented in Section 2.4.

2.1.1 The IMF and SFE of Molecular Clouds

Simulations of molecular cloud dynamics are valuable tools in understanding the condi-

tions in the ISM. Typically these simulations adopt idealized initial conditions similar to those in

observed clouds: a gas cloud∼ 1−10 pc in size supported against gravity by a turbulent velocity

field such that the initial virial ratio, i.e. the ratio of the kinetic energy to the potential energy of

the cloud, is ≲ 0.5. One model involves injecting turbulence into a volume of gas in the initial

conditions and allowing it to decay over time. This can be done by either using smoothed parti-

cle hydrodynamics (SPH, e.g. Klessen, 2001; Bonnell et al., 2006) or grid-based methods (e.g.,

Gammie et al., 2003). Another model involves adding turbulence continuously over time, simu-

lating the effect of momentum injection from outside flows or energy from massive stars inside

the cloud (Vazquez-Semadeni et al., 1997; Ballesteros-Paredes et al., 2006; Padoan et al., 2007).

Many of these models adopted an isothermal equation of state, while others have included self-

22



consistent cooling and heating functions (e.g. Koyama & Inutsuka, 2004; Audit & Hennebelle,

2005) and molecular chemistry (e.g. Glover et al., 2010).

The fragmentation of molecular clouds into stars is a long-standing problem. Observational

studies (Salpeter, 1955; Kroupa, 2002; Chabrier, 2005) have found that the masses of stars follow

a “Initial Mass Function” (IMF) with a power law dN/ d logM ∝ M−Γ at the high-mass end

(Salpeter 1955 calculate Γ ≈ 1.35). Various theoretical models have been constructed to explain

this (Padoan & Nordlund, 2002; Mac Low & Klessen, 2004; Hennebelle & Chabrier, 2008; Hop-

kins, 2012a) based on gravoturbulent fragmentation of the host cloud. Radiative stellar feedback

has been invoked to explain the precise shape of the IMF, using both simulations (e.g., Bate,

2009b) and analytic models (e.g., Guszejnov & Hopkins, 2016). Early pioneering simulations

of cluster formation approached the problem of producing a well-defined IMF (e.g., Bate et al.,

2003; Bate & Bonnell, 2005; Klessen et al., 2008; Offner et al., 2008), but were often limited in

terms of statistics or resolution. More recent work, with increasing computing power, provided

more reliable statistics and IMF distributions (e.g., Bonnell et al., 2003; Bate, 2009b; Bonnell

et al., 2011; Girichidis et al., 2011; Krumholz et al., 2011; Bate, 2012; Ballesteros-Paredes et al.,

2015).

Simulations attempting to capture the stellar IMF require a high dynamic range to resolve

both brown dwarfs and OB stars. Most recently, Bate (2019) resolve in detail the mass spectrum

of brown dwarfs while only producing stars of up to 3 M⊙, finding that low metallicities do not

produce observable differences in the stellar IMF, while increasing fragmentation. Gavagnin et al.

(2017) have lower mass resolution but capture more massive stars that emit significant quantities

of ionising radiation, arguing that this alters the high mass end of the IMF. In the absence of

radiation and cooling, Lee & Hennebelle (2018b) and Lee & Hennebelle (2018a) study the early

23



formation of protostellar Larson cores, and find that the choice of equation of state (eos) has

a strong influence on the peak of the IMF. In general, these works are relatively successful at

reproducing not only the IMF but also stellar multiplicity and separation.

Previous authors have included ideal MHD in their simulations (Myers et al., 2013; Krumholz

et al., 2016; Cunningham et al., 2018). However, since these authors only form stars up to ∼20

M⊙, they neglect ionising radiation. Non-ideal MHD effects, while challenging to include in

resolving the IMF for reasons of computational cost, appear to affect the dynamics of protostar

formation on small scales (Masson et al., 2016; Vaytet et al., 2018). The physics that shapes

the IMF is complex, and a full treatment that covers non-ideal MHD, both low and high energy

radiation, chemistry and the full mass range of stars remains difficult with modern computational

resources.

As well as the stellar IMF, an important consideration is how many stars are formed out of

a given mass of gas, or the Star Formation Efficiency (SFE). The efficiency of conversion of gas

into stars is typically much lower than 100% since energetic processes from massive stars are able

to disperse the cloud in which a star cluster forms before all of the gas collapses into protostars.

These processes are widely termed “feedback” (see review by Dale et al., 2015). Recent work

favours ionising radiation as the main driver of molecular cloud dispersal (Dale et al., 2005;

Gritschneder et al., 2009; Peters et al., 2010; Walch et al., 2012; Dale et al., 2012), as opposed to

other effects such as stellar winds (Dale et al., 2014), although Howard et al. (2016) find that UV

photoionisation has little effect on the initial evolution of the SFE.

The relationship between gas properties and the SFE is a matter of ongoing study. Lada

et al. (2010) and Heiderman et al. (2010) argue for a constant ratio between SFE and cloud mass

for gas above a certain surface density, although this is still subject to discussion (Gutermuth

24



et al., 2011; Hony et al., 2015). There is no clear theoretical link between the SFE and projected

column density, although Clark & Glover (2014) argue that there may be a link between the

observed column density and the local density around the star. Geen et al. (2017) reproduce the

SFE observed by Lada et al. (2010), although they find that the result is likely to be dependent

on the average density of the neutral gas in the cloud. These simulations produce similar results

to the simulations of Colin et al. (2013). Geen et al. (2018) finds that SFE can change by up

to a factor of 4 by varying the initial velocity field of the cloud and the stellar IMF, although

relationships can be found between the early cloud state and the final SFE. Semi-analytic models

by Vázquez-Semadeni et al. (2018) also find considerable scatter in the SFE.

2.2 Numerical Simulations and Methods

We conduct our numerical simulations using the AMR radiative magneto-hydrodynamical

code RAMSES (Teyssier, 2002; Bleuler & Teyssier, 2014). Radiative transfer is implemented us-

ing a first-order moment method with M1 closure described in Rosdahl et al. (2013). Kim et al.

(2017) demonstrates that M1 closure method is inaccurate near sources only in regions where the

flux is about an order of magnitude smaller than the mean value (due to shielding), while it agrees

with adaptive ray-tracing methods (e.g., Wise et al., 2014; Hartley & Ricotti, 2016) both at larger

distances from individual sources and on global scales. M1 closure, however, is significantly

more computationally efficient than ray-tracing methods. The ionising photons interact with neu-

tral gas and we track the ionisation state and cooling/heating processed of hydrogen and helium

(see Geen et al., 2017, for details). Our simulations include magnetic fields in the initial condi-

tions, but we do not include the chemistry of molecular species (i.e., formation/dissociation). 3-D

25



‘zoom-in’ simulations of the chemical evolution of molecular clouds suggest that, for gas at solar

metallicity, the cloud is almost fully molecular with H2 fractions around 0.9 in the later stages of

transition to dense molecular phase (Seifried et al., 2017).

We simulate a set of isolated and turbulent molecular clouds that collapse due to their own

gravity. We explore a grid of simulations varying the initial gas mass and compactness (i.e., the

core density) of the clouds. In our simulations, dense proto-stellar cores collapsing below the

resolution limit of the simulations produce sink particles. These sinks may represent single stars

or multiple stars or even clusters of stars if the resolution is not sufficiently high. However, in all

our simulations we aim at reproducing a realistic high-mass end of the stellar IMF and therefore

realistic feedback from individual massive stars. To accomplish this goal, sink particles emit

hydrogen and helium ionising photons according to their mass as described in § 2.2.3. The gas

is ionised and heated by massive stars, producing over-pressurised bubbles that blow out the gas

they encounter. In our simulations low mass stars and proto-stellar cores do not produce any

feedback. In this work we do not include mechanical feedback from supernova (SN) explosions

and from stellar winds and we also neglect the effect of radiation pressure from infrared radiation.

However, with the exception of the two most massive clouds in the set of simulations representing

today’s molecular clouds (the lowest density set), all the simulations stop forming stars before the

explosion of the first SN. Therefore, neglecting SN feedback is well justified in these cases. We

find that in all simulations star formation has ceased before ∼ 5− 6 tff , which is the typical time

it takes for feedback to act. We stop simulations no earlier than this point. A few simulations are

continued beyond this time. This does not have an effect on the IMF since the mass function of

sink particles does not change after star formation has ceased.

For the simulations in the set in which SN explosions should occur while star formation is

26



ongoing, in order to compensate for this missing feedback, we do not shut down UV radiation

feedback from massive stars after the time the star should have exploded as SN. In the following

sections we provide some more details on the simulations set up.

2.2.1 Initial Conditions

We run a grid of 14 simulations of clouds with a range of central densities and initial gas

masses. We also run some additional simulations varying the initial gas metallicity and therefore

the gas cooling function. The magnetic field strength in the initial conditions is set such that

va = 0.2 σ3D, where va is Alfven wave velocity and σ3D is the turbulence velocity dispersion.

This νa is ∼ 2 times smaller than that measured in a group of molecular clouds by Crutcher

(2012) who finds νa ≈ 0.5σ3D.

The clouds have initially a spherically symmetric structure with density profile of a non-

singular isothermal sphere with core density nc. The cloud extends out to rgas = 3rc, where rc

is the core radius. Beyond rgas the cloud is embedded in a uniform density envelope that extends

to 6rc with a density 0.01 nc. Outside of the envelope the number density is constant at 1 cm−3.

The box length Lbox is set to 48rc in each simulation. The initial value of the (isothermal) sound

speed of the cloud is set to cs = 0.24 km/s, while the envelope and background densities are in

pressure equilibrium.

The initial density profile is perturbed with a turbulent velocity field, analogously to the set

up used in Geen et al. (2017). The initial turbulence of the clouds follows a Kolmogorov power

spectrum with random phases and has an amplitude such that the cloud is approximately in virial

27



equilibrium. All simulations have the same set of random phases. The initial cloud virial ratio

αvir =
5σ2

3DR

3GM
≈ 0.4, (2.1)

is kept constant in all the simulations. Therefore the ratio tff/tturb, where tturb ≡ R/σ3D is kept

constant in all the simulations. However, the sound crossing time tcr ≡ R/cs, where throughout

this chapter we assume cs = 10 km/s, is not constant. The virial parameter, αvir, is small enough

to ensure collapse and fragmentation, but sufficiently large to prevent a rapid radial collapse

of the cloud. Before allowing any star formation in the cloud we evolve these idealized initial

conditions for ∼ 3tff , so that the turbulent velocities develop into density perturbations and the

initial conditions relax into a quasi-equilibrium turbulent medium. If we do not allow the initial

conditions to relax before forming stars, the stars form mostly near the center of the cloud during

the transient relaxation phase.

28



Ta
bl

e
2.

1:
In

iti
al

co
nd

iti
on

s
of

ou
r

16
si

m
ul

at
io

ns
.

(a
)

M
ea

n
nu

m
be

r
de

ns
ity

of
th

e
cl

ou
d,

ex
cl

ud
in

g
th

e
en

ve
lo

pe
.

T
he

co
re

de
ns

ity
is

∼
5

tim
es

hi
gh

er
.(

b)
T

he
gl

ob
al

fr
ee

-f
al

lt
im

e
of

th
e

cl
ou

d
(t

f
f
≡

3√
3
π

3
2
G
ρ
c
≈

1.
3√

3
π

3
2
G
ρ
).

(c
)M

ax
im

um
le

ve
lo

fr
efi

ne
m

en
t.

(d
)I

ni
tia

l
cl

ou
d

m
as

s,
ex

cl
ud

in
g

th
e

en
ve

lo
pe

.(
e)

T
he

na
m

e
of

ea
ch

cl
ou

d
us

ed
th

ro
ug

ho
ut

th
e

ch
ap

te
r.

Se
e

Se
c.

2.
2.

1
on

ho
w

th
ey

ar
e

de
fin

ed
.(

f)
In

iti
al

cl
ou

d
ra

di
us

,e
xc

lu
di

ng
th

e
en

ve
lo

pe
.(

g)
M

ax
im

um
sp

at
ia

lr
es

ol
ut

io
n.

(h
)D

en
si

ty
th

re
sh

ol
d

fo
rs

in
k

fo
rm

at
io

n.
(i

)J
ea

ns
m

as
s

at
th

e
si

nk
de

ns
ity

th
re

sh
ol

d.
(j

)T
ur

bu
le

nc
e

M
ac

h
nu

m
be

r.
(k

)S
ou

nd
cr

os
si

ng
tim

e
r g

a
s/
c s

fo
rc

s
=

10
km

/s
.(

l)
M

et
al

lic
ity

of
th

e
ga

s
us

ed
in

th
e

co
ol

in
g

fu
nc

tio
n,
Z

=
[F

e/
H

].
(m

)T
hi

s
se

tu
p

ha
s

2
ex

tr
a

si
m

ul
at

io
ns

w
ith

lo
w

er
m

et
al

lic
iti

es
be

si
de

s
on

e
w

ith
sa

m
e

m
et

al
lic

ity
as

al
lo

th
er

on
es

.
Se

e
Se

c.
2.

3.
4.

(x
)T

he
m

ea
n

su
rf

ac
e

de
ns

ity
in

a
sq

ua
re

of
th

e
si

ze
of

th
e

cl
ou

d
ra

di
us

.
(x

)E
sc

ap
e

ve
lo

ci
ty

at
th

e
cl

ou
d

ra
di

us
of

th
e

in
iti

al
cl

ou
d.

m
g
a
s(
M

⊙
)

1.
0
×

10
3

3.
2
×
10

3
1.
0
×
10

4
3.
2
×
10

4
1.
0
×
10

5
3.
2
×
10

5

C
lo

ud
N

am
e

X
S-

F
S-

F
M

-F
L

-F
X

L
-F

n
g
a
s
=

1.
8
×
10

2
cm

−
3

r g
a
s

(p
c)

5.
0

7.
3

11
16

23
Σ

(M
⊙

pc
−
2
)

41
61

89
13

1
19

3
v e

sc
(k

m
/s

)
2.

3
3.

4
5.

1
7.

4
11

∆
x
m
in

(A
U

)
50

0
73

0
11

00
16

00
23

00
t ff

=
4.
4

M
yr

n
si
n
k

(c
m

−
3
)

1.
2
×
10

7
5.
6
×
10

6
2.
6
×
10

6
1.
2
×
10

6
5.
6
×
10

5

M
J
(M

⊙
)

0.
3

0.
4

0.
6

0.
9

1.
3

l m
a
x
=

15
M

4.
6

6.
8

10
15

22
t c

r
(M

yr
)

0.
5

0.
7

1.
1

1.
5

2.
3

Z
/Z

⊙
1

1
1

1
1

C
lo

ud
N

am
e

X
S-

C
S-

C
M

-C
L

-C
,L

-C
-l

m
,L

-C
-x

lm
n
g
a
s
=

1.
8
×
10

3
cm

−
3

r g
a
s

2.
3

3.
4

5.
0

7.
3

Σ
19

3
28

3
41

5
60

9
v e

sc
3.

4
5.

1
7.

4
11

∆
x
m
in

46
0

68
0

10
00

15
00

t ff
=

1.
4

M
yr

n
si
n
k

1.
4
×
10

7
6.
5
×
10

6
3.
0
×
10

6
1.
4
×
10

6

M
J

0.
3

0.
4

0.
6

0.
8

l m
a
x
=

14
M

6.
8

10
15

22
t c

r
0.

23
0.

33
0.

5
0.

7
Z

1
1

1
1,
1/
10

,1
/4
0

C
lo

ud
N

am
e

X
X

S-
V

C
X

S-
V

C
S-

V
C

M
-V

C
L

-V
C

n
g
a
s
=

1.
8
×
10

4
cm

−
3

r g
a
s

0.
7

1.
1

1.
6

2.
3

3.
4

Σ
60

9
89

4
13

12
19

25
28

27
v e

sc
3.

4
5.

1
7.

4
11

16
∆
x
m
in

15
0

22
0

32
0

46
0

68
0

t ff
=

0.
44

M
yr

n
si
n
k

1.
4
×

10
8

6.
5
×
10

7
3.
0
×
10

7
1.
4
×
10

7
6.
5
×
10

6

M
J

0.
08

0.
12

0.
17

0.
26

0.
38

l m
a
x
=

14
M

7
10

15
22

32
t c

r
0.

07
0.

10
0.

15
0.

23
0.

33
Z

1
1

1
1

1

29



XS-F S-F M-F L-F XL-F

XS-C S-C M-C L-C

XXS-VC XS-VC S-VC M-VC L-VC

3 4 5
log gas mass (M )

2

3

4

5

log
 g

as
 d

en
sit

y (
cm

3 ) Bonnell03

Bonnell11

Bate19
Ballesteros15

Bertelli16

Jones18

Jones18

Lee18

Lee18

Gavagnin17

Figure 2.1: Simulation parameters in this work (colored ovals) compared to previous works
(stars). The parameter space considered here is mass of the gas cloud (x-axis) versus mean
particle number density of the cloud (y-axis). The labels showing the previous work found in the
literature include: Bonnell et al. (2003, 2011); Ballesteros-Paredes et al. (2015); Bertelli Motta
et al. (2016); Gavagnin et al. (2017); Jones & Bate (2018); Lee & Hennebelle (2018a); Bate
(2019).

A detailed list of the parameters in our simulations is shown in Table 2.1. The clouds in

Table 2.1 are labelled with letters of two or three parts. The first part is either ‘XXS’ (extra-

extra-small), ‘XS’ (extra-small), ‘S’ (small), ‘M’ (medium), ‘L’ (large), or ‘XL’ (extra-large),

representing various initial gas masses of 103, 3.16×103, 104, 3.16×104, 105, and 3.16×105 M⊙,

respectively. The second part is either ‘F’ (fiducial, which are the most similar to clouds in the

solar neighbourhood), ‘C’(compact), or ‘VC’ (very compact), in order of increasing initial mean

gas density. The mean particle number density of the cloud, ngas = ρgas/(µmp), where µ = 1.4

is the mean molecular weight of the atomic gas, increases by a factor of ten between each set of

simulations from 1.8× 102 cm−3 to 1.8× 104 cm−3. The ‘L-C’ setup has two more simulations

that have a third part in the name, ‘lm’ and ‘xlm’, representing ‘low-metallicity’ and ‘extra-low-

metallicity’. A comparison of our setups with the literature is shown in Figure 2.1.

30



2.2.2 Resolution and Sink Formation

We use a Cartesian grid with an octree structure with cells that we subdivide into 23 child

cells as the simulation evolves (“adaptive refinement”). Our starting refinement level is ℓmin = 7

(corresponding to ∆x = Lbox/2
7) and maximum level of refinement is ℓmax = 15 for runs

with the lowest mean density and ℓmax = 14 for all the other runs. The resolution is therefore

∆xmin = Lbox/2
15 for the ”fiducial” runs, which corresponds to resolutions between 500 AU and

2300 AU. The ”compact” clouds have resolution between 460 AU and 1500 AU and the ”very

compact” clouds between 150 AU and 680 AU.

In order to resolve the Jeans length with N grid cells it is required that

λJ = cs

√
π

Gρ
> Nsink∆x. (2.2)

From Equation (2.2), the Jeans length is resolved with at least Nsink grid points if ρ < ρJ , where

ρJ =
πc2s

GN2∆x2
. (2.3)

In our simulations we enforce the refinement criterion that the Jeans length is resolved with at

least Nref = 10 cells. Hence, when the local density goes up and reaches a point where λJ

becomes smaller than Nref∆x, each cell is refined individually into eight new child cells. This

refinement condition is always true, up to the maximum refinement level (when ∆x = ∆xmin).

When the gas density exceeds ρmax
J = ρJ(∆x = ∆xmin, N) at the maximum refinement level,

we cannot continue to resolve the Jeans length with at least Nref cells. We therefore create sink

31



particles to trace material above these densities. We set ρsink = ρJ(∆x = ∆xmin, N = Nsink)

as critical density threshold to form sink particles. Sink particles are created on the fly using a

peak detection algorithm (see Bleuler & Teyssier, 2014, for details on sink particle formation in

RAMSES). We first detect density clumps above a density threshold fcρsink, with fc = 0.1. Then,

the algorithm performs a peak density check, a collapsing check (∇ · v = 0), and virial check

before forming a sink particle.

In order to avoid numerical fragmentation it is usually suggested that Nsink ≥ 4 (Truelove

et al., 1997). In our simulations we adopt Nsink = 5 for reasons detailed in Appendix A.1. With

an initial sound speed cs = 0.24 km/s the Jeans mass at the sink density threshold is

MJ =
4π

3
ρJ

(
λJ

2

)3

= 0.55M⊙

(
∆xmin

1000AU

)
, (2.4)

which results in MJ ∼ 0.08M⊙–0.8M⊙ for the compact and very compact clouds and∼ 0.3M⊙–

1.3M⊙ for the fiducial clouds.

The sink particles are then treated like point masses and accrete gas based on the mech-

anism described as ‘threshold accretion’ in Bleuler & Teyssier (2014). The dynamics of the

sink particles takes into account gravitational force from gas and stars and it is evolved using a

leap-frog integration scheme. The effect of gas dynamical friction is not included.

2.2.3 Feedback and Properties of UV source

In our simulations, ionising UV photons are emitted from sink particles from the time they

form to the end of the simulation. Massive stars have lifetime of few Myrs, shorter than the

duration of some of our simulations, and they may explode as SNe during the simulation. Since

32



we are not implementing SNe feedback, we keep the stars emitting radiation after their death

to compensate for the lack of SNe in the attempt of avoiding underestimating feedback effects.

While SN explosions produce a significant amount of mechanical energy (typically 1051 egs),

the energy associated with ionising radiation from massive stars integrated through their main-

sequence lifetime is comparable (or larger for more massive stars) and this feedback starts acting

earlier than SN feedback. For an O-star, more than half of the radiation is emitted in hydrogen

ionising photons. Typically ∼ 10% of a star’s hydrogen is burned in the nuclear fusion process,

with an energy efficiency of ∼ 0.7%. Thus the amount of energy radiated by a massive star

during its life time is ∼ 2× 10−3M∗, or ∼ 4× 1052 ergs for a 20 M⊙ star.

For each simulation we estimate the total hydrogen-ionising photon emission rate at a given

time as Scl(mcl) = 8.96× 1046 s−1 (mcl/M⊙) (see Geen et al., 2017), where mcl is the total mass

of the sink particles. This is calculated by Monte Carlo sampling a stellar population as described

in Geen et al. (2016) (See Sec. A.2). The fraction of the total hydrogen ionising photon emission

rate attributed to each sink particle is based on the following relation

q(mi) = V (0.3mi)

(
Scl(Σimi)

ΣiV (0.3mi)

)
, (2.5)

where V (m) is the hydrogen-ionising photon emission rate from a star with mass m, using the

fits from Vacca et al. (1996).

The factor 0.3 is an empirical factor to account for the scaling between the masses of sink

particles and those of massive stars, necessary because we do not fully resolve the fragmentation

of sink particles into proto-stars. See Section 2.3.1 for further discussion. The correction factor

X ≡ Scl(Σimi)/ΣiV (0.3mi) is very close to unity in most simulations in which we resolve

33



massive stars and it is introduced only to prevent overproducing ionising radiation in case massive

stars are poorly resolved. In all simulations we also impose X ≤ 1.

For the fiducial clouds we also include He0 and He+ ionising photons with total emission

rates being SHe0 = (1.178×1046 s−1)(M∗/M⊙) and SHe+ = (2.422×1043 s−1)(M∗/M⊙). These

rates are calculated using the same method as for hydrogen ionising photons described above us-

ing a Kroupa IMF (Kroupa, 2002). For the luminosity of individual stars we use Schaerer (2002)

fitting for Q(He0) and Q(He+) with extrapolations above 150 M⊙. This is contrasted by the

model of Gavagnin et al. (2017), who assume blackbody spectra for each star. See Appendix A.2

for details.

Various authors have concluded that UV photoionisation is typically the most important

process in regulating star formation on a cloud scale. Haworth et al. (2015) find that additional

processes beyond hydrogen photoionisation have a correcting factor of 10% at best. Radia-

tion pressure mainly becomes important at very high surface densities, which principally affects

smaller scales than the ones studied here - see Crocker et al. (2018) for idealized conditions

and Kim et al. (2018) for simulations with self-consistent star formation feedback. Dale et al.

(2014) further find that winds have a minimal effect on the star formation efficiency of molecular

clouds. We are thus justified in our choice to focus on UV photoionisation feedback in this work,

but discuss cases where this may not be sufficient later in the chapter.

2.2.4 Cooling

We use the radiative cooling function described in Geen et al. (2016). The cooling in

neutral gas is based on the prescription in Audit & Hennebelle (2005), which includes cooling

34



from carbon, oxygen and dust grains as well as the effect of the ambient UV background in

the ISM. For collisionally ionized gas at temperatures > 104 K we use Sutherland & Dopita

(1993) cooling function. The out-of-equilibrium cooling of photoionised hydrogen and helium

is treated as described in Rosdahl et al. (2013). Out-of-equilibrium cooling of photoionised

metals is treated with a piecewise fit to the cooling curve given in Ferland (2003). We assume a

uniform metallicity as listed in Table 2.1. For most simulations, this is solar metallicity, though

we perform some simulations at sub-solar metallicity. We do not implement out-of-equilibrium

molecular chemistry.

2.3 Results

In this section we present and discuss the results of our simulations. In § 2.3.1, we study

the mass function of cores (sink particles) and the IMF. In § 2.3.2 we focus on the star formation

efficiency and in § 2.3.3, on the star formation rate.

A representative sample of snapshots from our simulation set is shown in Figures 2.2 and

2.3. These figures show the time evolution of the density-weighted projections of the gas density

and temperature as well as the position of the sink particles along line of sight for three simula-

tions with different mean cloud densities (fiducial on the left column, compact, middle column,

and very-compact on the right column), and a cloud mass of 3.2× 104 M⊙.

We observe clearly that the star formation efficiency increases with increasing cloud den-

sity, and the stellar cluster that is formed at the end of the simulation remains more compact and

self-gravitating for the densest cloud. The effect of radiative feedback from massive stars is also

clearly visible in the density and temperature projections. H II regions break out of the dense

35



Figure 2.2: Line-of-sight projections of the (density-weighted) gas density for three simulations
with cloud mass 3.2× 104 M⊙. From left to right we show clouds with increasing mean density:
ngas ∼ 1.8 × 102 cm−3, representing our fiducial clouds in the local universe, ngas ∼ 1.8 ×
103 cm−3 and ngas ∼ 1.8 × 104 cm−3, respectively. From top to bottom we show the time
evolution of the clouds. Sink particles are displayed as cyan dots. The snapshots shown in the
top row represent the initial conditions of the turbulent cloud: no stars have formed at this time
because the highest density is below the threshold for star formation, but the cloud idealized
initial conditions have been already evolved for ∼ 3 tff in order to develop a turbulent density
field. In the bottom row snapshots, star formation has stopped and most of the gas has been
expelled as a result of UV feedback from massive stars.

36



Figure 2.3: Same as Fig. 2.2 but showing the density-weighted projection of the temperature.

37



101

102

103

104
XS-F

8.5 Myr, 124 M�
25 Myr, 381 M�

S-F
6.0 Myr, 77 M�
8.6 Myr, 419 M�
25 Myr, 508 M�

M-F
6.0 Myr, 240 M�
8.5 Myr, 905 M�
25 Myr, 1367 M�

L-F
6.0 Myr, 1072 M�
8.6 Myr, 3695 M�
21 Myr, 5696 M�

XL-F
6.0 Myr, 6330 M�
8.6 Myr, 13745 M�
22 Myr, 24565 M�

101

102

103

104

dN
/d

lo
g

M
(M
�

)

XS-C
8.0 Myr, 100 M�

S-C
1.8 Myr, 1 M�
2.6 Myr, 14 M�
8.3 Myr, 491 M�

M-C
1.8 Myr, 171 M�
2.6 Myr, 562 M�
8.3 Myr, 2976 M�

L-C
1.8 Myr, 1738 M�
2.6 Myr, 4717 M�
8.3 Myr, 13685 M�

10−1 100 101 102

101

102

103

104
XS-VC

0.6 Myr, 5 M�
0.9 Myr, 19 M�
2.5 Myr, 478 M�

10−1 100 101 102

S-VC
0.6 Myr, 131 M�
0.8 Myr, 358 M�
2.4 Myr, 3221 M�

10−1 100 101 102

M (M�)

M-VC
0.6 Myr, 1636 M�
0.9 Myr, 5680 M�
2.5 Myr, 14416 M�

10−1 100 101 102 103

L-VC
0.6 Myr, 9771 M�
0.9 Myr, 27352 M�

10−1 100 101 102 103

XXS-VC
0.9 Myr, 1 M�
2.5 Myr, 97 M�

0

2

4

6

t/
t ff

Figure 2.4: Temporal evolution of the mass functions of the stars, obtained by multiplying by 0.4
the masses of the sink particles (see text). The dashed lines are analytic Kroupa IMF for systems
normalised to the total mass of the sink particles at the corresponding time. The time and total
mass in sink particles are shown in the legend. We see good agreement between the shifted sink
particles mass function and the analytic mass-normalised Kroupa IMF at the high-mass end, both
in terms of the power-law slope and normalisation.

filaments destroying them and reducing the overall mass in dense gas in which stars are formed.

2.3.1 Stellar Initial Mass Function

2.3.1.1 Cores Fragmentation and Initial Mass Function

Maps in the continuum of cluster regions and larger areas in star-forming systems allow us

to construct a core mass function (CMF), that is the mass function of high-density gas concen-

trations (starless cores) with mass sufficiently large to be identified given the resolution of the

observations (e.g., Motte et al., 1998). Observations of a well-resolved CMF in the Pipe nebula

show a striking similarity to the stellar IMF, but shifted to higher masses by a factor of a few,

38



which suggests that the IMF is the direct product of the CMF with a roughly constant core-to-star

conversion efficiency ∼ 30% (e.g., Matzner & McKee, 2000; Alves et al., 2007)

Previous works on star formation in molecular clouds which adopted sink particles (like

in the present chapter) have investigated the mapping between the masses of pre-stellar cores

at the time they become self-gravitating and the final masses of the stars that form within them

(e.g., Padoan et al., 2001; Smith et al., 2009). For instance, Smith et al. (2009), using SPH

simulations, find that at early times the relationship between stellar masses and the parent cores

can be reproduced within a modest statistical dispersion with the star being about one-third of the

parent core mass.

We find results in agreement with these previous studies. Figure 2.4 (and Figure 2.7) show

the stellar mass function for our grid of simulations obtained assuming that the sink particles

are about a factor of 2.5 more massive than the corresponding massive stars they produce. This

means that the number of stars of a given mass (> 1 M⊙) is given by a Kroupa IMF for a star

cluster with total mass equal to the total mass of the sinks. In Figure 2.7 the IMF is shown at the

end of the simulation when star formation has stopped, while in Figure 2.4 we also show the time

evolution of the IMF, together with the Kroupa IMF for a cluster with total mass equal to the total

mass in sink particles (dashed lines). This means that we are assuming nearly 100% efficiency

of star formation in the cores (i.e., the cores fragment into stars) or a lower efficiency but the gas

expelled by feedback is later transformed into low-mass stars. However, we think that this second

model is less physically motivated. We can see that the shifted sink mass function (SMF) matches

the Kroupa IMF at the high-mass end, both in terms of the slope and normalisation at any time

during the formation of the star cluster. The figure suggests that the birth of stars in a cluster

follows the same random sampling of the universal mass function throughout the star-formation

39



process. In other words, it appears that there is not a bias toward formation of high mass-stars or

low-mass stars during the early times when the cluster is in the formation process.

As discussed above the sink particles can be interpreted as pre-stellar cores, and each sink

particle converts ∼ 40% of its mass to a single massive star, with the rest of the mass ending

up in low-mass stars, filling up the lower mass end of the IMF. The flattening or cut-off of the

IMF at the low-mass end observed in our simulations is likely due to insufficient spatial and mass

resolution to capture the formation of low mass cores or the fragmentation of more massive pre-

stellar cores. To further clarify, we note that the interpretation that only 40% of the core mass is

converted into a star and the remaining 60% is returned to the gas phase, never to participate in

star formation, would not produce the correct normalisation of the IMF. This is because in this

scenario, although the stellar masses are ∼ 40% of the core masses, the total stellar mass and the

star formation efficiency would be reduced by a factor of ∼ 2.5, lowering the expected number

of massive stars below the value found in the simulations.

So far we have assumed that all the gas in the cores fragments into stars with η = 100%

efficiency. In this case we find a conversion factor ε = 0.4 between the CMF and the IMF (such

that the normalisation of the IMF agrees with the observed one). However, it is possible to match

the observed IMF also in models in which η < 1. Note that in this case the TSFE shown in all our

plots should be re-scaled by a factor η. In models with η < 1, the conversion factor ε that matches

the mass-normalised empirical IMF, is ε = 0.4η1/Γ. For η = 1, 0.69 and 0.4, ε = 0.4, 0.3 and

0.2, respectively.

The results of this section justify our assumptions to model radiative feedback in Sec-

tion 2.2.3 and it further implies that our simulations are self-consistently treating the formation

of individual massive stars, their feedback effects, and can be used reliably to estimate of the

40



10−1 100 101 102

101

102

103
dN
/d

lo
g(

M
/
M
�

)