Magnetic Fields and Turbulence in Protoplanetary Disks
byDaniel. A. Gole
B.A., State University of New York, College at Geneseo, 2013
M.S., University of Colorado at Boulder, 2017

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillmentof the requirements for the degree ofDoctor of Philosophy
Department of Astrophysical and Planetary Sciences
2019This thesis entitled:
Magnetic Fields and Turbulence in Protoplanetary Disks
written by Daniel. A. Gole
has been approved for the Department of Astrophysical and Planetary Sciences

Prof. Philip J. Armitage

Dr. Jacob B. Simon

Prof. Mitchell C. Begelman

Prof. Steven R. Cranmer

Prof. Jeffrey B. Weiss

Date
The final copy of this thesis has been examined by the signatories, and we find that both thecontent and the form meet acceptable presentation standards of scholarly work in the above
mentioned discipline.

iii
Gole, Daniel. A. (Ph.D., Astrophysics)
Magnetic Fields and Turbulence in Protoplanetary Disks
Thesis directed by Prof. Philip J. Armitage

Protoplanetary disks (PPDs) are disks of gas and dust that form around young stars andevolve into planetary systems over the course of a few million years. PPDs are poorly ionized,with only the innermost radii being hot enough to be thermally ionized, and radiation only able
to ionize the surface layers of the optically thick disk. As a result, the disk will have turbulencein it’s surface layers driven by the magnetorotational instability (MRI), but the mid-plane will bedominated by non-ideal magnetohydrodynamic (MHD) effects that damp the MRI, leading to a“dead zone” (DZ). This thesis explores several topics relating to magnetic fields, non-ideal MHD
effects, and turbulence in PPDs.
We quantify how forced turbulence in Ohmic DZs varies with the ratio between the columnof magnetically active and inactive gas. We find that the Reynolds stress near the mid-plane drops
rapidly with increasing DZ thickness, becoming negligible for DZs with an active to dead mass
ratio less than a few percent. The DZ’s fluid motions are dominated by “r-modes”, a oscillatory
circulation pattern that is inefficient at transporting angular momentum.
We present a 1D model that couples the net-magnetic-flux and surface density evolution of
the disk. We use this model to investigate the interaction between the stellar magnetic cycle andthe disk, showing that some disks may exhibit cyclical or outburst-like accretion behaviour due tothe interaction between the net-flux in the disk and a Hall-effect dominated DZ.
We study turbulence in the outer regions of protoplanetary disks, where ambipolar diffusion
is the dominant non-ideal effect. We show that the properties of the MRI are not determined purely
locally, with transport of toroidal magnetic flux from the active layers to the mid-plane playing asignificant role.
We investigate the streaming instability’s ability to form planetesimals in a turbulent en-

ivvironment. We find that turbulence does not greatly change the shape of the mass-spectrum ofplanetesimals that form. However, no planetesimals are able to form at all in the presence of turbulence with α ∼ 10−3 , implying that DZs are essential if planetesimals are to be formed by thismechanism.

v

Acknowledgements

I would like to thank the following people and groups of people that contributed to this thesis,
both directly and indirectly:• My advisors for this thesis: Phil Armitage and Jake Simon. I appreciate your patience withme over these 5 years, and your willingness to allow me to take projects in the directions
that I was most excited about, while still gently steering me in the right direction.
• The members of my thesis committee at-large: Mitch Begelman, Steve Cranmer, and Jeff
Weiss. I know this is a time consuming task for you all and I greatly appreciate yourwillingness to participate in the process.
• Steve Lubow, Xeuning Bai, and Rixin Li for useful conversations, help with codes, andfeedback on my work.
• My research and teaching mentors during my time as an undergraduate: Aaron Steinhauer,David Meisel, and Charlie Freeman at SUNY Geneseo and Heidi Newburg at Rensselaer
Polytechnic Institute.
• My parents, Tom and Lisa, who were my first science teachers and have encouraged andsupported me through 23 years of school. I would also like to thank my brother and lifelong
best friend, Mike.
• The graduate students of the APS department, in particular those in my cohort who helped
make my 6 years in Boulder very enjoyable: Nick, Addi, Becky, Hillary, Ryan, Daniel, Erika,

viMarcus and Drake. I feel very lucky to have ended up at the same graduate program as allof you.
• My band-mates: Marek, Tristan, Addi, and Bethany for allowing my to keep music a partof my life all these years. I would also like to thank my many ultimate frisbee teammatesduring my time here for providing a welcoming and supportive community in which to both
improve and have fun.
I acknowledge support from NASA under awards NNX13AI58G, NNX14AB42G and
1280NSSC18K0640 from the Origins of Solar Systems and Astrophysics Theory programs. I also acknowledge support from NASA through grants NNX13AI58G, NNX14AB42G and NNX16AB42G(P.J.A), from the NSF through grant AST 1313021 (P.J.A.), and from grant HST-AR-12814 (P.J.A.)
awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. J.B.S.’s support wasprovided in part under contract with the California Institute of Technology (Caltech) and the Jet
Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executedby the NASA Exoplanet Science Institute.
This work utilized the Janus supercomputer, which is supported by the National Science
Foundation (award number CNS-0821794), the University of Colorado Boulder, the University ofColorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer isoperated by the University of Colorado Boulder. S.H.L. acknowledges support from NASA grant
NNX11AK61G. This work also utilized the RMACC Summit supercomputer, which is supportedby the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University ofColorado Boulder, and Colorado State University. The Summit supercomputer is a joint effortof the University of Colorado Boulder and Colorado State University. Finally, I acknowledge theTexas Advanced Computing Center (TACC) at The University of Texas at Austin for providing
HPC resources that have contributed to the research results reported within this thesis. URL:
http://www.tacc.utexas.edu.

Contents

Chapter
1 Introduction
1.1

1.2

1.3

1.4

1

Protoplanetary Disk Physics and Phenomena . . . . . . . . . . . . . . . . . . . . . .

1

1.1.1

Formation and Basic Structure . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1.2

Accretion, Angular Momentum Transport, and Turbulence . . . . . . . . . .

2

1.1.3

The Magneto-Rotational-Instability . . . . . . . . . . . . . . . . . . . . . . .

3

1.1.4

The Ionization State of PPDs . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.1.5

Non-Ideal MHD Effects and Dead Zones in PPDs . . . . . . . . . . . . . . . .

6

1.1.6

Observations of PPD Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.1.7

Winds from PPDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.1.8

Episodic Accretion in PPDs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Solids in Disks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.1

Vertical Settling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.2.2

Radial Drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

An Overview of Planet Formation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1.3.1

Coagulation at Small Size-Scales . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.3.2

Planetesimal Formation via the Streaming Instability . . . . . . . . . . . . . 19

1.3.3

Planetesimals to Terrestrial Planets . . . . . . . . . . . . . . . . . . . . . . . 20

1.3.4

Turbulence, Dead Zones, and Planet Formation . . . . . . . . . . . . . . . . . 21

Organization of This Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

viii
2 Turbulence, Transport and Circulation in Ohmic Dead Zones

23

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.2

Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.3

2.2.1

Numerical Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.2.2

Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1

General, Large-scale Properties . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.3.2

Saturation Time-scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.3.3

Stresses and energies in the saturated state dead zone . . . . . . . . . . . . . 37

2.3.4

Transport Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

2.3.5

Density Perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.3.6

Accretion Rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.3.7

Structure of the Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

2.3.8

Large scale velocity structures in the dead zone . . . . . . . . . . . . . . . . . 47

2.4

Implications for Particle Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

2.5

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3 The Interaction of the Stellar Magnetic Cycle with Hall Dead Zones in Protoplanetary Disks 60
3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.2

The Model

3.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

3.2.1

Static Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

3.2.2

Coupling MRI Turbulence to the Net Magnetic Field

3.2.3

Surface Density and Net Magnetic Field Evolution . . . . . . . . . . . . . . . 66

3.2.4

Evolving the Model Self-Consistently . . . . . . . . . . . . . . . . . . . . . . . 67

3.2.5

Summary of Disk Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.2.6

Stellar Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

. . . . . . . . . . . . . 64

Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

ix
3.3.1

Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3.3.2

Boundary Conditions and Truncation Radius . . . . . . . . . . . . . . . . . . 70

3.3.3

Smoothing α . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

3.3.4

Advancing in Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

3.4

Modeling the Influence of the Stellar Magnetic Cycle on the disk . . . . . . . . . . . 72

3.5

A Criterion for Outbursts from Hall Dead Zones . . . . . . . . . . . . . . . . . . . . 75
3.5.1

Outburst Mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

3.5.2

Outburst Criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

3.6

Results of Time-Dependent Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

3.7

Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4 The Nature of Turbulence in the Outer Regions of Protoplanetary Disks

87

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.2

Numerical Algorithm, Parameters, and Initial Conditions . . . . . . . . . . . . . . . 90

4.3

4.4

4.2.1

Resistivity Profile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

4.2.2

Resolution Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

4.2.3

List of Simulations and Varied Parameters . . . . . . . . . . . . . . . . . . . . 96

Stratified Ambipolar Simulations with no FUV . . . . . . . . . . . . . . . . . . . . . 96
4.3.1

Expectations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

4.3.2

Magnetic Field Strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

4.3.3

Turbulence and Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

Stratified Ambipolar Simulations with FUV . . . . . . . . . . . . . . . . . . . . . . . 105
4.4.1

Mid-plane field strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

4.4.2

Hydrodynamic Turbulent Properties . . . . . . . . . . . . . . . . . . . . . . . 107

4.4.3

Stress and Magnetic field perturbations . . . . . . . . . . . . . . . . . . . . . 107

4.5

Vertical transport of magnetic field in the FUV simulations . . . . . . . . . . . . . . 111

4.6

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

x
5 Forming Planetesimals in a Disk with Driven Turbulence

118

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

5.2

Methods, Setup, and Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

5.3

5.4

5.5

5.6

5.2.1

Fluid Algorithm and Parameters . . . . . . . . . . . . . . . . . . . . . . . . . 120

5.2.2

Particle Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

5.2.3

Particle Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

5.2.4

Stirring Turbulence: Spectrum of Velocity Perturbations . . . . . . . . . . . . 125

5.2.5

Stirring Turbulence: Magnitude of Velocity Perturbations . . . . . . . . . . . 126

5.2.6

Characteristics of Driven Turbulence . . . . . . . . . . . . . . . . . . . . . . . 129

5.2.7

Varied Parameters and Table of Simulations Performed . . . . . . . . . . . . 133

Analysis Techniques for Clumps and Particles . . . . . . . . . . . . . . . . . . . . . . 133
5.3.1

Finding Gravitationally-Bound Particle Clumps . . . . . . . . . . . . . . . . . 133

5.3.2

Statistics of the Planetesimal Mass Spectrum . . . . . . . . . . . . . . . . . . 135

5.3.3

When to Measure the Clumps . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

Results I: Snapshots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
5.4.1

The Turbulence Threshold for Planetesimal Formation . . . . . . . . . . . . . 142

5.4.2

Mass Fraction in Planetesimals . . . . . . . . . . . . . . . . . . . . . . . . . . 143

5.4.3

Differential Mass Spectra and Simple Power Law Fits . . . . . . . . . . . . . 143

5.4.4

More Sophisticated Fits to the Cumulative-Mass-Spectrum . . . . . . . . . . 146

Results II: Clump Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
5.5.1

Methods and Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147

5.5.2

Results from Clump Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . 151

Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

6 Summary and Looking Forward

158

6.1

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158

6.2

Looking Forward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

xi
Bibliography

162

xii

TablesTable
2.1

Dead Zone Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.2

Summary of time and volume averaged values in the steady state. Times are given
in units of orbits. In the rightmost two columns, ”T.E.” refers to the transport
efficiency. We compare the transport efficiency in the dead and active zones, as
defined in equation 2.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

2.3

Parameters of the circulation modes, vertical particle diffusion parameters, and critical τs values for each run.

3.1

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Parameter choices for our disk model. Given these parameters, a specification of M∗(which determines Ω), and ṁ, our static disk model determines all other quantities
as a function of radius. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

3.2

Assumed stellar parameters for several stellar types. . . . . . . . . . . . . . . . . . . 69

4.1

Simulation Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.1

Simulation Parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

5.2High level summary of clumping behaviour in high resolution simulations. tfirst
denotes the time of the formation of the first clump. tsnap shows the time at which
we measure the values for the next 3 columns. . . . . . . . . . . . . . . . . . . . . . . 145

xiii
5.3

Summary of fits to the cumulative mass spectrum of planetesimals based on single
snapshots from our high resolution simulations. ∆IC denotes the average of the two
information criteria considered. We give the parameters for the STPL and BCPL
for all runs regardless of ∆IC as a way of comparing from run-to-run and because
the interpretation of their parameters is intuitive.

5.4. . . . . . . . . . . . . . . . . . . 148

Summary of fits to the cumulative mass spectrum of planetesimals based on the initial
masses of each clump from our clump-tracking analysis. We give the parameters for
the STPL and BCPL for all runs regardless of ∆IC as a way of comparing from runto-run and because the interpretation of their parameters is intuitive. For all cases
here, no fits other than the BPL had ∆IC < 10 For reference, the SPL fits to thedifferential spectra for this method are 1.372 ± 0.023, 1.357 ± 0.035, and 1.380 ± 0.047
for the 3 runs shown from top to bottom, respectively.

. . . . . . . . . . . . . . . . 155

xiv

FiguresFigure1.1

A cartoon illustration of the linear regime of the ideal MRI. The two dots represent fluid elements and the “spring” represents a magnetic tension force. Due todifferential rotation in the disk, the fluid elements have different velocities, causing
them to drift apart over time and stretch the spring. This causes a prograde force
on the outer element which will gain angular momentum and drift outward, and a
retrograde force on the inner element which will loose angular momentum and drift
inward. This drift stretches the spring further and reinforces the instability. . . . . .

1.2

5A cartoon “map of the disk”, emphasizing the most relevant non-ideal MHD effects
in each region. We also point out the various regions and phenomena studied byeach chapter of this thesis. Credit to Armitage (2007), and Simon et al. (2018) for
the general idea for figures of this type. A few purely notational choices that requireclarification: 1. There are not actually breaks in the particle layer as shown, this isjust to allow room for text. 2. Planetesimal formation is not necessarily specific to
the radii shown here, it will generally happen throughout the radial extent of the disk.

2.1

8

Resistivity η as a function of z for each of our simulations, as labeled in the legend.
The mass in the active zone as a fraction of the total mass varies between about 0.58
(for simulation “S”) and 0.029 (for simulation “L”). . . . . . . . . . . . . . . . . . . 30

xv
2.2

The dependence of the Elsasser number, Λ on z for each of our simulations, as labeled
via the legend. The solid black horizontal line shows where Λ = 1. The dead zone isdefined here as the region where Λ < 1. . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.3

The log of the velocity magnitude for the medium sized dead zone run, plotted withdarker colors indicating higher velocity. The outer regions are clearly turbulent,
transitioning into more laminar flow at the mid-plane. . . . . . . . . . . . . . . . . . 34

2.4

Space-time plots of the Maxwell stress = −Bx By (first and third plots) and Reynolds
stress = ρvx δvy (2nd and 4th plots) for the small (top 2) and medium (bottom 2)
dead zone sizes. The quantities are averaged in the x and y directions. Note thatthe color scales are kept constant for like-quantities across the runs. . . . . . . . . . 35

2.5Time evolution of the volume averaged Kinetic Energy, Reynolds Stress, and Maxwell
Stress (from top to bottom), for the medium sized dead zone. The average over thedead zone is taken from |z/H| = 0 to 0.5, and for the active zone from |z| = 2.0 to
3.3. In addition, the plots were smoothed by time averaging over 10 orbits. . . . . . 36

2.6

The saturation time as a function of the ratio of the dead zone mass to the active
zone mass.

2.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

Time averaged vertical profiles of the Reynolds stress (left) and Maxwell stress (right)
for the different sized dead zones. The sharp transition seen near |z| = 3 in the
Reynolds stress for the ML and L runs corresponds to where the r modes (discussed
in detail in §2.3.8) shut off. The r modes have a larger spatial extent in the ML and
L runs due to the larger dead zone. Therefore the location of the transition feature
happens in a region of lower Reynolds stress, making the features most prominent
in these runs.

2.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

The average Reynolds stress near the mid-plane as a function of the ratio of theactive zone mass to the dead zone mass. Note that the mass ratio on this x-axis is
inverted relative to figures 6 and 11 to show the linear relationship.

2.9. . . . . . . . . 40

The density perturbation profiles for the 5 runs, averaged over the saturated state. . 43

xvi
2.10 The non-dimensional accretion rate as a function of z for the different dead zone
sizes, averaged over the saturated state. . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.11 The ratio of the accretion rate in the dead zone to the accretion rate in the active
zone as a function of the mass ratio of the two zones. The accretion rate was averaged
over the region defined to be the dead zone in table 2.1, which changes from run torun. Time averaging was done over the saturated state. . . . . . . . . . . . . . . . . 44
2.12 Autocorrelation functions for the density perturbation for the medium dead zone(top) and large dead zone (bottom). Moving from left to right, the ACF is calculatedusing planes |z| = 0 to 0.5, 1.0 to 1.5, 2.0 to 2.5 and 3.0 to 3.5. Each plot is 2H on aside. The plots on the left show that perturbations near the mid-plane get shearedout by the background flow, while perturbations in the upper layer are smaller andless aligned with the y-axis.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

2.13 Autocorrelation functions for the vertical velocity for the medium dead zone (top)and large dead zone (bottom). Moving from left to right, the ACF is calculated using
planes |z| = 0 to 0.5, 1.0 to 1.5, 2.0 to 2.5 and 3.0 to 3.5. Each plot is 2H on a side.
The large scale circulation modes can clearly be seen in the leftmost plots, while theplots on the right corresponding to the upper layers are consistent with the structureof MRI-driven turbulence in the ideal MHD limit, as seen in other quantities andprevious studies (e.g., Guan et al., 2009). . . . . . . . . . . . . . . . . . . . . . . . . 46
2.14 The angle of the density perturbation ACF away from the y-axis for each of the deadzone runs as a function of z. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.15 A snapshot from the medium-sized dead zone run of the vertical velocity in the xz plane, averaged over y. The figure clearly depicts the presence of two verticallycoherent updrafts / downdrafts, stretching throughout the dead zone region.

. . . . 50

2.16 The vertical velocity plotted for the x-y plane at the mid-plane of the medium-sizeddead zone run. Time advances from top to bottom by 3.5 Ω−1 each panel, showing
half of a cycle of this standing wave pattern. . . . . . . . . . . . . . . . . . . . . . . . 51

xvii
2.17 The average auto-correlation function of the random vertical velocity for all points−0.5 < z/H < 0.5 as a function of the lag time. The vertical diffusion parameter isthe area under this curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.1

A qualitative representation of the static magnetic field geometry. The black circle
in the center represents the star and the gray shaded area is the disk. The disk isthreaded by open field lines shown in red. The star produces closed dipole-like fieldlines, shown in blue. In isolation this would be a true dipole, but the interaction
with the open line region causes the dipole lines to be compressed inward, while theopen lines are pushed outward. This results in a region where there are concentrated
closed field lines from the star overlapping the inner edge of the disk.

. . . . . . . . 74

3.2

A cartoon illustration of the outburst mechanism proposed by Armitage (2016).

. . 77

3.3

The power laws for the diffusion time at the radius where the Hall DZ begins plottedfor several stellar masses and two α values. The α values here refer to α in the active
region. The The horizontal line shows the solar cycle period for the sun, roughly 20
years. For any given ṁ and M∗ , if tcycle is above or roughly equal to the diffusion
time, outbursts or cyclic behavior would be expected. . . . . . . . . . . . . . . . . . 80

3.4

A space-time diagram of the net-field in the disk for a G star with tcycle = 20yr andṁ = 10−8 M /yr. The horizontal dashed line represents the location of the inneredge of the Hall DZ. In this case, the diffusion fronts of aligned flux do not reachthe inner edge of the dead zone before the stellar field changes sign, so no episodic
or cyclic accretion occurs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

3.5

Space-time diagrams and the disk luminosity for a M4 star with tcycle = 20yr andṁ = 10−10 M /yr. The top 3 panels show net-field, local accretion rate, and surface
density respectively, with the bottom plot showing the luminosity calculated via
equation 3.35.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

xviii
3.6

Light curves for all 5 stellar types considered for an accretion rate of ṁ = 10−9 M /yr.
The stellar types decline in mass from top to bottom. The stellar cycle time increases
from left to right: 6, 20, and 40 years respectively. Note that the scale of the x-axis
is not the same in all plots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

4.1

The critical and most unstable MRI wavelengths at the mid-plane and at z/H = 1.5)
in terms of our fiducial cell size δx = H/30. β0 = 104 at the mid-plane, whichcorresponds to β ∼ 103 at z/H = 1.5 due to the decrease in gas pressure. . . . . . . 97

4.2

Temporally and spatially averaged β within 0.5H of the mid-plane. The black line
shows βmin (see Bai and Stone, 2011). βy generally follows βmin . . . . . . . . . . . . 101

4.3

β distributions for simulations with Am = 0.1, 1, 10 in blue, green, and red respectively. The dashed vertical lines corresponds to βmin for each simulation. In general,βmin gives a reasonably good description of the typical field strength.

4.4. . . . . . . . 101

Hydrodynamic diagnostics of turbulence δρ and δv time and space averaged (within
0.5H of the mid-plane) for simulations SLN Am0.032 through SLN Am10. The levelof turbulence increases as the damping decreases, consistent with expectations. . . . 103

4.5

The velocity perturbation profiles for our series of simulations with no FUV layer(constant Am). The profiles are shown in a gradient from blue to red, corresponding
to Am = 0.032, 0.1, 0.32, 1.0, 3.2, and 10.0. The observed limits on TW Hya andHD163296 (Flaherty et al., 2015, 2017, 2018) are shown in the dotted and dashedlines respectively. (Note that in the case of TW Hya, this limit only applies to|z| ∼ 2–3H.) Even in a disk with no FUV-ionized layers, for β0 = 104 it would takea diffusivity of at least Am = 0.032, and likely higher, to explain the observations of
weak turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

xix
4.6

The magnetic stress calculated from perturbations (black) and from laminar fields(blue). Note that at Am = 1, the points are directly on top of one another. The
stress clearly scales with Am. For Am < 1, the stress is dominated by small-scale
perturbations.

4.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

The time-averaged β parameter plotted versus Am for SLF Am0.01 S0.1 throughSLF Am100 S0.1 for 3 spatial regions: the mid-plane, above the transition, andthe corona. The black line shows the minimum β for the mid-plane MRI (given inequation 4.10). Because the field strength exceeds what would be expected from thelocal MRI for Am < 1, we expect that this field is sourced elsewhere and transported
to the mid-plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

4.8

Distribution of β taken over the mid-plane region of SLF Am0.032 through SLF Am1.0,
with the mid-plane Am shown on each plot. The vertical dashed lines indicate theminimum beta that permits the MRI for each Am parameter. While the MRI mightoccasionally be active in small patches, it is unlikely to dominate the dynamics ofthe mid-plane when Am < 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

4.9

Velocity and density perturbations for SLF Am0.01 through SLF Am100. The midplane quantities do not scale with Am for Am < 1 indicating that these perturbations
are not generated by local MRI. In addition, the mid-plane turbulence does not seemto affect the level of turbulence in the corona. . . . . . . . . . . . . . . . . . . . . . . 108

4.10 Magnetic field perturbations versus Am within different vertical regions. The midplane perturbations scale roughly as Am1/3 , indicated by the black dashed line.
Despite the mid-plane scaling, the mid-plane does not influence the other regions.

. 110

xx
4.11 The magnetic stress calculated from perturbations (black) and from laminar fields(blue). The red line indicates the predicted upper limit of stress from Bai and Stone(2011). The black dashed line indicates the expected power law scaling of the stressfrom perturbations: Am2/3 for Am < 1 (based on the scaling seen in figure 4.10),
and Am1.2 for 1 ∼ Am (based on equations 4.10 and 4.13). The mid-plane stress isdominated by the laminar component for very diffusive mid-planes (Am < 1). . . . . 110
4.12 Space-time diagram of sign(By )By2 /2P0 for SLF Am0.1. Here, P0 is the mid-planepressure. The horizontal dashed lines indicate zFUV . . . . . . . . . . . . . . . . . . . 113
4.13 Time shift of maximum correlation relative to the mid-plane as a function of z for the
SLF Am0.1 simulation. The dashed vertical lines indicate the center of the transition
between FUV ionization and no FUV ionization. We see that the mid-plane lags
behind the transition region, implying that diffusion towards the mid-plane is able
to beat out the effect of buoyancy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.1

The relationship between the driving amplitude dE/dt and the resulting averagevelocity perturbation δv. Each black point is generated by averaging over a lowresolution simulation. The grey lines show an linear interpolation between the black
points in log-space. The red lines demonstrate how the appropriate driving rate isfound to produce a simulation with δv = 10−2 , ie. α = 10−4 .

5.2. . . . . . . . . . . . . 128

A demonstration that our stirring algorithm produces a box of homogeneous, isotropic
turbulence. Data is taken from simulations hr moderate np. Time averages are takenfrom 5Ω−1 to 10Ω−1 . The dashed horizontal line shows the expected δv given thestrength of the driving.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

xxi
5.3

A demonstration that our stirring algorithm produces the roughly same spectrum
for 3 different levels of driving, and that the spectrum roughly follows a p = −5/3
power law over the inertial range. The power plotted here is normalized such thatthe mean is zero. The power is also compensated with a k ∝ −5/3 dependence,such that a −5/3 power law would appear perfectly horizontal on this plot. Timeaverages are taken from 10Ω−1 to 50Ω−1 . . . . . . . . . . . . . . . . . . . . . . . . . 132

5.4

Mid-plane slices of the particle density at the onset of the streaming instability in asimulation with no driven turbulence (hr control). . . . . . . . . . . . . . . . . . . . 136

5.5

Typical differential and cumulative mass-histograms of the planetesimals. Datashown is from hr control at t = 27Ω. The black dashed lines show a simply powerlaw fit to the data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

5.6

An analysis of our control run: hr control, which does not have any driven turbulence.
The top 3 plots from top to bottom show: the time evolution of the slope of thesimple power law fit to the differential histogram (with a 1σ error bar), the massfraction of planetesimals, and the number of gravitationally bound clumps. The
horizontal dashed lines roughly show the range of power law indices determined byprevious studies. The dashed blue line shows the time at which we choose to measure
the spectrum, with that spectrum shown at the bottom in blue. Similarly, the redline shows where the the spectrum shown in red is taken. The oval denotes theintermittent population.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

5.7

An analogous figure to the previous one for our fiducial run: hr moderate.

. . . . . 145

5.8

Cumulative spectra and models for snapshots of hr control, hr weak, and hr moderatefrom top to bottom. Fit parameters to the best models can be found in 5.3. The
square points show the characteristic mass scale relevant to each fit. The bottom
histogram shows the true values, with each plot above offset by a factor of 10 for
convenient visualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

xxii
5.9

An illustration of our particle tracking scheme. The example shown here is intended
to be demonstrative of our scheme and does not represent actual data from a simulation. Each circle represents a gravitationally-bound clump. Red corresponds to
a new clump, gray to a previously found clump that was identified again, white toa clump that was not found (but will continue to be looked for in future frames),
and blue to a new clump that formed from a split. Clumps are consistently labeledfrom frame to frame. Clumps that merge are given 2 labels and are continued to be
tracked as both clumps. In taking statistics, we use the masses of all of the red clumps.152

5.10 An analogous figure to figure 5.8, again for hr control, hr weak, and hr moderatefrom top to bottom. In this figure, only the initial mass of each clump is used, as determined by our clump tracking algorithm. The square points show the characteristic
mass scale relevant to each fit.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Chapter 1

Introduction

1.1

Protoplanetary Disk Physics and Phenomena
Accretion disks are ubiquitous throughout the universe, found rotating about a wide array ofdifferent objects including super massive black holes (AGN), white dwarves (dwarf novae), neutronstars/stellar-mass black holes (x-ray binaries), and young stars (circumstellar or protoplanetarydisks (PPDs)). While these disks exist across many orders of magnitude in size-scale and energy,
they all share a very general common origin: they form as the result of gravitational collapse of asystem that has net-angular-momentum. This net-angular-momentum defines an axis about which
the resulting system will rotate as the matter in the disk structure accretes inwards towards thecentral object.

1.1.1

Formation and Basic Structure

PPDs form as a part of the star formation process. When a large cloud of gas becomes
sufficiently cool to no longer provide enough pressure support to balance it’s self-gravity, it collapses. The resulting system is called a young stellar object (YSO). After an initial phase during
which the YSO accretes it’s surrounding envelope, the system arrives at the configuration we willconsider throughout this thesis: an (optically visible) star surrounded by a disk of gas and dust.
Observationally, these objects are referred to as a Class II source (Lada and Wilking, 1984), or a“classical T Tauri Star” (Joy, 1945).
Disks extend from ∼ 0.1 to ∼ 100 AU in the radial direction and are very thin in the vertical2direction, typically having aspect ratios (h/r) on the order of a few percent. Disk masses are on theorder of 1 to 10% of the stellar mass. To first order, the structure of the disk is dominated by thestar’s gravity and the disk’s self-gravity can be ignored. As a result, disks generally rotate withina few percent of the Keplerian orbital velocity (more on this in §1.2). In the vertical direction,
the disk’s pressure profile balances the vertical component of the star’s gravity. For a thin disk
(which is always the case for a PPD), the pressure profile in the vertical direction works out to be a
Gaussian with scale height h = cs /Ω, where cs is the sound speed in the fluid and Ω is the angularorbital velocity.

1.1.2

Accretion, Angular Momentum Transport, and Turbulence

The most fundamental question in protoplanetary disk theory is quite simple: how do they
accrete? Based on observations of the disk-fractions of stellar populations in star forming regions(eg. Haisch et al. (2001)), the observed lifetime of a typical disk is on the order of a few million
years. This gives a time-scale during which the disk must accrete onto the central star. The specific√angular momentum of a fluid element in an accretion disk is given by l = GM∗ r, which is notably
an increasing function of radius. So in order for accretion to occur, angular momentum must betransported outward to compensate for the inward moving mass. A logical source of this angularmomentum transport is viscosity, which allows adjacent annuli within the disk to transfer energy
and angular momentum to one another. Outer annuli will rotate more slowly and tug backward on
inner annuli, causing the inner annuli to lose angular momentum (or rather, transport it outwards).
So while viscosity allows angular momentum to be transported, the source of this viscosity
in an accretion disk is not obvious. The molecular viscosity for a typical protoplanetary disk
is ∼ 2.5 ∗ 107 cm2 s−1 . Taking a typical length scale to be 10 AU yields a viscous time-scale of∼ 3 × 1013 years. This is over 2000 times the life of the universe, never-mind the few million year
observed evolution time of a typical protoplanetary disk. The molecular viscosity is far too small to
explain the observed lifetime, so an effective viscosity must arise from local turbulence (Balbus andPapaloizou, 1999). Turbulence then generates hydrodynamic and magnetic stresses that transport

3
angular momentum. Considering turbulent eddies on scale h (the pressure scale height of thedisk) circulating at the sound speed cs , the turbulent viscosity is conventionally parameterized asνturb = αcs h, where α is an unknown dimensionless parameter that scales the strength of turbulence(Shakura and Sunyaev, 1973). Thus, if an α value is specified then the evolution of the surface
density of the disk can be determined. While this is a convenient convention, we still have notspecified the physical source of this turbulence or the value of α.

1.1.3

The Magneto-Rotational-Instability

The leading candidate for driving turbulence in accretion disks is the magnetorotational
instability (MRI) (Balbus and Hawley, 1991, 1998). This instability arises from a magnetic field
coupling to deferentially rotating fluid. As the fluid shears, field lines must stay frozen into the fluid,causing the field to shear and stretch. Fluid parcels at larger radii feel a force in the retrograde
direction, causing them to move farther outward. In turn, the field stretches more and exerts astronger force. A cartoon illustration of this instability is shown in figure 1.1.
Under ideal magnetohydrodynamics (MHD), a disk with a vertical magnetic field and differential rotation has a local, linear instability if
Bz2
3ρc2 d≤ 2 s , (Ω2 ) < 0,
8π
π
dr

(1.1)

where Bz is the vertical magnetic field, ρ is the density, cs is the sound speed, r is the radial
coordinate, and Ω is the rotation rate. Note that this is a weak-field instability, making this criteria
fairly trivially satisfied. Under these conditions, the linear growth rate of the fastest growing MRI
mode is ∼ Ω. This growth rate is very fast, so the instability will quickly saturate and enter the
non-linear regime.
In the case of a vertical seed-field, the wavelength of the most unstable mode of the MRI−1/2

is given by λm /h = 9.18β0, where β0 is defined to be the standard plasma β parameter: theratio of the gas pressure to the magnetic pressure (Hawley et al., 1995). For a reasonable choiceof β0 = 104 , this characteristic MRI wavelength is about one tenth of the scale height. Hence,

4
if one wants to resolve the MRI’s behaviour and study the non-linear regime numerically, a local
simulation is usually necessary.
This introduces the “shearing box” approximation: a small patch of the disk is simulated in
the co-rotating frame (Hawley et al., 1995). In this context, Cartesian coordinates can be assumed,
radial gradients can be ignored, and a shearing background velocity must be applied such that the
inner edge of the box has a positive azimuthal velocity and the outer edge of the box has a negative
azimuthal velocity. Early numerical studies showed that the nonlinear state of the MRI sustains
turbulence with an effective α parameter of . 10−2 (Brandenburg et al., 1995; Stone et al., 1996),
with most angular momentum transport occurring due to magnetic stresses. This result was very
exciting for theorists, as this α parameter predicts viscous evolution time-scales on the order ofmillions of years, matching the expectation from observational constraints on disk lifetimes. Morerecent studies have largely confirmed this level of turbulence in a sufficiently ionized fluid, and havealso showed that the MRI can act even without a net-vertical magnetic field (Simon et al., 2012;
Davis et al., 2010).
While the ideal MRI described here may explain angular momentum transport in high energy
disks like those around black holes, we have not yet addressed the fact that PPDs are generally
not well ionized. In the next few sections, we consider the ionization state of the disk and the roleof non-ideal MHD effects in determining the the behaviour of the MRI and angular momentum
transport.

1.1.4

The Ionization State of PPDs

The degree of coupling between the magnetic field and the gas is crucial to the efficacy of
the MRI in transporting angular momentum. In fully ionized disks, such as those around black
holes, the coupling is very strong and the MRI will drive turbulence at a level of α = 10−2 , or evengreater in the presence of a strong net-vertical field (Simon et al., 2013a; Salvesen et al., 2016).
However, the colder environments present in protoplanetary disks lead to low ionization levels; theresult is substantially weakened gas-field coupling such that non-ideal MHD effects are relevant,5

Figure 1.1 A cartoon illustration of the linear regime of the ideal MRI. The two dots represent fluid
elements and the “spring” represents a magnetic tension force. Due to differential rotation in thedisk, the fluid elements have different velocities, causing them to drift apart over time and stretch
the spring. This causes a prograde force on the outer element which will gain angular momentum
and drift outward, and a retrograde force on the inner element which will loose angular momentum
and drift inward. This drift stretches the spring further and reinforces the instability.

6
and the MRI cannot operate at full capacity throughout the disk.
First, we will consider thermal ionization. Typical temperatures at the mid-plane of thedisk vary from up to a few thousand K at the innermost radii (inside ∼ 1AU), to tens of K in
the very outer regions (∼ 100AU). Even at the innermost radii, this temperature is less than thethermal ionization temperature for Hydrogen (∼ 104 K), and even the most easily ionized trace
species like alkali metals (∼ 5000K). However, it turns out that the disk does not need to beparticularly well ionized in order for the field to be coupled to the gas, with the critical ionization
fraction lying around 10−13 (Gammie, 1996). Based on chemical equilibrium modeling done by
Umebayashi (1983), only a relatively modest temperature of 800-1000 K is required to achieve thiscritical ionization fraction. As a result, we expect that the inner radii of PPDs will be sufficiently
thermally ionized to be in the ideal MHD regime.
In the outer disk, thermal ionization can largely be neglected, however radiative ionizationfrom cosmic-rays (Gammie, 1996), X-rays (Igea and Glassgold, 1999), and FUV photons (PerezBecker and Chiang, 2011a) will play a role. Gammie (1996) put forth a layered accretion model inwhich the surface layers of the disk are well-ionized by radiation and exhibit ideal MRI behaviour,while the mid-plane of the disk is shielded from radiation and is very weakly ionized. This leads
to a large magnetic diffusivity in the mid-plane region that quenches the MRI there, creating aso-called “dead zone”. While the layered accretion model provides a great starting point, there are
3 distinct non-ideal MHD effects that will be relevant in different regions of the disk, and each has
a different effect on the MRI. We will discuss these effects in the following sections.

1.1.5

Non-Ideal MHD Effects and Dead Zones in PPDs

Within a PPD, the gas density, ionization fraction, and field strengths all vary by ordersof magnitude. As a consequence, we have no expectation of universality in the properties of disk
turbulence sourced from MHD processes like the MRI, whose characteristics will vary with theimportance of the 3 non-ideal MHD terms: Ohmic diffusion, ambipolar diffusion, and the Hall
effect (Armitage, 2011; Turner et al., 2014; Balbus and Terquem, 2001; Kunz and Balbus, 2004; Xu

7and Bai, 2016). These 3 non-ideal effects of MHD break the ”frozen in” condition of the magnetic
field that is true for ideal MHD. The non-ideal induction equation is given by∂B
J×B J×B×B+= ∇ × v × B − η∇B −∂tenecγρi ρ

(1.2)

where B is the magnetic field, v is the velocity, J is the current, ρ is the total density, ρi is the ion
density, η is the resistivity, ne is the electron number density, and γ is the drag coefficient, related
to the ion-neutral collision rate. The 2nd, 3rd, and 4th terms on the right hand side correspond toohmic diffusion, the Hall effect, and ambipolar diffusion respectively.
While the relative importance of these 3 terms is itself an interesting and relevant topic
of study (Lesur et al., 2014; Kunz and Balbus, 2004), as a general rule ambipolar diffusion will
dominate at low densities, ohmic diffusion will dominate at high densities, with the Hall effectdominating in an intermediate regime, but also being relevant at high densities. Translating toregions of the disk, large radii (∼ 30 − 100 AU) will be dominated by ambipolar diffusion, (e.g.,
Desch, 2004; Kunz and Balbus, 2004; Simon et al., 2013b,a), while at AU-scales both Ohmic
diffusion and the Hall effect modify the mid-plane dynamics (Lesur et al., 2014; Bai, 2014; Simon
et al., 2015a), with ambipolar diffusion playing an important role in the atmosphere (Desch, 2004;
Kunz and Balbus, 2004; Wardle, 2007; Bai and Stone, 2013a). Accurately capturing the full nonideal physics is computationally difficult, and the numerical expense limits the feasible resolution
and / or run length, so it is common for studies to focus on a regime where only one or two effects
dominate, as we will do throughout this thesis. Figure 1.2 shows a “map of the disk”, which gives arough idea which effects are most relevant in each disk location, as well as pointing out the different
regions and effects studied in this thesis.

1.1.5.1

Ohmic Term

Physically, the Ohmic term arises due to collisions between charge-carrying species and neutrals, which causes the currents that generate magnetic fields to dissipate. Gammie (1996) considers
the damping of the MRI by Ohmic diffusion in a disk with scale height h, Alfvén speed vA and8

Figure 1.2 A cartoon “map of the disk”, emphasizing the most relevant non-ideal MHD effects ineach region. We also point out the various regions and phenomena studied by each chapter of thisthesis. Credit to Armitage (2007), and Simon et al. (2018) for the general idea for figures of thistype. A few purely notational choices that require clarification: 1. There are not actually breaks
in the particle layer as shown, this is just to allow room for text. 2. Planetesimal formation is notnecessarily specific to the radii shown here, it will generally happen throughout the radial extent
of the disk.

9resistivity η. By equating the minimum Ohmic damping rate, η/h2 , with the linear MRI growthrate on the same scale, vA /h, he argued that local turbulence would damp for resistivities exceedingηcrit ∼ hvA . As discussed above, this corresponds to an ionization fraction on the order of 10−13 ,which is not achieved via thermal nor plausible non-thermal ionization sources. Consequently,
η > ηcrit at the mid-plane on AU-scales. This damping effect at the mid-plane lead to a layered accretion model: MRI-active surface layers with efficient angular momentum transport, and a “dead”
mid-plane region with no MRI activity.
Modern studies usually quantify the Ohmic damping of the MRI in terms of a dimensionlessquantity called the Elsasser number:
Λ=

2
vaz
,
ηΩ

(1.3)

where the MRI is suppressed for Λ < 1. While the behaviour of the other two non-ideal terms
involves quite a bit of subtlety, the Ohmic term is more straight forward and binary: either the
MRI is permitted or it is suppressed.
The effect of Ohmic diffusion has since been well quantified by numerical studies. Simulations by Fleming and Stone (2003) use a model that included only Ohmic diffusion and verify
Gammie’s qualitative disk structure. However, they also showed that waves generated within the
MRI-active zone could leak into the dead zone, producing a non-zero level of turbulence and purelyhydrodynamic stresses. This leaking effect is the primary topic of study in chapter 2.

1.1.5.2

Hall Term

The Hall effect can create an electromotive force on a plasma due to a difference between
the electron velocity and the mean fluid velocity. Balbus and Terquem (2001) look at the effectsof the Hall term on protoplanetary disks and find that Ohmic diffusion dominates over the Hall
effect if the neutral density exceeds 1018 cm−3 or if the Alfvén speed is much less than the soundspeed. Dead zones dominated by the Hall Effect are thought to exist in the inner region of thedisc (∼ 1 AU), and have been shown by several studies to have an interesting behaviour withrespect to direction of the vertical magnetic field (Wardle, 1999; Simon et al., 2015b; Lesur et al.,102014; Bai, 2014, 2015). When the field is aligned with the spin axis of the disk, the Hall Effectis able to “revive” the dead zone: the MRI is replaced with the ”Hall-shear” instability (Lesuret al., 2014; Simon et al., 2015a; Kunz, 2008), resulting in stresses similar to the MRI (α ∼ 10−2 ).
However, when the vertical field is anti-aligned with the spin of the disk, the Hall Effect suppresses
the MRI a and the dead zone is largely true to it’s name. In this case, there is still a low levelof intermittent turbulence caused by the non-axisymmetric variant of the Hall-shear instability(Simon et al., 2015a), but the overall result is stress on the order of α ∼ 10−4 . This dichotomous
behaviour could have interesting consequences in disks. In chapter 3 we consider the ability of Hall
dominated dead zones to produce episodic accretion events.

1.1.5.3

Ambipolar Term

Ambipolar diffusion (AD) is caused by the breaking of the quasi-neutral assumption of ideal
MHD on small scales. The field drifts with the charged particles, causing it to drift relative to theneutral fluid. The ambipolar MRI can be suppressed in a similar fashion to the Ohmic MRI, butwith some added complications. Analytic work has shown that the relevant control parameter (atleast in the linear limit) for the growth of the MRI is the ratio of the rate of ion-neutral collisions
to the orbital frequency:
Am =

γρi
,
Ω

(1.4)

where γ is the drag coefficient (related to the ion-neutral collision rate) and ρi is the mass-density ofions. The numerator as a whole represents the ion-neutral collision rate. In contrast to the Elsasser
number for Ohmic diffusion, Am does not represent a strict criteria for MRI activity/inactivity.
Generally speaking, the MRI will be more damped at lower values of Am. However, it is important
to note that the strength of AD depends on the magnetic field strength itself, meaning that unstablemodes will always exists for a sufficiently weak field (Blaes and Balbus, 1994; Kunz and Balbus,
2004).
This damping has been the topic of several numerical studies as well, seeking to understandthe full, non-linear evolution of the MRI under these conditions. Hawley and Stone (1998) studied

11
this problem treating the ions and neutrals as separate fluids coupled with a collisional drag term.
They establish limits on the relevance of AD, showing that the fluid behaves essentially the same as
the ideal limit for collision frequencies on the order of 100Ω, while for collision frequencies . 0.01Ω
the neutral fluid is quiescent, unable to be affected by the motion of the ions. Bai and Stone (2011)
study the non-linear evolution of the ambiploar MRI in an unstratified shearing box using a single
fluid method. In contrast to the two-fluid approach, in the single fluid MHD limit, the ion-neutral
coupling is strong (which is shown to be appropriate for protoplanetary disks Bai 2011a), the ion
density is a fixed fraction of the neutral density instead of obeying it’s own continuity equation,and the ion-density makes up a negligible fraction of the total mass-density. Bai and Stone (2011)
showed that for any given Am value, one can always find a magnetic field strength at which the
MRI is present in its full nonlinear state. For smaller Am, weaker magnetic field is required, which
in turn leads to a lower level of turbulence. AD’s dependence on the field strength is defining
difference between Ohmic and Ambipolar diffusion (at least in terms of their effect on the MRI).
The first work to self-consistently model the AD-dominated outer disk with vertical gravity
included was carried out by Simon et al. (2013b) and Simon et al. (2013a). These authors carried out
local and vertically stratified calculations, in which AD was quite strong near the mid-plane. This
work included ionization of upper disk layers by FUV photons (Perez-Becker and Chiang, 2011a)
leading to a scenario similar to the Gammie (1996) model; active regions vertically sandwich a
region of weak or no turbulence. Furthermore, these works determined that a net vertical field is
required to predict high enough accretion rates to match the expected life-time of disks aroundyoung stars (assuming that disk accretion is at least approximately in steady-state). In chapter
4 we further explor the physics of disks with ambipolar damping zones, primarily focusing on therole of vertical transport of energy and magnetic fields, which has not been addressed previously.

1.1.6

Observations of PPD Turbulence

In recent years, it has become possible to put direct observational constraints on turbulencein the outer regions of disks with ALMA. Molecular line broadening measurements of HD163296

12
place an upper limit on the velocity dispersion (δv) of the turbulence throughout the disk at 0.05cs(Flaherty et al., 2015, 2017). Similar measurements of the TW Hydra system also show low levelsof turbulence: δv . 0.08cs (Flaherty et al., 2018). These measurements mainly measure the outerregion of the disk (tens - 100 AU), where ambipolar diffusion is the dominant non-ideal effect.
Simulations relevant to this region of the disk predict characteristic velocity perturbations in therange of 0.1 − 1cs in the surface layers of the disk, definitively inconsistent with the observations(Simon et al. (2015a), chapter 4). These observations have been somewhat baffling to theorists,
and this discrepancy has become known as the “weak turbulence problem”.
These observations motivate a more precise understanding of the physics of the outer disk
regions for two primary reasons. First, the presence of weak turbulence has only been established
for a few systems; if disks with strong turbulence do exist, it behooves us to explore this turbulencein more detail. Second, if indeed all disks show significantly weak turbulence, then we need tounderstand precisely why the MRI is prevented from operating. It is a very robust and vigorous
instability. From a theory standpoint, it is much easier to explain its presence than its absence.
Two basic scenarios are possible: 1. The surface layers of the disk are somehow prevented frombeing radiatively ionized and the magnetic field is weak (Simon et al., 2018). 2. Some unknown
process is actively arresting the MRI in the surface layers.

1.1.7

Winds from PPDs

Following the focus of this thesis as a whole, this introduction has mainly focused on MRI
turbulence as the driver of angular momentum transport and accretion in PPDs. However, disk
winds also merit mention. Whereas turbulence provides a mechanism for redistributing angularmomentum within the disk, winds are able to effectively remove angular momentum from thesystem entirely. It is likely that both winds and turbulence play an role in the angular momentum
evolution in PPDs.
The basis for most wind models in disks starts with Blandford and Payne (1982), who derive
the MHD solutions for a wind launched from a Keplerian disk with a poloidal magnetic field and13infinite conductivity. These solutions show that the magnetic field lines, which are anchored in the
conducting disk, act as rigid wires and will centrifugally accelerate gas away from the disk as longas the field lines are inclined at least 30 degrees away from the spin axis. For a very simplistic
illustration, see figure 1.2. This wind then exerts a torque on the disk, removing angular momentum
and allowing the disk to accrete. While this mechanism may apply fairly directly to high-energy
accretion disks, winds from PPDs require further consideration due to their poorly-ionized nature.
Bai et al. (2016) constructs a model for magneto-thermal winds from poorly-ionized PPDs,
integrating results from local simulations which include non-ideal MHD effects. While the basic
principles of the Blandford-Payne mechanism still apply, in this case the wind is launched solely
from the ionized surface layers of the disk, with the mid-plane playing essentially no role. Thismodel demonstrates that winds can drive reasonable accretion rates in the inner disk (∼ 10−8 M /yr
at 1 AU). Furthermore, mass loss due to this wind is quite high, on the order of 0.1 − 1 times theaccretion rate itself. However, the strength of this wind scales strongly with the net-vertical fieldin the disk, which is largely unconstrained.

1.1.8

Episodic Accretion in PPDs

While most studies of PPDs are designed to have relevance to either planet formation or
angular momentum transport, another active area of research is episodic accretion. Observationsof young stellar objects in the optical and infrared have revealed a wide range of out-bursting
behaviours in which the disk’s apparent brightness increases for a brief time. These events vary
greatly both in amplitude and time-scale (Herbig, 1977; Hartmann and Kenyon, 1996). Historically,
these events have been placed into two categories: FU Orionis objects (FUors) and EX Lupi
objects (EXors) (Herbig, 1989). FUor outbursts generally occur for several to tens of years withaccretion rates increasing by several orders of magnitude compared to the quiescent state, althoughproperties of individual objects vary greatly. In contrast, EXor outbursts are shorter and weaker,
typically bursting for months to years. They also occur more frequently, typically having a few years
between bursts. As the sample of observed out-bursting young stars has expanded, the wide range

14of parameters for both categories has caused the distinction between the two cases to become lessclear. Several observations of “in between” events have led to the idea that these events may be partof a continuous spectrum instead of a dichotomy (Audard et al., 2014). While this classification isuseful from an observational standpoint, it is not clear that FUors and EXors are caused by distinct
physical mechanisms, or that all events of each type are caused by the same physical mechanism.
Modeling of the spectral energy distribution of FUor events, in combination with basic timescale arguments, has led theorists to conclude that these events are most likely from the innerdisc (Zhu et al., 2007, 2008). The viscous accretion time of the disc is roughly Ω−1 α−1 (h/r)−2 .
Plugging in some typical numbers (α = 0.1, h/r = 0.1) yields a viscous time of about 160 years at
1 AU, which is roughly consistent with the time-scale of FUors. Any explanation of these events
invoking phenomena on larger size-scales would be difficult to justify. One proposed mechanism
for triggering FUors is called the gravo-magneto instability (Martin and Lubow, 2011). In this
scenario, the inner edge of the dead zone builds up enough mass to become unstable to self gravity,
which causes turbulent heating, which increases the ionization fraction of the gas enough to trigger
the MRI. This model relies on the presence of a true dead zone with little turbulence, and could besensitive to turbulence induced in the dead zone via hydrodynamic mixing from the active layers(the effect studied in chapter 2).
EXors, occurring on even shorter time-scales than FUors, would naively be expected to
be triggered by something in the very inner disc as well, or by magnetospheric accretion effects.
Armitage (2016) suggests that the dichotomous behaviour of Hall dead zones (described in §1.1.5.2)
could lead to episodic and/or cyclic accretion with properties similar to EXor events. While themagnetic field is aligned with the spin axis of the disk, the dead zone will be “active”, and theturbulent stress will drive a relatively high accretion rate. While the field is anti-aligned, the deadzone will be “dead”, and the lack of turbulence will cause the accretion rate to be much lower.
This change in net-field polarity could be caused by the stellar magnetic field diffusing though theinner disc, with the stellar field changing sign on some characteristic timescale. We explore thismechanism in detail in chapter 3.

15
Now that we have given an overview of the physics involved for the gaseous component of thedisk, we can move on to take a look at the behaviour of solids in the disk, and most importantly
the interaction between the gas and solid components.

1.2

Solids in Disks
Solid objects (or “particles”) interact with the fluid component of the disk via aerodynamic

forces, with the gas providing a drag force on particles that have a velocity relative to the local
velocity of the gas. Broadly speaking, the drag exerted on a particle with size s will be proportional
to it’s surface area (scaling as s2 ), while the particle’s inertia (ie. mass) will be proportional toit’s volume (scaling as s3 ). So as a particle’s size increases, the relative strength of aerodynamic
forces will decrease. Dust grains will be very easily blown around by the fluid’s movement, while
object’s on the order of kilometers will not be greatly affected by the fluid. The strength of the
aerodynamic force on particles is typically quantified by the friction time:
tfric =

mv
,
FD

(1.5)

where m is the particle mass, v is the relative velocity, and FD is the drag force. This is alsofrequently parameterized in terms of the “dimensionless stopping time”, given by τs = Ωtfric . Small
particles, which are very well coupled to the fluid, will have τs  1. In contrast, large objects
that do not particularly care about the fluid will have τs  1. With a basic framework for talking
about drag in place, we will now briefly discuss two important particle phenomena in disks: verticalsettling and radial drift.

1.2.1

Vertical Settling

Unlike the fluid, which has a vertical density distribution determined by balancing the verticalcomponent of the star’s gravity with the fluid’s pressure gradient, particles do not feel pressure
gradients. As a result, if we consider a non-turbulent disk, we would expect to find all solid particles
in a very thin layer at the mid-plane of the disk, feeling only the effect of gravity (Weidenschilling,161980). The presence of turbulence complicates the situation, providing an effective diffusion of
the particle layer in the vertical direction that will depend on the stopping time of the particles
considered and the magnitude of the turbulent velocities in the disk. If the turbulent diffusion
coefficient for the particles is assumed to be the same as the turbulent diffusion coefficient for
momentum (ie. viscosity), and the particles are relatively well coupled to the fluid (τs < 1), thena simple expression can be derived for the scale height of the particle layer hp :
hp
=
h

r

α
.
τs

(1.6)

In other words, the settling changes the scale height of the particle layer relative to the gas by apfactor of α/τs . The density of the particle layer is then enhanced by the inverse of this factor.
For example, a τs = 10−2 particle layer in a turbulent disk with α ∼ 10−2 will not have enhanced
density. However, a dead mid-plane can have orders of magnitude less turbulence, creating a muchthinner dust layer and enhancing the particle density by a factor ∼ 100 for a stress of α ∼ 10−6 .
This density enhancement is crucial for many planet formation models as it allows for much fasterparticle coagulation, and is another reason that understanding the level of turbulence at the midplane of PPDs is essential.

1.2.2

Radial Drift

While settling arose due to particles not feeling pressure gradients in the vertical direction,
another important effect called radial drift arises due to particles not feeling pressure gradients in
the radial direction. If we choose to work in the frame of a particle orbiting the central star, theparticle experiences only 2 forces which balance each other: gravity and the centrifugal force. This,pof course, leads to the particle orbiting the star at the familiar Keplerian velocity: vk = GM∗ /r.
In contrast, if we consider a fluid element in the same frame, it will experience 3 forces: gravity,
the centrifugal force, and an inward pressure gradient force due to the outward-decreasing pressure
profile of the disk (Weidenschilling, 1977). If the pressure profile of the disk is given by a power17
law: P ∝ r−n , then the orbital velocity of the fluid is given by 2 1/2
cs
vφ = vk 1 − n
.
vk

(1.7)

Notably, this orbital velocity is sub-Keplerian, with the difference between the two velocities usually
being on the order of 1%. However, even this small deviation has very interesting consequences:
particles in the disk are constantly experiencing a “headwind”. They are trying to orbit more
quickly than the fluid they are immersed in, and as a result they loose angular momentum and willdrift radially inward. This effect is felt by particles of all sizes, but is most prominent for objects
in an intermediate regime: τs ∼ 1 (Weidenschilling, 1977; Takeuchi and Lin, 2002). Physically, thiscorresponds to objects roughly 1 meter in size, which will drift inwards to the star on a remarkably
short time-scale of ∼ 1000 years. This is the basis of the “1-meter barrier” for planet formation.
Considering a zeroth-order approximation of planet formation wherein small particles coagulate
over time to form larger and larger objects, as soon as an object approaches roughly 1 meter insize, it will fall into the star and halt any further coagulation unless it is able to rapidly grow to asize-scale on the order of kilometers. This poses probably the hardest problem in planet formation:
theorists need a mechanism in this intermediate regime that is able to grow objects from cm to km
scales very quickly. We will revisit this problem in the following section.

1.3

An Overview of Planet Formation
To put it simply: explaining the entire planet formation process is an incredibly hard problem.

One needs to start with dust grains on the scale of microns and end up with a giant rock (and/orball of gas) on the scale of thousands of kilometers or larger. That’s a span of at least 12 orders ofmagnitude in size with different physical processes dominating the formation process at different
size-scales. In this section, we give a very brief overview of the planet formation process, with anemphasis on planetesimal formation (the phase directly studied in this work).

18
1.3.1

Coagulation at Small Size-Scales

At micron scales, the evolution of particles (often called “grains” at this scale) is dominatedby collision physics: do grains stick together or fragment when they collide?. The physics ofthe sticking efficiency of grains is complex and dependent on many factors including size, shape,
composition, and collision velocity. However, this is a rare area of astrophysics where laboratory
experiments have been incredibly helpful in determining this behaviour. Dominik et al. (2007) gives
a review of these experiments in combination with coagulation models (eg. Dullemond and Dominik(2005); Schmitt et al. (1997)) to give a useful picture of the early stages of planet formation. In
growing from micron to mm scale, surface electromagnetic forces are able to stick grains together
very efficiently, even for relatively large collision velocities (∼ 1m/s for spheres, or up to ∼ 100 m/s
for irregular grains). Once mm-cm scales are reached, the expected relative velocity of collisions
increases and composition and grain structure become more important, with fragmentation being
overall a more likely outcome of collisions (Zsom et al., 2010). If grains are not particularly compact,
then more collision energy can be dissipated. Wurm et al. (2005) shows that, on this size-scale,collision velocities of up to 10 m/s or higher can still lead to net-growth of large grains while also
producing many small fragmented grains.
After particles have coagulated to cm-size scales, we approach the 1-meter barrier posed by
the radial drift effect described in the previous section. The next growth mechanism must growobjects very quickly (thousands of years) to kilometer sizes. This phase of formation is known asplanetesimal formation, where a planetesimal is defined to be an object with relevant self-gravity
whose dynamics are no longer greatly affected by gas drag. One promising theory for planetesimalformation that will be discussed in this thesis involves a combination of a phenomenon called thestreaming instability and particle self-gravity.

19
1.3.2

Planetesimal Formation via the Streaming Instability

The streaming instability (SI) is a local, linear instability in Keplerian disks. The SI is a resultof the aerodynamic coupling of solids to the fluid and the momentum-feedback of the solids onto thefluid. The SI was first quantified by Youdin and Goodman (2005), who use a 2-fluid approach toconsider the interaction between the gas and solids in a unstratified disk with uniform particle sizes,
no turbulence, and no self-gravity. While there is no particularly obvious intuitive interpretationof the instability, it generally concentrates solids within the disk into streams and clumps. A very
non-physical analogy would be that of birds flying in flocks to reduce wind resistance. Youdin
and Goodman (2005) show that the instability is robust, existing across a wide range of particle
parameters. Growth rates are maximized when the stopping time of the particles is near the
dynamical time (τs ∼ 1) and the solid-to-gas ratio (ie. metallicity) is near 1. If we are considering
the disk as a whole, this is a drastically super-physical value for the metallicity, which would be
expected to be near the solar value of ∼ 0.01. However, we must keep in mind that vertical settling
can concentrate solids to the mid-plane of the disk, creating a local region with a greatly enhanced
metallicity. Growth rates under these optimal conditions can be up to ∼ 0.1Ω, which is slower thanthe dynamical time but very importantly faster than the time-scale for radial drift. The abilityof the SI to concentrate solids into clumps on a quick time-scale excited theorists, giving a new
potential mechanism for rapid planetesimal formation.
Since the initial formulation of the SI, considerable numerical work has been done to study
the non-linear regime. Johansen and Youdin (2007) demonstrate that the SI is capable of producing
particle over-densities of greater than 100 times the background particle density, and also drives
turbulence at a very low level (α ∼ 10−5 ). While these particle over-densities are promising for
initiating the collapse of gravitationally bound planetesimals, there are non-trivial constraints on
the parameters that lead to this clumping. In particular, the degree of clumping is a very strong
function of the solid-to-gas ratio: changing the metallicity from solar to twice solar can produce aten-fold increase in the typical particle over-density (Johansen et al., 2009).

20
To carry the planetesimal formation process one step further, several recent studies haveincluded particle-self-gravity (Simon et al., 2016; Abod et al., 2018; Schäfer et al., 2017). With
this included, particle clumps that are concentrated to sufficient densities by the SI become boundby their own self gravity and form a bonafide “planetesimal”. These studies then characterize
the mass-spectrum of the planetesimals, generally finding that the spectrum is well-described by apower law of the formdN
∝ Mp−p ,
dMp

(1.8)

where the power law index p is measured to be ∼ 1.6 ± 0.1. We save a detailed discussion of these
spectra and SI clumping for chapter 5.

1.3.3

Planetesimals to Terrestrial Planets

Once planetesimals have formed, interactions with the gas component of the disk do stillplay a role: gravitational interactions with turbulent density fluctuations will excite the planetesimal’s orbits (more on this in §1.3.4), while aerodynamic interactions will damp orbits. However,the primary driver of further formation in this phase is mutual self-gravity between planetesimals.
Gravitational focusing will enhance the effective cross-section for collisions of planetesimals by a fac2 /σ 2 ), where v
tor of (1+vesc
esc is the escape velocity from the surface of the planetesimal and σ is thevelocity dispersion of planetesimals. This enhanced cross section makes planetesimal-planetesimalcollisions quite frequent. Lissauer and Stewart (1993) derive growth rates for planetesimals in this
scenario. Most notably, the growth rate is proportional to the surface density of solids in the disk,
and the effectiveness of this type of growth is dependent on a fairly low value of σ (ie. planetesimals
are generally on circular orbits). If this velocity dispersion is sufficiently low, the largest planetesimals in the disk will exhibit runaway growth as their gravitational focusing factor continues toincrease as they accrete more planetesimals, making accretion more efficient as they grow. Aftertens of millions of years, the end result is a terrestrial-like planet (or gas-giant core).

21
1.3.4

Turbulence, Dead Zones, and Planet Formation

In this section we will briefly emphasize several areas in which a detailed understanding of
PPD turbulence plays an important role in understanding planet formation, further motivating thestudy of non-ideal MHD effects and dead zones.
In the early coagulation phase, the effective diffusivity generated by turbulent velocity fluctuations within the dead zone sets the thickness, and hence the density, of settled particle layers.
In addition, turbulence plays the dominant role in determining the collision velocity of small grains.
The presence of dead zones is friendly to the formation of large grains in this regime, as a lessturbulent mid-plane will lead to a denser particle layer, more frequent collisions, and collisions thatare more likely to result in coagulation rather than fragmentation.
At planetesimal size scales, surface density fluctuations in the disk due to turbulence generate stochastic gravitational torques on planetesimals (Ida et al., 2008; Gressel et al., 2011, 2012;
Okuzumi and Ormel, 2013). The net-result of these torques is a more excited orbit (higher eccentricity and inclination), which leads to larger relative velocities between planetesimals (the σparameter in the previous section). The excitation of the orbits scales roughly proportional to thedensity perturbation δρ. In a fully turbulent disc, the relative velocities of the planetesimals areabove the threshold for catastrophic disruption for bodies on the order of 1 km. For this reason,
planet formation models that rely on runaway growth of bodies of this size fail in a fully turbulent
disc. Gressel et al. (2011) find that the presence of a largely non-turbulent dead zone near the
mid-plane reduces the excitation of the planetesimals by a factor of 20, putting the relative velocity
of km sized bodies in between the catastrophic disruption threshold for weak and strong materials(Stewart and Leinhardt, 2009). This allows planetesimals composed of strong materials to survive
collisions and aggregate gradually.

22

1.4

Organization of This Thesis
This thesis consists of 6 chapters: this introduction, 4 chapters describing distinct projects

on topics related to turbulence and magnetic fields in protoplanetary disks, and finally a summaryand discussion of future work. Each of the 4 middle chapters will have a scope roughly equivalent
to that of a typical journal article. These chapters were largely written to each be understandable
as a stand-alone work, and as a result there may be some small overlap from chapter-to-chapter(in particular in the introduction and methods sections).
This first chapter is an introduction that provides a general background of protoplanetary
disk physics and phenomena, focusing on areas relevant to this thesis. In addition to this overall
introduction, there will be an introduction to each chapter that will re-emphasize some of therelevant physics already described here, as well as put the chapter in context relative to other work.
The 3 chapters following the introduction each explore one of the 3 non-ideal (MHD) effects as it
relates to disk turbulence and phenomena that stem from this turbulence: In chapter 2, we use 3D
shearing box simulations to study Ohmic dead zones, exploring the effect that the fluid motions inthe surface layers of the disk have on the dead zone. In chapter 3, we use a 1D disk model to study
the stellar magnetic field’s interaction with the disk, showing that the stellar magnetic cycle may
be able to produce episodic accretion events from Hall dead zones. In chapter 4, we use 3D shearing
box simulations to study ambipolar dominated dead zones, again emphasizing the importance of
the interaction between the active and dead zones and demonstrating the importance of verticaltransport processes in disks. In chapter 5, we again use 3D shearing boxes explore the effect that
turbulence has on planetesimal formation via the streaming instability. Finally, in chapter 6, wesummarize the most important conclusions of this thesis and brief discussion of future work.
As a note to the reader, I will be using the pronoun “we” to describe in the first-person what
I have done throughout this thesis, as is standard practice in writing multi-author journal articles.

Chapter 2

Turbulence, Transport and Circulation in Ohmic Dead Zones

2.1

Introduction
Turbulence and angular momentum transport in well-ionized accretion disks are generated by

the magnetorotational instability (MRI; Balbus and Hawley, 1991, 1998). In weakly ionized disks,
however, the efficiency of the MRI is modified by non-ideal magnetohydrodynamic effects — Ohmic
diffusion, ambipolar diffusion, and the Hall effect (Blaes and Balbus, 1994). These processes are
predicted to impact the dynamics of accretion in the cool outskirts of cataclysmic variable (Gammieand Menou, 1998) and Active Galactic Nuclei disks (Menou and Quataert, 2001), and are dominant
across most radii in protoplanetary disks (for reviews, see e.g. Armitage, 2011; Turner et al., 2014).
The simplest explicit model for the effect of the non-ideal MHD terms was introduced by
Gammie (1996), who considered the damping of the MRI by Ohmic diffusion in a disk with scaleheight h, Alfvén speed vA and resistivity η. By equating the minimum Ohmic damping rate, η/h2 ,with the linear MRI growth rate on the same scale, vA /h, he argued that local turbulence would be
damped for resistivities exceeding ηcrit ∼ hvA . For AU-scales in protoplanetary disks this resistivity
corresponds to an ionization fraction of the order of 10−13 , which is larger than that expected at
the disk mid-plane via thermal or plausible non-thermal ionization sources. These considerations
led to the prediction of an Ohmic dead zone around the mid-plane in which η > ηcrit . In thismodel, accretion preferentially occurs in non-thermally ionized surface layers in which the MRI
drives turbulence. Simulations by Fleming and Stone (2003), using a similarly simplified modelthat included only Ohmic diffusion, verified Gammie’s qualitative disk structure. However, they

24
also showed that waves generated within the MRI-active zones could leak into the dead zone,
producing turbulence and purely hydrodynamic stresses.
The presence of a dead zone can qualitatively change the nature of protostellar accretion andplanet formation. For accretion, the dependence of the vertically integrated stress α (Shakura andSunyaev, 1973) on the surface density Σ determines whether the disk, at some radius r, is able
to attain a steady-state for a specified accretion rate Ṁ . Simple models of Ohmic dead zones do
not in general admit steady-states for all Ṁ (Gammie, 1996), and are potentially vulnerable to a
gravo-magneto limit cycle instability (Martin and Lubow, 2011) that has been conjectured to be
the origin of episodic accretion in Young Stellar Objects (Armitage et al., 2001; Zhu et al., 2009).
For planet formation, the effective diffusivity D generated by turbulent velocity fluctuations within
the dead zone sets the thickness, and hence the density, of settled particle layers. Surface density
fluctuations, on the other hand, generate stochastic gravitational torques that may dominate theexcitation (and limit the growth) of planetesimals (Ida et al., 2008; Gressel et al., 2011, 2012;
Okuzumi and Ormel, 2013).
Although the relative importance of the non-ideal MHD terms varies with density, temperature, and degree of magnetization (Balbus and Terquem, 2001; Kunz and Balbus, 2004; Xu and Bai,
2016), in protoplanetary disks it is often the case that multiple non-ideal effects are simultaneously
important. At large disk radii (∼ 30 − 100 AU) ambipolar diffusion is dominant (e.g., Desch, 2004;
Kunz and Balbus, 2004; Simon et al., 2013b,a), but on AU-scales both Ohmic diffusion and theHall effect modify the mid-plane dynamics (Desch, 2004; Kunz and Balbus, 2004; Wardle, 2007;
Lesur et al., 2014; Bai, 2014; Simon et al., 2015a), while ambipolar diffusion is important in theatmosphere (Desch, 2004; Kunz and Balbus, 2004; Wardle, 2007; Bai and Stone, 2013a). Accurately capturing the full non-ideal physics is computationally difficult, and the numerical expense
limits the feasible resolution and / or run length. Here, following Fleming and Stone (2003) andOkuzumi and Hirose (2011), we instead study the model problem of a purely Ohmic dead zone,
which we further simplify by imposing a sharp vertical transition between the active and deadlayers. By examining the interaction of dead zone and active zone gas in this controlled manner,

25
we are able to better understand the physical effects that are observed in the simulations. We areparticularly interested in determining the functional form of α(Σ), which is the critical quantity
for disk stability, and in studying the nature of dead zone fluid motions and their implications forsettled particle layers.
The plan of this chapter is as follows. In §2.2 we describe the numerical methods and initialconditions used. §2.3 presents the results of the simulations, including the scaling of dead zone
properties with dead zone thickness and the nature of turbulence and large scale fluid flows in
the simulated disks. §2.4 discusses the implications for particle settling and growth if large-scale
circulatory flows, seen in our simulations, are present in real disks. We wrap up with our conclusionsin §2.5.

2.2

Methods

2.2.1

Numerical Algorithm

We compute the structure of Ohmic dead zones using Athena. Athena is a second-orderaccurate Godunov flux-conservative code that uses constrained transport (CT; Evans and Hawley,
1988a) to enforce the ∇ · B = 0 constraint, the third-order in space piecewise parabolic method(PPM) of Colella and Woodward (1984) for spatial reconstruction, and the HLLD Riemann solver(Miyoshi and Kusano, 2005; Mignone, 2007) to calculate numerical fluxes. A full description of theAthena algorithm along with results showing the code’s performance on various test problems canbe found in Gardiner and Stone (2005), Gardiner and Stone (2008), and Stone et al. (2008).
We employ the shearing box approximation to simulate a small, co-rotating patch of an
accretion disk. Taking the size of the shearing box to be small relative to the distance from thecentral object, we define Cartesian coordinates (x, y, z) in terms of the cylindrical coordinates(R, φ, z 0 ) such that x = R − R0 , y = R0 φ, and z = z 0 . The box co-rotates about the central object
with angular velocity Ω, corresponding to the Keplerian angular velocity at R0 .
The equations of resistive MHD in the shearing box approximation are,26∂ρ
+ ∇ · (ρv) = 0,∂t∂(ρv)
ĪP *) = 2qρΩ2 xî − ρΩ2 z k̂ − 2Ωk̂ × ρv,+ ∇ · (ρvv − BB + ¯∂t∂B− ∇ × (v × B) − η∇2 B = 0.
∂t

(2.1a)(2.1b)(2.1c)

Here ρ is the gas density, B is the magnetic field, P * is the total pressure (related to the magneticĪ is the identity tensor, v is the velocity, ηfield and the gas pressure, P , via P * = P + B · B/2), ¯
is the Ohmic resistivity, and q is the shear parameter defined by q = −d ln Ω/d ln R, which is 3/2
in this case of a Keplerian disk. An isothermal equation of state with a constant sound speed cs isused,
P = ρc2s .

(2.2)

We use a standard set of boundary conditions appropriate for the shearing box approximation.
At the azimuthal boundaries we exploit the symmetry of the disk and apply periodic boundary
conditions. At the radial boundaries we apply shearing periodic boundaries, such that quantities
that are remapped from one radial boundary to the other are shifted in the azimuthal direction by
the distance the boundaries have sheared relative to each other at the given time (Hawley et al.,1995). In the vertical direction, we employ the modified outflow boundary condition of Simon
et al. (2013b), which are the standard outflow boundaries but with the gas density extrapolated
exponentially into the ghost zones; this method reduces the artificial build-up of toroidal magnetic
flux near the vertical boundaries. Such outflow boundaries can lead to significant mass loss, particularly in the case of strong vertical magnetic flux (Simon et al., 2013a). However, such drastic
mass loss is not observed in our simulations, consistent with expectations from our magnetic field
geometry as described below. In general, these boundary conditions can break conservation, so
methods described in Stone and Gardiner (2010) are used in order to keep the correct variables
conserved. Furthermore, we employ Crank-Nicholson differencing to conserve epicyclic energy tomachine precision. A full discussion of these methods and Athena’s shearing box algorithm can

27
be found in Stone and Gardiner (2010).

2.2.2

Initial Conditions

We run all simulations using a domain size of Lx × Ly × Lz = 2H × 4H × 8H, where H isthe scale height, defined here as,
√
2cs
H=
.
Ω

(2.3)

We assume that the disk is initially in hydrostatic equilibrium. Balancing the vertical componentof gravity from the star with the gas pressure gives a Gaussian density profile:ρ(z) = ρ0 exp

−z 2
H2


.

(2.4)

We set ρ0 = 1 as the initial density at the mid-plane. A density floor of 10−5 is applied to keep
the Alfvén speed from getting too high (and thus the time step too small) in the region far fromthe mid-plane. The density floor also prevents β (defined below) from getting too low, which can√cause numerical problems. We set the constant isothermal sound speed cs = 1/ 2, and choose theangular velocity to be Ω = 1, which gives H = 1.
We initialize the magnetic field to be purely toroidal with constant β parameter of β = 400,
where β is the ratio of the gas pressure to the magnetic pressure,

β=

2P
.
B2

(2.5)

Note that a factor of 4π has been subsumed into the definition of the magnetic pressure, following
the standard definition of units in Athena (Stone et al., 2008).
The saturation level of MRI-driven turbulence is known to be a function of the net verticalmagnetic flux (Hawley et al., 1995; Salvesen et al., 2016), which can additionally give rise tooutflows (Bai and Stone, 2013b) and to qualitatively distinct Hall physics (Lesur et al., 2014; Bai,

28
2014; Simon et al., 2015a). Here, our primary goal is to study the relationship between active layer
turbulence and motions induced by this active layer in an idealized model of the dead zone, forwhich the actual saturation level of the MRI is of secondary importance. For simplicity, and to
reduce the parameter space that needs to be studied, we therefore choose a field geometry with nonet vertical magnetic flux. Consistent with this choice, we do not observe strong MHD outflows
from our simulations.
In order to seed the MRI, we apply random perturbations to the gas pressure and velocity.
Following the procedure used by Hawley et al. (1995), the perturbations are uniformly distributed
throughout the box and have zero mean. Pressure perturbations are a maximum of 2.5% of thelocal pressure, and velocity perturbations are a maximum of 5 × 10−3 cs . We employ a uniform
numerical resolution of 60 grid cells per H, giving a total resolution of 120 × 240 × 480 zones per
simulation.

2.2.2.1

Resistivity Profile

Five simulations were done, varying the ratio of the surface density of the active zone Σactive
to the total surface density Σtotal . The surface density of the active layer is physically set by abalance between recombination and ionization from sources such as cosmic rays (Gammie, 1996),
X-rays (Igea and Glassgold, 1999), and FUV photons (Perez-Becker and Chiang, 2011c). In mostcircumstances this balance is achieved on time scales that are shorter than the turbulent mixing time
between the surface and the interior of the disk (Bai, 2011a). It is then a reasonable approximation
to ignore time-dependent ionization effects, and take the column of ionized gas to be fixed in space
and time. Varying the active to total surface density ratio is thus equivalent to changing the sizeof the dead zone.
Physically, the resistivity varies with column density as a smooth function whose form is set
by the details of the ionization-recombination balance. To more accurately study the scaling ofdead zone properties with dead zone size, however, we ignore this complexity and instead createa simple dead zone by imposing a sharp vertical transition in the Ohmic resistivity using an error

29function. Specifically, we adopt,η = MIN(η1 , η2 )

(2.6)





ηmid
z + zcrit
η1 =
ERF
+1 ,
2
hη

(2.7)




ηmid
z − zcritη2 =1 − ERF
.
2
hη

(2.8)

where

and

Here ηmid = 0.005 is the resistivity of the mid-plane, hη = 0.2H controls how quickly the resistivity
transitions from the dead zone at the mid-plane to the active region, and zcrit controls the height at
which this transition happens. Note that zcrit is specified by the initial conditions and remains fixed
throughout the simulation, it is not adjusted on-the-fly to account for changes in the instantaneousvertical density distribution.
Depending upon the ionization source (cosmic rays or X-rays) the ratio of active to total
column at 1 AU in a purely Ohmic dead zone model is expected to be of the order of 0.1 to10−2 , if the underlying disk has a surface density comparable to the Minimum Mass Solar Nebula.
Substantially smaller values of the order of 10−3 are possible if mass accumulates in the deadzone, as occurs in some models for episodic accretion (Armitage et al., 2001; Zhu et al., 2009;
Martin and Lubow, 2011). As we show later, however, modeling extremely thin active layers isnumerically challenging, because the time scale for kinetic energy in the dead zone to saturate
becomes prohibitively long. With this constraint in mind, we are able to reliably model dead zones
whose active to total column varies between ∼ 0.5 and ∼ 0.03. The specific values of this ratio,
along with the corresponding values of zcrit , are shown in Table 2.1. The corresponding resistivity
profiles are plotted in Figure 2.1.

30

Run
S
MS
M
ML
L

zcrit /H
0.64
0.93
1.16
1.37
1.55

Active-to-Dead Mass Ratio
0.58
0.23
0.11
0.054
0.029

Table 2.1 Dead Zone Size

Figure 2.1 Resistivity η as a function of z for each of our simulations, as labeled in the legend. The
mass in the active zone as a fraction of the total mass varies between about 0.58 (for simulation“S”) and 0.029 (for simulation “L”).

31
The degree to which the resistivity suppresses the MRI can be quantified by the Elsasser
number, which compares the MRI growth rate to the resistive damping rate,

Λ=

va2
.
ηΩ

(2.9)

We show Λ as a function of z in Fig. 2.2. Where Λ > 1, the most unstable mode of the MRI is
unaffected by ohmic diffusion. Where Λ < 1, the growth rate of the MRI goes as ΛΩ (Sano andMiyama, 1999). The initial resistivity and magnetic field strength were chosen such that Λ < 1 at
the mid-plane to keep the MRI suppressed in the dead zone, while maintaining a sufficiently small
resistivity such that the diffusive time scale does not limit the time step.
The local wavelength of the most unstable MRI mode is the greater of the ideal MRI wavelength and the resistive MRI wavelength, given as (Sano and Miyama, 1999),

2πvaΩ
λidealλresistive =.
Λ
λideal =(2.10)(2.11)

With a resolution of 60 zones per H, these MRI modes are well resolved in the initial conditions, with
20 zones per MRI wavelength at minimum in the active region and about 40 zones per wavelength
in the dead zone.

2.3

Results

2.3.1

General, Large-scale Properties

Illustrative properties of the simulations are shown in Figure 2.3 and Figure 2.4. Figure 2.3
shows a rendering of the velocity fluctuations in the simulation domain, while Figure 2.4 shows
space-time diagrams of the Maxwell and Reynolds stresses for the small and medium dead zone
cases (note that the stress is not normalized to the mean local density, and hence is not equivalent
to a local α value). It is clear from either figure that the imposed resistivity profile delivers the

32

Figure 2.2 The dependence of the Elsasser number, Λ on z for each of our simulations, as labeled
via the legend. The solid black horizontal line shows where Λ = 1. The dead zone is defined here
as the region where Λ < 1.

33
expected morphology of an idealized Ohmic dead zone. The Maxwell stress saturates rapidly to adominant level in the active zone, and is damped to very low levels in the dead zone. The Reynolds
stress starts low in the dead zone and then builds up over time. Fluctuations in the Reynolds stressin the dead zone ultimately build up to be as high or higher than in the active region.

2.3.2

Saturation Time-scales

Figure 2.5 shows the time evolution of the kinetic energy, Reynolds stress and Maxwell stressin the active and dead zone regions, for the run with a medium-sized dead zone. The stresses
in the active zone saturate on a time scale of 10-20 orbits, consistent with previous studies of
MRI-driven turbulence (e.g., Simon et al., 2009). The properties of turbulent and non-turbulent
fluid motions in the dead zone, conversely, evolve in a substantially different way. The dead zone
Reynolds stress appears to initially saturate on a rapid time scale, but continuing the run reveals
long (∼ 102 orbit) time scale, large amplitude fluctuations. The kinetic energy, meanwhile, grows
essentially monotonically and saturates only after approximately 300 orbits of evolution. (Thedead zone Maxwell stress also continues to evolve for a long time, but at such a low level as to bephysically uninteresting.) As we will discuss later, the bulk of the kinetic energy is associated with
non-turbulent fluid motions, and hence the kinetic energy evolves in a distinctly different manner tothe Reynolds stress. To quantify how the saturation time scale depends on the dead zone thickness,
we take the most conservative (slowest saturating) measure and define the saturation time to bewhen the average mid-plane kinetic energy did not increase by more than 3 percent over 20 orbits.
The saturation time increases as the mass of the dead zone increases relative to the active zone, as
shown in Figure 2.6 and Table 2.2. This is qualitatively consistent with the expectation that energy
and turbulence are being generated in the active layer and injected into the dead zone, so that alarger dead zone takes longer to saturate. However, the detailed evolution of the kinetic energy as
a function of time (both for the dead zone shown in Figure 2.5 and the other runs) proved hard
to fit with a simple quantitative model of energy injection plus damping. Clearly, though, somephysically realistic dead zones — for which the ratio of the active to dead mass might be as low

34

Figure 2.3 The log of the velocity magnitude for the medium sized dead zone run, plotted withdarker colors indicating higher velocity. The outer regions are clearly turbulent, transitioning into
more laminar flow at the mid-plane.

35

Figure 2.4 Space-time plots of the Maxwell stress = −Bx By (first and third plots) and Reynolds
stress = ρvx δvy (2nd and 4th plots) for the small (top 2) and medium (bottom 2) dead zone sizes.
The quantities are averaged in the x and y directions. Note that the color scales are kept constantfor like-quantities across the runs.

36

Figure 2.5 Time evolution of the volume averaged Kinetic Energy, Reynolds Stress, and Maxwell
Stress (from top to bottom), for the medium sized dead zone. The average over the dead zone istaken from |z/H| = 0 to 0.5, and for the active zone from |z| = 2.0 to 3.3. In addition, the plots
were smoothed by time averaging over 10 orbits.

37
as 10−3 — would require very long simulations, with well-controlled numerical dissipation — tomodel faithfully.

2.3.3

Stresses and energies in the saturated state dead zone

Having defined the saturation time for the different dead zone sizes, we compute timeaveraged vertical profiles of the Reynolds and Maxwell stresses using data from the remaining
portion of the runs. These stress profiles are plotted in Figure 2.7. The mid-plane Maxwell andReynolds stresses both decline as the thickness of the dead zone increases, but the Maxwell stressis always negligible. Only the Reynolds stress is interesting for physical purposes.
If Reynolds stress in the dead zone derives exclusively as a consequence of energy transported
there from the active zone, then its magnitude ought to be set by an equilibrium between energy
injection from the active zone and damping (which might be turbulent dissipation in the deadzone, or leakage of dead zone energy back into the active zone). We therefore expect a scaling of
the mid-plane Reynolds stress with the ratio of the mass in the active zone to that in the deadzone. The derived scaling is shown in Figure 2.8. A linear fit to the hρvx δvy i-maz /mdz relation
provides a good description of the numerical results at the disk mid-plane. As can be seen in
Figure 2.7, however, the mid-plane stress becomes extremely small for the thickest dead zone, andthis is reflected in the linear fit predicting that the mid-plane stress goes to zero at finite dead zone
thickness. We do not expect this to be strictly correct, but the results suggest that Ohmic deadzones have physically negligible mid-plane stresses once the active zone column becomes less thanor of the order of 1% of the total disk column.

2.3.4

Transport Efficiency

The ratio of Reynolds stress to kinetic energy (KE) is indicative of the efficiency of angularmomentum transport. The ratio is defined as,hρvr vφ ihreynolds stressi=
hKEi
hρv 2 /2i

(2.12)

38

Figure 2.6 The saturation time as a function of the ratio of the dead zone mass to the active zone
mass.

39

Figure 2.7 Time averaged vertical profiles of the Reynolds stress (left) and Maxwell stress (right)
for the different sized dead zones. The sharp transition seen near |z| = 3 in the Reynolds stress for
the ML and L runs corresponds to where the r modes (discussed in detail in §2.3.8) shut off. The r
modes have a larger spatial extent in the ML and L runs due to the larger dead zone. Therefore thelocation of the transition feature happens in a region of lower Reynolds stress, making the features
most prominent in these runs.

40

Figure 2.8 The average Reynolds stress near the mid-plane as a function of the ratio of the active
zone mass to the dead zone mass. Note that the mass ratio on this x-axis is inverted relative tofigures 6 and 11 to show the linear relationship.

41
where the brackets indicate an average over a given spatial volume. Table 2.2 contains this efficiency
diagnostic for the different runs, time averaged in the saturated state. The dead zone spatial average
is taken from from |z/H| = 0 to 0.5, while the active zone is from |z/H| = 2 to 3.3. These limits
were maintained across the runs to compare like quantities.
The efficiency in the active zone is about 0.17, independent of the size of the dead zone, while
the efficiency in the dead zone scales strongly with the size of the dead zone. The (hydrodynamic)transport efficiency for the medium sized dead zone is about an order of magnitude less than in
the active region, where ideal MHD conditions apply. There is a substantial further decline for thelarger dead zones. The large kinetic energy content of the dead zone is thus largely irrelevant as
far as angular momentum transport goes, implying some combination of (a) non-turbulent motionsand (b) turbulence that at the mid-plane is inefficient at transport.

2.3.5

Density Perturbations

The magnitude of density fluctuations is of interest, especially in the mid-plane where stochastic gravitational forces from such fluctuations excite the velocity dispersion of planetesimals (Idaet al., 2008; Gressel et al., 2011; Okuzumi and Ormel, 2013). We define the amplitude of density
perturbations as,ρ(x, y, z, t) − hρixy (z, t)δρ(x, y, z, t) =.
ρhρixy (z, t)(2.13)The time and volume averaged density perturbation is,δρδρ(z) = h| (x, y, z, t)|ixy ,ρρ(2.14)

where the over-bar represents a time average and the brackets represent a volume average over the
given dimensions.
The vertical profile of the density perturbations in the saturated state can be seen in Figure 2.9
for each different run. The perturbation in the active layer is relatively consistent across the different

42runs, peaking at about ' 0.25 (i.e. moderately compressible). However, the perturbation at the
mid-plane decreases as the dead zone size increases, being as high as 0.02 for the smallest deadzone and 3 × 10−3 for the largest. This diagnostic again reveals that, while there is a large amountof kinetic energy in the dead zone, the fluid motions are essentially incompressible and are notturbulent — even hydrodynamically — in the same way as the active zone.

2.3.6

Accretion Rates

Disks with dead zones are potentially subject to outbursts as the dead zone gains mass thatcan rapidly accrete once mid-plane temperatures are high enough for the MRI to set in (Armitage
et al., 2001; Zhu et al., 2009; Martin and Lubow, 2011). These outbursts can be avoided if the deadzone is sufficiently active, such as in the steady state disk models of Terquem (2008). However,such models generally require that the mass accretion rate through the dead zone be larger than in
the active region. The simulations provide a means of estimating the accretion rates through theRactive and dead zones by performing integrals ∝ ρα(z)dz. The dimensionless accretion rate froma single, constant z plane is calculated as follows,hρvx δvy − Bx By ixy
ṁ(z)
3πc2s
=
hρixy
∆z,
mtot Ω
mtot Ω2
hP ixy(2.15)

where mtot is the total mass of the shearing box. In Figure 2.10, we plot this accretion rate profile
for each run. We then integrate Equation 2.15 over z to calculate the ratio of the mass accretion
rate through the dead zone to the mass accretion in the active zone — parameter f in equation
11 of Martin and Lubow (2014) — as a function of the ratio of the dead to active zone masses.
This ratio is shown in Figure 2.11. The result is a non-monotic relationship in which f is less thanunity for all plotted cases. For the three smallest dead zones, the dominant effect is that adding
dead zone mass increases f simply because there is more mass present in the dead zone to accrete.
However, for the largest two dead zones, the dominant effect is that the stress is very low, whichcauses f to be very small. Since f is not large, we expect then that outbursts will generally occur(see Fig. 5 of Martin and Lubow (2014)).

43

DZ Size
S
MS
M
ML
L

Sat. Time
120
230
290
430
490

Time Averages
150-230230-300
250-450
400-700
500-800

Reynolds active
dead
40.1
67.5
114
2250
4550

KE active
dead
11.6
10.2
13.8
51
106

T.E. Dead
0.048
0.023
0.02
0.00980.0041

T.E. Active
0.17
0.15
0.16
0.18
0.18

Table 2.2 Summary of time and volume averaged values in the steady state. Times are given inunits of orbits. In the rightmost two columns, ”T.E.” refers to the transport efficiency. We compare
the transport efficiency in the dead and active zones, as defined in equation 2.12

Figure 2.9 The density perturbation profiles for the 5 runs, averaged over the saturated state.

44

Figure 2.10 The non-dimensional accretion rate as a function of z for the different dead zone sizes,
averaged over the saturated state.

Figure 2.11 The ratio of the accretion rate in the dead zone to the accretion rate in the active zone
as a function of the mass ratio of the two zones. The accretion rate was averaged over the region
defined to be the dead zone in table 2.1, which changes from run to run. Time averaging was done
over the saturated state.

45
2.3.7

Structure of the Turbulence

To further quantify the turbulent structure as a function of height, we use autocorrelation
functions (ACF) to characterize the degree of non-axisymmetry present in the turbulence. The 2D
ACF in the x-y plane for a quantity f is defined as,

R
ACF(f (∆~x )) =f (t, ~x )f (t, ~x + ∆~x )d2 ~xR
f (t, ~x )2 d2 ~x
,(2.16)

where ~x is the position in the x-y plane, ∆~x is the shift in this position, and the over bar indicatesan average over time. With this definition, the ACF is normalized by its maximum value at ∆~x = 0.
The time averages are performed during the saturated state for each run (see table 2.2). An ACF
was calculated in the x-y plane for each z and time, and then averaged over the desired range of z
and time-interval. Examples of time and volume averaged ACFs for the density perturbation andvz can be seen in Figures 2.12 and 2.13 respectively.
From the figures it is clear that the structure of the turbulence (or lack thereof) is a strong
function of z, as expected for a disk with a dead zone. In the active zone the density perturbation decorrelates over a short distance scale, and the ACF is a small, significantly tilted ellipse, consistent
with previous studies (Guan et al., 2009; Simon et al., 2012). Near the mid-plane the density
fluctuations are almost axisymmetric, consistent with the dynamics in this region being dominatedby spiral density waves. In an intermediate region near the lower parts of the active zone, thepresence of the dead zone causes the turbulence not to look like ideal MRI even where Λ  1.
The shape of the ACFs was quantified by fitting an ellipse to the density perturbation ACF.
The ellipse was fit such that it had the minimum possible area while containing at least 10 percentof the total value of the entire ACF. In the case of the ideal MRI, prior shearing box simulations
show that the expected angle of the density perturbation ACF is about −7◦ offset from the y-axis(Guan et al., 2009). This angle is plotted as a function of height for each run in Figure 2.14. In theactive zone (higher than z/H = 2.0 for all runs), the angle is near the expected value for ideal-MHD
MRI turbulence, while in the dead zone the turbulent fluctuations are more closely aligned along46

Figure 2.12 Autocorrelation functions for the density perturbation for the medium dead zone (top)and large dead zone (bottom). Moving from left to right, the ACF is calculated using planes |z| =0 to 0.5, 1.0 to 1.5, 2.0 to 2.5 and 3.0 to 3.5. Each plot is 2H on a side. The plots on the left showthat perturbations near the mid-plane get sheared out by the background flow, while perturbationsin the upper layer are smaller and less aligned with the y-axis.

Figure 2.13 Autocorrelation functions for the vertical velocity for the medium dead zone (top) andlarge dead zone (bottom). Moving from left to right, the ACF is calculated using planes |z| = 0
to 0.5, 1.0 to 1.5, 2.0 to 2.5 and 3.0 to 3.5. Each plot is 2H on a side. The large scale circulation
modes can clearly be seen in the leftmost plots, while the plots on the right corresponding to theupper layers are consistent with the structure of MRI-driven turbulence in the ideal MHD limit, asseen in other quantities and previous studies (e.g., Guan et al., 2009).

47
the y-axis. The degree of alignment generally increases as the dead zone thickness increases. We
also note that the angle measured at the very outer part of the disc (|z/H| > 3.0) is discrepant ascompared to ideal-MHD expectations. This could be due to the influence of boundary conditions
near the edge of the box, or because the magnetic pressure dominates the gas pressure in thisregion, quenching the MRI and thus changing the expected structure of turbulence in this region.
The ACFs of the magnetic field components are consistent with no magnetic activity inthe dead zone and the expectations of the ideal MRI from Guan et al. (2009) in the active zone.
However, the ACF of the vertical velocity displays distinctly different behavior. In the active zone,
the vertical velocity ACF has the modestly non-axisymmetric structure that is characteristic of
turbulence in the active zone, as seen in the ACF of other quantities. However, near the deadzone mid-plane, there are coherent, large-scale structures in the vertical velocity. We analyze theproperties of these structures further in the next section.

2.3.8

Large scale velocity structures in the dead zone

The structure of the large scale velocity field in the dead zone is shown in Figure 2.15 andFigure 2.16. Figure 2.15 shows a snapshot of the vertical velocity in the x-z plane (averaged over y)for the medium sized dead zone simulation. Figure 2.16 shows the vertical velocity in the x-y planeat the mid-plane at a number of different time slices. In the example shown, the vertical velocity isdominated by a mode with a horizontal wavelength of H, which is coherent vertically between aboutz = ±2H. Two of these modes fit within the radial extent of the shearing box. An inspection of
the time sequence of snapshots at z = 0 shows that the large scale velocity structures take the formof a low frequency oscillatory meridional circulation. A Fourier analysis shows that the frequency
is approximately 0.036Ω, implying a period of about 4.5 orbits. This dominant (primary) mode is
not the only such structure present. Figure 2.16 shows half a cycle of the oscillation. When thedominant mode is at a minimum, a weaker (secondary) mode with wavelength 2H/7 is apparent.
We have examined how the kinetic energy is partitioned between small-scale and large-scale fluid
motions as a function of height. As is clear visually, small scale kinetic energy dominates by a

48

Figure 2.14 The angle of the density perturbation ACF away from the y-axis for each of the deadzone runs as a function of z.

49
factor of a few in the active zone, while conversely large scale kinetic energy dominates by a factor
of a few in the dead zone.
The spontaneous development of large scale velocity structure in the simulated dead zone issomewhat surprising, and we have made various checks to assess whether similar structures would
be present in real disks. We first note that similar features were already seen in the shearing box
simulations of Oishi and Mac Low (2009), which were run with a different code (Pencil). Thoseauthors also found evidence for similar modes (though with reduced amplitude) in simulations
where the isothermal equation of state was changed to ideal. We are therefore confident that thelarge scale circulation is not a numerical artifact of the Athena code, though as we will discuss
below, the manifestation of disk modes in the shearing box approximation can differ from those
expected in a physical disk. We have also verified that the mid-plane of the simulated disk is stable
to within a grid cell over the course of a vertical oscillation cycle. This implies that our use ofa resistivity profile that is fixed in z, rather than being specified by column density directly, isunlikely to drive the modes unphysically.
We interpret the standing wave pattern in Figure 2.16 in terms of disk modes. We ignore theeffects of magnetic fields that play a key role in the active zone. We apply the approximation thatthe dead zone occupies the full vertical extent of the disk in computing the mode structure. Thisapproximation should approximately hold, since the active zone involves only about 11% of the disk
mass and lies well outside the disk mid-plane region where the dead zone resides. In the notation of
Lubow and Pringle (1993) the implied wave angular frequency from Figure 2.16 is ω ' 0.22Ω. Sinceplots in Figure 2.16 are nearly independent of the azimuthal coordinate y, we take these modes tobe axisymmetric. The only axisymmetric modes that have ω < Ω are r modes (Lubow and Pringle,
1993; Ogilvie, 1998). These modes derive their support from inertial and rotation forces. For disksthat behave adiabatically in the presence of vertical entropy gradients, buoyancy forces also playa role in confining the mode near the disk mid-plane (see Figure 11 of Lubow and Pringle (1993),
who referred to these modes as g modes). In the simulations of this chapter, buoyancy plays norole, since the disk is strictly isothermal.

50

Figure 2.15 A snapshot from the medium-sized dead zone run of the vertical velocity in the x-z
plane, averaged over y. The figure clearly depicts the presence of two vertically coherent updrafts/ downdrafts, stretching throughout the dead zone region.

51

Figure 2.16 The vertical velocity plotted for the x-y plane at the mid-plane of the medium-sizeddead zone run. Time advances from top to bottom by 3.5 Ω−1 each panel, showing half of a cycleof this standing wave pattern.

52
We suggest that the wave pattern in Figure 2.16 is the lowest order (in terms of verticalstructure) r mode. This mode has no vertical nodes in the vertical velocity. Its linear dispersion
relation for an isothermal Keplerian disk follows from equation (54) of LP93 with n = 0 (no verticalnodes in the vertical velocity), κ = Ω (which is equivalent to κLP93 = 1), and γ = 1 (isothermal
perturbations in an isothermal disk) and is given by
√
ω=

2Ω
2

q
p
2 + K 2 − |K| 4 + K 2 ,(2.17)√
where K = kx cs /Ω = kx H/ 2 is the scaled radial wavenumber.
The perturbations caused by the linear wave depend on a scale factor that we take to be thevertical velocity amplitude W . This amplitude depends on the strength of forcing from the active
layer. The standing wave properties are given by√ ωz 
2W
sin (kx x) cos (ωt),ΩH√2zvy (x, z, t) = −W sin (kx x) sin (ωt),
2Hvx (x, z, t) =(2.18)(2.19)

vz (x, z, t) = W cos (kx x) cos (ωt),
W ωzp(x, z, t) = 2
p0 exp (−z 2 /H 2 )H 2 Ω2(2.20)

× cos (kx x) sin (ωt),(2.21)

where p0 is the unperturbed disk pressure at mid-plane. Note that we have selected a phasing inwhich vz (0, z, 0) = W . In the low frequency limit for fixed W , the pressure perturbation vanishes,
while the wavenumber diverges in lowest order as K = Ω/ω. The velocity divergence does notvanish in this limit and is finite.
The radial boundary conditions used in the simulations imply that there is a discrete set ofradial wave numbers in the disk with
√
2K
2π`
kx =
=
H
Lx

(2.22)

for integer ` and Lx = 2H in the simulations presented here. The discreteness is due to the limited
box size. The wave structures given by Equations (2.18) - (2.21) satisfy this boundary condition

53
for the radial wavenumbers found in the simulations. The primary mode with radial wavelength H
and the secondary mode with radial wavelength 2H/7 identified in Figure 2.16 correspond to ` = 2
and 7 respectively. The level of excitation of these particular modes is determined by the overlapof their space and time dependences with these dependences for the turbulence.
For the primary mode that has ` = 2, wavelength H, and the wavenumber K =√2 π, ther mode dispersion relation given by Equation (2.17) evaluates to ω = 0.2147Ω, in good agreementwith the angular frequency ω ' 0.22Ω that is determined by the time changes in Figure 2.16 forthe same wavelength H disturbance. The agreement gives us confidence that we have properly
identified this dead zone disturbance as being due to the lowest order r mode.
We expect that the other (secondary) wave mode with radial wavelength 2H/7 to be roughly
similar to the primary r mode that has radial wavelength H and no vertical nodes in the verticalvelocity, n = 0, since both modes are excited by the active layer turbulence. The n = 1 r mode has
zero vertical velocity at the disk midplane at all times and therefore cannot explain the presence of avertical velocity in Figure 2.16. (The n = 1 mode may exist, but is is not detectable in Figure 2.16.)
A possible candidate for this secondary mode seen in Figure 2.16 is the n = 2 r mode which has
a nonzero vertical velocity at the mid-plane. The dispersion relation of Equation (2.17) applied tothis mode predicts a lower angular frequency that is about half the frequency of the n = 0 mode.
Due to the substantial frequency difference between the primary and secondary modes, we expect a
substantial phase difference between the modes when the primary mode has zero vertical velocity.
We then expect the secondary mode to have a nonzero vertical velocity at this time, as seen in the
middle panel of Figure 2.16.
We speculate on how the shearing box simulations might carry over to a full disk. For a fulldisk, the radial wave numbers would be continuous rather than discrete. Furthermore, the waves
would not generally be standing waves. A mix of inwardly and outwardly propagating r mode
waves would be excited by turbulence in the active layer. An outwardly propagating r mode withn = 0 becomes an outwardly propagating p (pressure) mode when it propagates past its corotation
radius that occurs where Ω(r) = ω (Lubow and Pringle, 1993). An inwardly propagating wave

54
remains as an r mode. An outwardly propagating r mode with n > 0 reflects at corotation andbecomes an inwardly propagating r mode. These waves would be subject to damping by shocks
and interactions with turbulence.
The flow pattern of an r mode is quite different for a disk with buoyancy. Buoyancy acts
to vertically confine r modes near the disk mid-plane. The confinement condition on z is that
N (z) < ω, where N is the vertical buoyancy frequency. For γ = 1.4, the r mode extends overheight |z| < 1.87ωcs /Ω2 , which for ω = 0.2Ω, implies |z| < 0.3H. The extent of driving the midplane motions could also be affected, if the base of the active region lies above this height. For any
reasonably sized dead zone, we would then expect a major change in the vertical flow properties ofa disk that undergoes adiabatic perturbations.

2.4

Implications for Particle Dynamics
Proceeding under the assumption that the wave modes seen in the simulations are present

in disks, we now ask what impact they would have on the distribution of solid particles. As Oishi
and Mac Low (2009) noted, the low frequency and non-turbulent nature of the vertical velocity
field means that there is no direct effect on particle settling for the small (dimensionless stopping
time τ  1) particles generally assumed as progenitors for planetesimal formation. However, a
low frequency vertical flow could have the effect of “buckling” a vertically settled particle layer,such that the mid-plane of the particle layer coincides with the mid-plane of the gas disk only fora fraction of an oscillation cycle. For realistic amplitudes, it is not clear that buckling of a midplane particle layer would have directly observable effects. A buckled particle layer would, however,suppress the rate of pebble accretion (Lambrechts and Johansen, 2012) on to larger bodies that are
aerodynamically or gravitational damped to the mid-plane of the gas disk.
To estimate the potential importance of this effect, we first determine the expected amplitudeof the vertical oscillations that would be experienced by a settled particle layer. To first order, thecirculation modes discussed in §2.3.8 have a sinusoidal time dependence at a given point in the
mid-plane,55vz = vz,0 cos(Ωc t)(2.23)

where Ωc is the oscillation frequency of the circulation mode and vz,0 is the peak vertical velocityof the mode at the mid-plane. Under this assumption, the maximum displacement of a particle
that is perfectly coupled to the gas is given by

∆zc =

vz,0
.
Ωc

(2.24)

The circulation mode parameters and this estimate of the displacement are given in Table 2.3 foreach different dead zone size. Because of the low frequencies involved, the modest velocities (of theorder of 0.1cs ) are able to induce vertical displacements that are of the order of H for all of thedead zone sizes considered.
The buckling of the particle layer will affect the efficiency of pebble accretion if the amplitude
exceeds the settled thickness of the layer. For particles with dimensionless stopping time τs the
thickness of this particle layer can be written as,

Hd
=
H

s

Dd,z
τs

(2.25)

where Hd is the scale height of the layer and Dd,z is the dimensionless vertical dust diffusion
coefficient. For particles that are well coupled to the gas (ts  Ω−1 ), this diffusion coefficient is
the same as the vertical gas diffusion coefficient Dg,z , given by,Ω
Dd,z = Dg,z = 2csZ ∞hvz (t)vz (0)idt.

(2.26)

0

The integrand is the auto-correlation function of the vertical gas velocity (Zhu et al., 2015). In orderto calculate this parameter, a high frequency (10 snapshots per Ω−1 ) 3D sample of data was output
from the medium-sized dead zone run over 5 orbits. To calculate the diffusion coefficient without
contamination from the circulation modes themselves, the circulation mode must be subtracted

56
from the velocity. Given that the mode is roughly constant across y and z near the mid-plane, thecirculation velocity vz,c was defined to be,vz,c (x, t) =hρ(x, y, z, t)vz (x, y, z, t)iy,z,ρavg(2.27)

where ρavg is the average density over the mid-plane region being considered. The random velocity
vz,random entering into equation (2.26) is then taken to be,vz,random = vz − vz,c .

(2.28)

Using this random vertical velocity, the average auto-correlation function over all points near
the mid-plane (−0.5 < z/H < 0.5) was calculated. This function, shown in Figure 2.17, is wellbehaved, and when integrated yields a dimensionless vertical diffusion parameter of 7.2 × 10−3 forthe middle sized dead zone. We roughly estimate this parameter for the other runs by assumingthat Dd,z scales with δv 2 at the mid-plane. Using equation 2.25, the scale height of the particle
layer can now be determined for a given stopping time τs .
If the displacement due to the circulation mode is comparable to the scale height of theparticle layer, the mode will have a significant effect on the structure of the layer. This criterion
can be used to define a critical stopping time τs,crit at which Hd and ∆zc are equal,τs,crit =Dd,z(∆zc /H)2(2.29)

These critical stopping time values are shown in Table 2.3. The results indicate that, in general,
particle layers with τs ≥ 10−3 will be significantly perturbed by these circulation modes.

2.5

Conclusions
We have used local numerical simulations to study the structure of Ohmic dead zones, that

may be relevant to accretion in protoplanetary disks (Gammie, 1996) and in the outer regions ofdwarf novae (Gammie and Menou, 1998). By setting up a simplified, yet well-defined model —57

Figure 2.17 The average auto-correlation function of the random vertical velocity for all points−0.5 < z/H < 0.5 as a function of the lag time. The vertical diffusion parameter is the area under
this curve.

DZ Size
S
MS
M
ML
L

vz,0
0.10
0.08
0.06
0.03
0.014

Ωc
0.039
0.038
0.036
0.024
0.024

∆zc /H
2.56
2.11
1.67
1.25
0.58

Dd,z
1.4 × 10−2
1.1 × 10−2
7.2 × 10−3
2.8 × 10−3
4.6 × 10−4τs,crit
2.2 × 10−3
2.6 × 10−3
2.6 × 10−3
1.8 × 10−3
1.4 × 10−3

Table 2.3 Parameters of the circulation modes, vertical particle diffusion parameters, and criticalτs values for each run.

58
that ignores the role of other non-ideal processes and has a sharp transition between magnetically
active and inactive zones — we were able to run long duration simulations at high resolution inorder to precisely quantify how dead zone properties scale with the thickness of the dead zone. We
found that:(i) The Reynolds stress, which dominates transport near the mid-plane, scales strongly withthe thickness of the dead zone. For modest dead zones, in which the ratio of magnetically
active to inactive column is of the order of 0.1 or higher, we find α ∼ 10−4 . Thicker
dead zones, in which the active column is only of the order of 1% of the total mass, havenegligible mid-plane stresses. As a result, small dead zones can have accretion rates thatare a significant fraction of the active zone accretion rate (10 to 30%), while larger deadzones have accretion rates around 1% of the active zone rate. This scaling supports models
in which mass can accumulate in the dead zone region, building up a reservoir of material
that can potentially accrete rapidly (Armitage et al., 2001; Zhu et al., 2009; Martin and
Lubow, 2011).
(ii) The structure of turbulence in the MRI-active region, at high z, resembles that seen in
ideal MHD simulations of fully active disks (Guan et al., 2009; Simon et al., 2012). The
hydrodynamic component of the turbulence becomes more axisymmetric and less efficientat angular momentum transport towards the mid-plane, especially for thick dead zones.
(iii) The time scale for saturation of the kinetic energy in the dead zone is long (reaching∼ 500 orbits for the thickest dead zone studied). Most of the kinetic energy is contained inlarge-scale, non-turbulent fluid motions, which take the form of an oscillatory meridional
circulation that is extended in the vertical direction. Similar fluid motions were observed in
independent simulations by Oishi and Mac Low (2009). We identify this motion with thelowest order r mode expected in a purely hydrodynamic disk model (Lubow and Pringle,
1993).
(iv) We expect the properties of the large-scale fluid motions to depend upon the nature of

59
the perturbations (isothermal or adiabatic), and upon where vertical driving of the mode
occurs. Moreover, the standing waves observed in our shearing box simulations would likely
be replaced by inwardly and outwardly propagating waves in a global disk model. However,
if models similar to those we observe are excited to significant amplitudes in more physicalmodels of protoplanetary disks, the resulting buckling of settled particle layers could reduce
the efficiency of pebble accretion for small (τ  1) particles.
Although purely Ohmic dead zones are an idealization, the inclusion of non-ideal physics in
disk simulations frequently leads to vertical stratification of turbulent properties, and regions of
the disk where mid-plane MHD stresses are negligible, as we have seen here. Further simulations
will be needed to determine if our results for the scaling of Reynolds stresses, and for the excitationof large scale fluid flows, carry over under more realistic conditions.

Chapter 3

The Interaction of the Stellar Magnetic Cycle with Hall Dead Zones inProtoplanetary Disks

3.1

Introduction
Observations of young stellar objects in the optical and infrared have revealed a wide rangeof cyclic and out-bursting behaviours that vary greatly both in amplitude and time-scale (Herbig,
1977; Hartmann and Kenyon, 1996). Historically, these events have been placed into two categories:
FU Orionis objects (FUors) and EX Lupi objects (EXors) (Herbig, 1989). FUor outbursts generally
occur for several to tens of years with accretion rates increasing by several orders of magnitude
compared to the quiescent state, although properties of individual objects vary greatly. In contrast,
EXor outbursts are shorter and weaker, typically bursting for months to years. They also occur
more frequently, typically having a few years between bursts.
As the sample of observed out-bursting young stars has expanded, the wide range of parameters for both categories has caused the distinction between the two cases to become less clear.
Several observations of “in between” events have led to the idea that these events may be part ofa continuous spectrum instead of a dichotomy (Audard et al., 2014). While this classification isuseful from an observational standpoint, it is not clear that FUors and EXors are caused by distinct
physical mechanisms, or that all events of each type are caused by the same physical mechanism.
SED modeling as well as basic time-scale arguments has led theorists to conclude that FUors aremost likely from the inner disk (Zhu et al., 2007, 2008). EXors, on even shorter time-scales, arealso expected to be triggered by phenomena at the inner edge of the disk, or by by magnetospheric61
accretion effects.
In this chapter, we will consider the interaction between the stellar magnetic field and thesurrounding accretion disk as a potential mechanism for episodic accretion. This interaction has
been explored in several contexts including the spin-down of the star due to magnetic torques fromthe disk (Ghosh and Lamb, 1979; Illarionov and Kompaneets, 1990), and the geometry and observational appearance of magnetospheric flows (van Ballegooijen, 1994; Aly, 1980), but there remains
a gap in the literature on the influence of the star’s magnetic field on the disk’s magnetic field anddisk physics. This is particularly important to understand given how fundamental magnetic fields
have been shown to be in understanding disk physics.
Any discussion of magnetic fields in PPDs will need to include non-ideal MHD effects, whichhave been shown to play an important role in understanding protoplanetary disk phenomena. In
general, these effects can allow the magnetic field to decouple from the gas such that the magnetorotational instability (MRI) (Balbus and Hawley, 1991) is no longer able to drive turbulence.
These regions of little to no turbulence are generally referred to as “dead zones” or “damping zones”(Gammie, 1996). Dead zones dominated by the Hall Effect are thought to exist in the inner regionof the disk (≤ 1 AU), and have been shown by several studies to have an interesting behaviour
with respect to direction of the vertical magnetic field (Wardle, 1999; Simon et al., 2015b; Lesuret al., 2014; Bai, 2014, 2015). When the field is aligned with the spin axis, the Hall Effect damps
the MRI and replaces it with the “Hall-shear” instability, which can still generate stress at roughly
the same level as the MRI itself (α ∼ 10−2 ) (Kunz, 2008; Lesur et al., 2014; Simon et al., 2015a).
However, when the vertical field is anti-aligned with the spin of the disk the Hall Effect is still
able to suppresses the MRI, but the Hall-shear instability is not able to produce stresses at thesame level. A non-axisymmetric version of the Hall-shear instability acts to produce weak andintermittent stresses, but the dead zone is largely true to it’s name.
Armitage (2016) suggests that this dichotomous behaviour of Hall dead zones could lead toepisodic and/or cyclic accretion. While the field is aligned with the spin, the dead zone will be“active”, and the turbulent stress will drive a relatively high accretion rate. While the field is anti-

62aligned, the dead zone will be “dead”, and the lack of turbulence will cause the accretion rate to
be much lower. This change in net-field polarity could be caused by flux from the stellar magnetic
cycle diffusing outward though the inner disk (Lubow et al., 1994) and reaching radii where theHall effect is relevant.
In this chapter, we develop a 1D, time-dependent model that co-evolves the disk’s surface
density and net-magnetic field in a self-consistent way such that the accretion physics can dependon the magnetic field. We then establish a prescription for the stellar magnetic field’s influence on
the disk. We use our model to develop a criterion for the parameters of star/disk systems that
could generate outbursts or cyclic behaviour. Finally, we apply our time-evolution model to show
the characteristics of these outbursts. This chapter is organized as follows: In section 2 we describe
the disk model that we will use. In section 3 we describe the methods we use to discretize themodel and evolve it in time. In section 4 we describe the magnetic field geometry of the star/disk
system and a model for the star’s contribution to the disk field. In section 5 we establish a criterion
for outbursts. In sections 6 we use our time-dependent model to study these events, and in section
7 we summarize our results.

3.2

The Model
We develop a model using a mean-field approach that self-consistently evolves the net mag-netic field and surface density of a turbulent accretion disk in one dimension.

3.2.1

Static Equations

Here we will give and briefly justify the equations that we will use to solve for the state ofthe disk at any given point in time (Shakura and Sunyaev, 1973). The sound speed is given by
c2s =

kb
Tc ,
µmp

(3.1)

where Tc is the temperature at the mid-plane of the disk, kb the Boltzmann constant, µ the mean
molecular weight, mp the proton mass, and cs the sound speed. Assuming hydrostatic equilibrium,

63
the vertical pressure scale-height of the disk is determined by balancing the vertical component of
the central star’s gravity and the vertical pressure gradient. For a non-self-gravitating disk, thescale height is given by
h=

cs
,
Ω

(3.2)

where Ω is the angular orbital velocity.
The turbulence in the disk gives rise to an effective viscosity ν, which we will parameterizewith the usual α-modelν = αcs h,(3.3)

where α is a dimensionless parameter that describes the level of turbulence. This viscosity allows
angular momentum to be transported within the disk such that some gas loses angular momentum
and accretes inward. The accretion rate is given by ṁ and in the steady state is related to theviscosity via
νΣ =

ṁ
,
3π

(3.4)

where Σ is the surface density.
In order to model the disk’s temperature there are two heating mechanisms to consider:heating due to accretion and heating due to incident radiation. In general, the inner disk will bedominated by accretion heating and the outer disk, where the surface area is much larger, will bedominated by heating from incident radiation. For this application of the model to the inner disk,
we will assume a simple active disk model with no incident radiation. We will consider the disk to
be heated at the mid-plane due to the loss of gravitational potential energy as mass falls inward.
This energy will travel to the surface of the disk through an optical depth τ , and will be radiatedfrom the surface. Under these assumptions, the surface temperature Td and central temperature
Tc are then given by
9
σTd4 = νΣΩ2 ,83
Tc4 = τ Td4 .
4

(3.5)
(3.6)

64
The optical depth is given by1
τ = ΣκR2(3.7)

where κR is the mean opacity. κR will (very generally) be some function of the gas density ρ andTc :κR = power law of ρ and Tc .

(3.8)

We will be applying this model to the inner several AU of a protoplanetary disk, where temperatures
would be expected to be in the range of 100 to 1000K. Fortunately, the opacity at these temperatures
is roughly constant, so we will assume a constant opacity κR = 3cm2 /g (Semenov et al., 2003).
Lastly, ρ at the mid-plane can be written as
1 Σ
.
ρ= √
2π h

(3.9)

Equations 3.1 through 3.9 form a system of 9 equations and 11 unknowns: cs , Tc , h, ν, Σ, Td , τ, κR , ρ, α, ṁ.
In the familiar problem solved by Shakura and Sunyaev (1973), α and ṁ are specified a-priori andthese equations have analytic solutions for the 9 remaining unknowns. We will extend this systemof equations in the next section to include a relationship between the net-magnetic field and the αparameter.

3.2.2

Coupling MRI Turbulence to the Net Magnetic Field

For a disk with a Hall-effect-dominated dead zone, the turbulence in the disk is coupled
to the net-vertical magnetic field. Non-ideal MHD simulations of the MRI have shown that the
space-time averaged α value in the case of a net field aligned with the spin axis is ∼ 0.01 or slightly
higher. However, in the case that the field is anti-aligned with the spin axis the MRI is suppressedsignificantly and α ∼ 10−4 to 10−3 (Simon et al., 2015b).
We will model this behaviour in a simple fashion, defining α in the Hall dead zone as the65following function of the magnetic field:




α

βz Ω > 0αdead ,βz Ω < 0,active ,α(βz ) =(3.10)

where βz is the usual plasma β parameter calculated using the net vertical field and including its
sign:
!
Pmid
,
βz = sign(Bz )Bz2(3.11)

where Pmid is the gas pressure at the mid-plane. Throughout this chapter we adopt αactive = 10−2
and αdead = 3 × 10−4 unless otherwise noted.
The above equation for determining α based on the magnetic field applies only in the Hall
MHD regime. This regime occurs when the central disk temperature is not high enough to sufficiently ionize the gas, which will be the case for temperatures below roughly 1000K (Umebayashi,1983). We will define the radius where the temperature is equal to 1000K as rdz , the start of theHall dead zone. We will use our model to solve for this location when each simulation is initialized
and assume it to be static as time evolves. A more advanced model might dynamically calculate
the location where Tc =1000K and change the transport velocities accordingly, but this is difficult
to do self consistently. In addition, temperature by nature does not change as drastically as otherquantities, usually coming into equations as T 4 , so assuming this to be static is not unreasonable.
A detailed discussion of the calculation of this radius will come later in section 3.5.
Adding equations 3.10 and 3.11 to the 9 equations from the previous section creates a systemof 11 equations and 13 unknowns (now including βz and Bz ). We will close this system in the next
two sections with time-evolution equations for Σ and Bz that are only a function of already existing
variables.

66
3.2.3

Surface Density and Net Magnetic Field Evolution

The surface density profile of the disk is evolved using the usual equation for a thin Kepleriandisk evolving due to viscous torques:∂Σ3 ∂ 1/2 ∂  1/2 r=r νΣ .
∂tr ∂r∂r(3.12)

We will follow Lubow et al. (1994) and Guilet and Ogilvie (2014) to time-evolve the poloidal
and radial components of the magnetic field in a turbulent disk. We will give an outline of themethod here, leaving the reader to refer to the original papers for some of the finer points and more
involved methods. The magnetic potential ψ is defined by
Bz =

1 ∂ψ
r ∂r

(3.13)

where Bz is the vertical component of the net magnetic field. The magnetic flux is transported by
the inward accretion flow (advection) and by outward by turbulent diffusion. Under the influenceof these two effects, the potential’s time evolution is given by∂ψ= −r(vadv Bz + vdiff Brs )∂t(3.14)

where Brs is the radial magnetic at the surface of the disk, vadv is the advective velocity, andvdiff is the diffusive velocity. Given ψ, Brs can be determined by solving an integral equation, or
equivalently solving a linear system of equations in the case of the discretized problem. For a full
discussion of these equations and there discretization, we refer the reader to Guilet and Ogilvie(2014), equations 6-9, 26, and 27.
Following the α model of turbulence, the advective velocity and time-scale are given by

∂ 
1/2
νΣr
,
Σr−1/2 ∂rr.

vadv = −tadv =

vadv

3(3.15)(3.16)The turbulent magnetic diffusivity η is defined by assuming a constant magnetic prandtl number
P r = νη (expected and assumed to be unity). In other words, we assume ν = η. The diffusive speed

67
and time-scale are then given byηη  h −1
=
,
h
r r
r
.
tdiff =
vdiff

vdiff =

(3.17)(3.18)In this case, the relevant dimensionless parameter that describes the ratio of viscous diffusion tomagnetic diffusion is the effective prandtl number: Preff . In the steady state, where vadv = ν/r,
this is given by
Preff =

vadv
3h
=
Pr.
vdiff
2r

(3.19)

Note that diffusion is faster than advection by a factor of (h/r)−1 (a primary result from Lubow
et al. (1994)). Thin disks diffuse net-flux outward, while only very thick disks can efficiently accrete
net-flux inward.

3.2.4

Evolving the Model Self-Consistently

We now have time evolution equations for the surface density and magnetic field that dependon the statically-solvable variables. Our time-evolution scheme is as follows:(1) Initialize by solving the static equations (3.1-3.11) at each radii with a constant initial ṁand power law Bz .
(2) Advance Σ and Bz using the time-evolution equations given in section 3.2.3.
(3) Re-solve the static equations at each radii, now with Σ known and ṁ unknown.
(4) Continue to repeat steps 2 and 3.
The numerical methods used to implement this algorithm are explained in detail in the next section.

3.2.5

Summary of Disk Parameters

While this is not the most sophisticated disk model, it’s simplicity means that there are
relatively few assumed parameters that need to be taken. In general, specifying our two primary

68
control parameters for the disk: α and ṁ, determines the disk’s properties. Our fiducial choices of
disk parameters that will be used throughout the chapter are found in table 3.1.

3.2.6

Stellar Parameters

While the focus of this work is on the physics of the disk, we must briefly discuss our choiceof parameters for the stars we are considering. Disks form at roughly the same time as their hoststars and typically only have lifetimes of a few million years, so we are considering very youngstars. As a general rule, young stars are larger in radius, more rapidly rotating, and more strongly
magnetized than their main sequence counterparts. However, the properties of young stars varyquite widely in these parameters, in particular the rotation rate and magnetic field strength. Inorder to constrain the scope of this work, we will err on the side of adopting closer to solar values
for each stellar type. Table 3.2 shows our assumed stellar parameters, loosely based on Johns-Krull(2009); Böhm-Vitense (2007); Vidotto et al. (2013); Casagrande et al. (2008). In addition to theparameters in the table, we will generally take the period of the stellar magnetic cycle to be a freeparameter, but our fiducial choice will be 20 years.
While we consider stars with low stellar masses in this work, we do not attempt to fully
address the many interesting and erratic phenomena associated with M-dwarfs. An interesting
topic for future study might be to integrate some sort of stochastic modeling of M-dwarf behaviour
into the model given here, but that is outside of the scope of this work. In a sense, we are considering”solar-like-M-dwarfs”: their mass, rotation period, radius, and field strength will be appropriate
for an M-dwarf, but their magnetic field cycles will be assumed to be similar to our sun.

3.3

Methods

3.3.1

Spatial Discretization

We use a logarithmically spaced grid in the radial direction that extends from rmin = 0.005
AU to rmax = 5 AU and has 600 points. This grid represents the maximum inner and outer extent

69

Parameterαidealαdz, active
αdz, dead
κR
µ
Pr

Value
10−2
10−2
3 × 10−4
3 cm2 /g
2.7
1

Table 3.1 Parameter choices for our disk model. Given these parameters, a specification of M∗(which determines Ω), and ṁ, our static disk model determines all other quantities as a functionof radius.

Class / Reference Name
G
K
M1
M4
M6

M∗ /M
1
0.8
0.5
0.3
0.1

R∗ /R
2
1.6
1.0
0.6
0.2

Table 3.2 Assumed stellar parameters for several stellar types.

B∗,0 (G)
1000
1000
2000
2000
2000

P∗
7 days12 days
2 days2 days2 days70
that we will use in any simulation. For each given simulation, a sub-section of this grid will beselected as the grid for that specific problem based on a calculation of the inner truncation radiusof the disk for the given parameters. Our choice of truncation radius as a function of the star anddisk parameters is given in the following section. The outer radius is not intended to be the outeredge of the disk, but is simply chosen to be far enough away from the relevant phenomena andcyclic behaviour that the outer boundary condition does not have an effect on the phenomena. The
diffusive time-scale at 5 AU for a disk with α = 0.1, h/r = 0.1 is on the order of several hundredyears. This is significantly longer than reasonable stellar cycle times (∼ years to tens of years), so
the net-field will not be directly coupled to the stellar cycle at this radius.

3.3.2

Boundary Conditions and Truncation Radius

Our boundary conditions (BCs) for the magnetic field will follow a similar prescription to
Guilet and Ogilvie (2014). The transport velocities follow the equations given in the previous
section between rin and rout . These radii define buffer zones at the inner and outer edges of the
disk: rmin < r < rin and rout < r < rmax respectively. We define a smoothing function that goes tozero at the grid boundaries and has a vanishing derivative on both ends:
smooth(x) =3x3
1 + 2x3(3.20)

where x = (r − rmin )/(rin − rmin ) between rmin and rin and x = (r − rmax )/(rout − rmax ) between
rout and rmax . This smoothing function multiplicatively modifies the diffusion velocity in the inner
buffer and the advection velocity in the outer buffer. Our choice of rin is physically motivated andis explained below. Our choice of rout is simply 75% of the rmax .
For the surface density, the inner BC is already taken care of, as it follows the same advection
velocity as the magnetic field and we want to allow matter to freely pass through our inner boundary.
However, the advection velocity for the surface density does not go smoothly to zero at the outerboundary as is does for the magnetic field. We use a feeding condition for the outer boundary, setsuch that ṁ stays constant at the outer-most grid point. As emphasized in the previous section,

71
this outer boundary is far enough away from the relevant phenomena that our choice of outerboundary condition does not affect our results.
Next we will consider the location of the truncation radius rin . The co-rotation radius rco
is defined to be where the stellar angular velocity (and therefore the angular velocity of stellar
magnetic field lines) is equal to the Keplerian angular velocity. Interior to this radius, the disk
orbits faster than the star and stellar field lines exert a retrograde torque on the gas, causing it tolose angular momentum. If this angular momentum loss is faster than the viscous evolution, thisis effectively where the disk is truncated. The magnetosphere radius rms , defines where the stellar
field is strong enough to generate this super-viscous magnetospheric accretion. These two radii ofinterest are given by
rms =

rco =

!2/7
3B∗2 R∗6,√
2Ṁ GM∗!1/3
GM∗ P∗2,4π 2(3.21)(3.22)

where M∗ is the stellar mass and P∗ the stellar rotation period (Koenigl, 1991). Canonical parameters for sun-like T Tauri stars (B∗ = 1kG, R∗ = 2R , Ṁ = 10−8 M yr−1 , M∗ = M , P∗ = 7 days)
yield rms ' 0.08AU and rco ' 0.07AU (Armitage, 2009). We will take the truncation radius to be
the lesser of these two radii. In table 3.2 we give the assumed stellar parameters that will be usedin this calculation for several types of stars that we will be consider later in the chapter.

3.3.3

Smoothing α

It is reasonable to assume that the turbulence is not perfectly local, and that it spreads fromone radii to another on a scale of roughly ∼ h. We use an smoothing procedure such that alpha isslightly non-local on this scale.
αsmooth,i = mean(αlocal,j , j = i − k, j = i + k)

(3.23)

where i and j are grid indices, and k is the number of grid cells that α is being smoothed over. We
take k=5, which is equivalent to a length scale of 0.025r for all r, given our spatial discretization

72
scheme. This length is on the order of the scale height, which will generally be between 0.01r and
0.1r.

3.3.4

Advancing in TimeThe model is evolved in time with a Runge-Kutta 4th order algorithm. Because there are two,
coupled, time-evolution equations, this requires that the full state of the disk (i.e. the staticallysolvable quantities) be solved for 4 times each time-step. The time step is dynamic and is calculated
by applying the CFL condition for the 2 transport terms:!
∆r 2∆rdt = C × min,,
vadv πvdiff(3.24)

where the safety factor, C, is chosen to be 0.1.
We now have an implementation of a self consistent model that co-evolves the surface density
and net-field of an accretion disk in 1 dimension. While this could be useful to study a numberof topics (for example the fate of primordial net-fields in disks), for the remainder of the chapter
we will seek to understand the influence of the magnetic field of a central star on the surrounding
accretion disk.

3.4

Modeling the Influence of the Stellar Magnetic Cycle on the disk
We will assume that a star in isolation will produce a dipole magnetic field given byB∗ (r) = B∗,0

r
R∗

−3
,

(3.25)

where B∗ is the vertical component of the dipole field in the equatorial plane, R∗ is the radius of
the star, B∗,0 is the field strength at the star’s surface, and r is the distance from the star in it’sequatorial plane. The strength of the field will vary with the stellar cycle, so we define the field
above to be when the field is at it’s strongest.
Next, we will consider the static magnetic field geometry in the case of a central star surrounded by a disk. We will assume that the disk is threaded by a net-field in the vertical direction

73(ie. open field lines). While the exact geometry will be complicated and depend on the relative
field strengths, we reason that there will be an inner “compressed-dipole” region that overlaps theinner edge of the disk, and an outer region with open field lines that will apply for most of the disk.
This is qualitatively similar to the geometries determined by Ostriker and Shu (1995), Muzerolle
et al. (2004), and Fendt et al. (1995). A qualitative representation of this field geometry is shownin figure 3.1.
Near the inner edge of the disk, closed lines will push the open lines threading the disk
outward, creating a small region where the disk will be threaded by a concentrated dipole field. At
these radii the disk will likely be thermally ionized throughout the whole column, but at the very
least would have ionized surface layers capable of carrying currents. Thus, we must now consider
the interaction of a stellar dipole magnetic field with a rotating, conducting accretion disk.
Lovelace et al. (1995) and Lynden-Bell and Boily (1994) have explored very similar problems
analytically: a stationary dipole interacting with conducting disk that is rotating as a solid-bodyand centered on the dipole. The magnetic field geometry that this situation produces is quite
complex, but we are primarily interested in the net flux contributed to the disk. Differential
rotation between the disk and star causes stellar field lines to twist. After the field lines are twisted√
more than 1/ 3 turns, they will reconnect and result in open field lines threading the disk. Because
this critical twist is small, we will assume that a closed field line threading the disk will become an
open field line after roughly one differential-rotation time. While there may be interesting effects
associated with these re-connection events that occur on the the time-scale of days (ie. the rotation
period of the star) (van Ballegooijen, 1994), we here we are consider the longer term impact of thisnet-flux contribution to the disk (time scale on the order of years). Thus, we assume that at the
inner edge of the disk closed stellar field lines will become open nearly instantly, and will contributenet-flux to the disk.
The stellar magnetic cycle causes the dipole strength (B∗,0 ) to vary in a roughly sinusoidal
way. We will model this by adding flux to the inner edge of the disk at a constant rate during thepositive-rising and negative-falling parts of the sinusoidal curve, with no flux being added during

74

Figure 3.1 A qualitative representation of the static magnetic field geometry. The black circle inthe center represents the star and the gray shaded area is the disk. The disk is threaded by open
field lines shown in red. The star produces closed dipole-like field lines, shown in blue. In isolation
this would be a true dipole, but the interaction with the open line region causes the dipole lines
to be compressed inward, while the open lines are pushed outward. This results in a region where
there are concentrated closed field lines from the star overlapping the inner edge of the disk.

75
the positive-falling and negative-rising quadrants of the cycle. The rate of flux added to the disk
is then given by
∂ΨB
4
=±
∂t
tcycleZ ∞2πrB∗ (r)dr,(3.26)

rin

where rin is the radius at which the disk is truncated and tcycle is the period of the stellar magnetic
cycle. This rate is defined such that the flux added to the inner edge of the disk during a quarterstellar-cycle is equal to the flux of the isolated stellar dipole threading radii from rin to infinity.

3.5

A Criterion for Outbursts from Hall Dead Zones
In this section we will summarize the outburst mechanism proposed in Armitage (2016) andconsider an analytic criterion for generating outbursts.

3.5.1

Outburst Mechanism

Armitage (2016) proposes a mechanism by which outbursts events can be triggered in adisk by the interaction of the stellar magnetic cycle with the disk’s Hall (DZ). This mechanism issummarized by figure 3.2.
Consider an initial configuration where the stellar and disk fields are anti-aligned with thespin-axis of the star. In this case, the Hall DZ is in the dead state: the MRI does not drive
turbulence to the same degree as ideal MHD or a Hall DZ with an aligned net-field. The turbulent
viscosity and advection velocity in the DZ are reduced by a factor of αactive /αdead compared tothe active state. Continuity demands that, in a steady state, ṁ must be constant as a function ofradius even if the α parameter is not constant. ṁ in the steady state is given simply by νΣ, andaccording to our α-prescription, ν ∝ α. Therefore, Σ in the DZ is enhanced by the inverse of thefactor that the transport speed is reduced: αdead /αactive .
At some point in the stellar magnetic cycle, the star’s field will change sign to be alignedwith the spin-axis of the star. At this point, the star will start contributing net-flux to the inneredge of the disk via the mechanism described in section 3.4. Once in the disk, this flux will diffuse
outward (Lubow et al., 1994) according to the diffusion velocity given in equation 3.17. When this76
wave-front of aligned flux is able to reach the Hall DZ, the portion of the DZ that ends up withan aligned net-field will transition from inactive to active, increasing the turbulence and transport
velocities. This causes the accretion rate in the active part of the DZ to increase by a factor of
αactive /αdead , and the excess mass that was stored in the slowly accreting dead zone is released andflows inward. This burst of released mass causes ṁ to increase throughout the inner disk, not justin the DZ region.

3.5.2

Outburst Criterion

In order for this mechanism to occur, the time it takes for the diffusion front to diffuse
through the inner disk and reach the Hall DZ must be on the order of the period of the stellar cycleor shorter (ie. tdiff (rdz ) . tcycle ). If tcycle  tdiff (rdz ) then the diffusion fronts of stellar flux will
combine together interior to the Hall DZ and the net-field in the DZ will effectively be decoupledfrom the stellar cycle.
While the statement of this condition for outbursts is relatively straightforward, there arequite a few parameters than need to be accounted for in calculating the diffusion time at rdz .
Combining equations 3.18 and 3.19 and again assuming Pr = 1 yields an expression for the diffusion
time in terms of the advection velocity:  
3
h −1
tdiff = r
v
2
r adv

(3.27)The expression for vadv that is used in the time-dependent model is given in equation 3.15, however
for the purposes of this simple analytic calculation we will use the steady-state approximation:
vadv = ν/r (eq. 3.4). We will again follow the α-prescription for viscosity: ν = αcs h, and will
substitute cs = hΩ (eq. 3.2). Combining these relations gives a concise expression for the diffusion
time:

−1

tdiff = Ω

 −1
h
α−1 ,
r

(3.28)

77

Figure 3.2 A cartoon illustration of the outburst mechanism proposed by Armitage (2016).

78
where all quantities are a function of radius except α. In this calculation we are considering field
diffusing through the active (ionized) region, so the appropriate α parameter will be αactive .
We will now reduce this expression further to be only a function of radius and our primarycontrol parameters: α, ṁ, and M∗ . For now we will ignore pre-factors and worry only about thepower law dependencies of each parameter, we will re-normalize the final result. Ω is defined as
usual:
1/2Ω ∝ M∗ r−3/2 .

(3.29)

We can use the static-disk model given in section 3.2.1 to solve for h/r for any given α and ṁ, andalso account for Ω’s dependence on M∗ :h
−7/20 1/5 −1/10∝ r1/20 M∗ṁ αr(3.30)

Plugging these scalings to equation 3.28 and re-normalizing gives
tdiff
=
308yr



α
10−2

−9/10 M∗M−3/20 ṁ
−7
10 M /yr−1/5 r1AU

29/20(3.31)

We need to evaluate the above expression at the radius where the disk transitions fromionized to neutral: rdz . This occurs when the central temperature of the disk (Tc ) drops below
about 1000K. Following the same procedure as for h/r above, the model given in section 3.2 canbe solved to determine the central temperature as a function of the control parameters, yielding
the following result:
Tc
=
767K



α
10−2

−1/5 M∗M

3/10 ṁ10−7 M /yr

2/5 

r1AU

−9/10(3.32)

Setting Tc = 1000K and solving for rdz givesrdz=
0.744AUα10−2

−2/9 M∗M

1/3 

ṁ
10−7 M /yr

4/9(3.33)

Finally, evaluating equation 3.31 at rdz yields the condition for outbursts as only a function of thecontrol parameters:−11/9  
4/9αM∗ 1/3ṁ
tdiff (rdz ) = 201yr. tcycle .
10−2
M10−7 M /yr(3.34)

79
From this relation we can see that outbursts will favor disks with high α, low ṁ, and stars
with low-mass and long stellar cycles. We will visualize this condition by plotting the power law
for tdiff (rdz ) as a function of ṁ, for an assumed α values of 0.01 and 0.03 We plot this for several
stellar masses. If tcycle lies above the power law for a given ṁ, then outbursts are likely. If tcycle
lies significantly below the power law, then outbursts are not likely.
Inspecting figure 3.3, we can see that for a sun-like solar cycle period, disks with α = 0.01
will only have outbursts for ṁ . 10−8 , even in the case of the smallest stars (M∗ ∼ 0.1M ). The
values plotted for ṁ represent fairly low, yet not completely unreasonable values for the accretion
rates of disks: 10−7 − 10−10 M /yr. In general the accretion rate is expected to decrease over the
lifetime of the disk. Consequently, based on this analysis, these sort of events might be expectedin the later stages of the evolution of the disk.

3.6

Results of Time-Dependent ModelingIn this section, we will finally apply our time-dependent model to the proposed outburst

mechanism. Our tests show that the stellar cycle’s contribution to the net field of the disk isgenerally dominant over any background field in the disk on the scales that trigger these outburst(tenths of an AU). As a result, the initialization of the background magnetic field is not a particularly
important factor. We choose to initialize the disk with anti-aligned flux following a power law:Bz ∝ r−1/2 (Okuzumi et al., 2014). The field strength is normalized such that βz = 105 at 5AU.
These choices give a fairly direct test of the qualitative mechanism proposed by Armitage (2016).
Furthermore, while this model does not always enter a perfectly cyclic and well behaved pattern
initially, we generally try to plot sections of each test that are typical of the overall behaviour.
First, we demonstrate the diffusion of the stellar field through the disk in cases which bothdo and do not produce outbursts. Figure 3.4 shows a space-time diagrams of the vertical magnetic
field in the disk, also denoting the radii at which the Hall DZ begins. We see that in the case ofa solar-like star with an accretion rate of ṁ = 10−8 M /yr, outbursts do not occur because thediffusion fronts of aligned flux do not reach the Hall DZ. In contrast, figure 3.5 shows an out-

80

Figure 3.3 The power laws for the diffusion time at the radius where the Hall DZ begins plottedfor several stellar masses and two α values. The α values here refer to α in the active region. The
The horizontal line shows the solar cycle period for the sun, roughly 20 years. For any given ṁ
and M∗ , if tcycle is above or roughly equal to the diffusion time, outbursts or cyclic behavior would
be expected.

81bursting cycle for an M4 star with a very low accretion rate of ṁ = 10−10 M /yr. In this case theHall DZ starts very near the inner truncation radius and the stellar flux is easily able to reach the
DZ on the time-scale set by the stellar cycle (20 years in this example).
In order to quantify the magnitude of the episodic accretion events, we calculate a simple
multi-temperature-blackbody luminosity of the disk. While this does not directly correlate to theapparent brightness of the star/disk system, it is nonetheless a useful measure of the relative changein accretion energy. The total luminosity of the disk at any point in time is given by:
Z routL=

2πrσTd4 ,(3.35)

rin

where Td is the surface temperature of the disk, which is self-consistently evolved in our model.
This luminosity is then normalized by the luminosity of the “steady-dead-state” in which the disk
is initialized. We plot an array of light-curves determined by this calculation, in figure 3.6, as wellas in the bottom panel of figure 3.5.
Naively, if the full Hall DZ was changed from dead to active, we would expect the accretion
rate (and hence luminosity) in the DZ to jump by a factor of αactive /αdead (∼ 30 for our parameters).
However, the inner (always active) portion of the disk will also exhibit an enhanced accretion rate
as the excess mass that was trapped in the DZ accretes inward. Consequently, the accretion
luminosity could, in principle, increase by more than this factor. In reality, the entire dead zone
does not convert to active at once, so we would expect αactive /αdead to be a rough upper limit on
the change in luminosity. The light curves we calculate generally confirm this naive result, showing
characteristic changes in accretion luminosity between a factor of a few and a factor of ∼ 30. For
outbursts that just barely meet the criterion from §3.5, the change in luminosity is smaller.
The time-scale of the outbursts is generally about a quarter of the stellar magnetic cycleor less, although the light curves are more complicated than a simple on/off paradigm. Weakoutbursts tend to also have shorter time-scales. Physically, this corresponds to the scenario where
the aligned-flux is just barely able to reach the inner edge of the Hall DZ before being cancelledout by the stellar magnetic field changing direction (see row 1 column 2 and row 3 column 2 in

82

Figure 3.4 A space-time diagram of the net-field in the disk for a G star with tcycle = 20yr andṁ = 10−8 M /yr. The horizontal dashed line represents the location of the inner edge of the Hall
DZ. In this case, the diffusion fronts of aligned flux do not reach the inner edge of the dead zone
before the stellar field changes sign, so no episodic or cyclic accretion occurs.

Figure 3.5 Space-time diagrams and the disk luminosity for a M4 star with tcycle = 20yr andṁ = 10−10 M /yr. The top 3 panels show net-field, local accretion rate, and surface density
respectively, with the bottom plot showing the luminosity calculated via equation 3.35.

83

Figure 3.6 Light curves for all 5 stellar types considered for an accretion rate of ṁ = 10−9 M /yr.
The stellar types decline in mass from top to bottom. The stellar cycle time increases from left toright: 6, 20, and 40 years respectively. Note that the scale of the x-axis is not the same in all plots.

84
figure 3.6).
These tests demonstrate that the simple time-scale criterion given in the previous section,
which implies that smaller stars are friendlier to these outbursts, is not all that needs to be considered. The criterion did not account for any stellar parameters other than mass, but the rotation
rate, field strength, and stellar radius also play an important role here in determining the truncation radius of the disk. K stars, which are physically larger and generally rotate more slowly thanM stars, will truncate their disks at a larger radius. This truncation radius can end up being quite
close to the start of the Hall DZ, giving the stellar flux a shorter distance to travel to trigger an
outburst. In figure 3.6, we see that K stars are actually the only class considered that produces
outbursts for tcycle = 6 years at the given accretion rate.

3.7

Conclusions and Discussion
We propose a model for a new mechanism for outbursts from protoplanetary disks. Thismechanism is a result of several features of these disks and their host stars:(1) Hall dead zones exist in the inner disk (∼ 0.1 AU to ∼ 1 AU) and suppress MRI turbulence
only if the vertical magnetic field is anti-aligned with the spin axis.
(2) The host star undergoes magnetic cycles in which the dipole-like field changes direction
periodically. This stellar field contributes flux to the very inner part of the disk (∼ 0.01−0.1
AU) significantly enough to dominate over any existing background field.
(3) Thin disks efficiently diffuse magnetic field outward, allowing this stellar field to reach theHall dead zone and influence the net-field of the disk on larger scales (∼ 0.1 − 1 AU).
We develop a 1D model to quantitatively describe the features of these outbursts and generate
expected light curves. Our primary conclusions are as follows:(1) On the scale of tenths of AU, the net field of the disk is dominated by the stellar magnetic
field and it’s cycle.

85(2) The mechanism is most effective for disks with low accretion rates and/or stars with longmagnetic cycles, with an accretion rate of 10−8 M /year or less required to get outbursts
from a star with a sun-like magnetic cycle in a disk with αideal = 10−2 . If stars were to
exhibit outbursts of this type, they would likely be class II sources, which have accreted
most of their mass already and have lower accretion rates.
(3) Outbursts due to this mechanism type are cyclic, assuming the stellar magnetic cycle iscyclic. Over the course of the cycle, the integrated luminosity changes by a factor of up to
αactive /αdead .
(4) The time-scale for outbursts is a quarter of the magnetic cycle time or less, with weakeroutbursts also being shorter. This is a useful, testable prediction.
This mechanism is able to produce events that roughly match some important features ofExors:(1) Exor outbursts have typical duration’s of 1-2 years (Herbig, 2008) with several years between events. Our model would predict this for a stellar cycle time of 5-10 years, which isshorter than the Sun’s but certainly still within reason.
(2) Exors typically have changes in luminosity of a few magnitudes, which agrees well with our
prediction of luminosity changes on the order of αactive /αdead .
(3) When observed in quiescence, Exor objects typically have absorption features typical of
K or M dwarfs and emission features typical of type II T-Tauri objects (Parsamian andMujica, 2004; Herbig, 2008). This agrees reasonably well with our prediction that thismechanism operates best around solar or sub-solar mass stars that are accreting slowly.
While these results are a useful proof of concept and are somewhat promising, the accretion
rates required by our outburst criteria are quite low, even for disks in the later stages of their evolution. Measured accretion rates for EXor objects are typically on the order of 10−7 M /yr (Audard

86
et al., 2014). In addition, there is more work to be done in considering a more detailed descriptionof the stellar field’s contribution to the disk’s net-field. While our prescription is reasonable, thisis a complicated magnetic field configuration that deserves further modeling. Furthermore, if there
is an Ohmic-dominated DZ interior to the Hall-dominated DZ, the low transport velocities in the
Ohmic DZ would act as a barrier to diffusion, making this mechanism less likely to operate. Finally,
the physics of young M-dwarfs merit more consideration than we have been able to give in thiswork. If the erratic behaviour of their magnetic activity was accounted for, we would generally
expect erratic accretion behaviour from the disk as well.

Chapter 4

The Nature of Turbulence in the Outer Regions of Protoplanetary Disks

4.1

Introduction
One of the most basic and fundamental processes at work in accretion disks is the removalof angular momentum from the orbiting gas, which allows this gas to accrete onto the central
object. A commonly accepted way for this transport to occur is turbulence, which allows adjacentannuli to exchange angular momentum on a macroscopic level, creating an effective enhanced
viscosity or “turbulent viscosity” (Shakura and Sunyaev, 1973). The leading candidate for sourcing
this turbulence is the magnetorotational instability (MRI; Balbus and Hawley, 1991, 1998), which
generates vigorous turbulence from an initially weak (compared to thermal pressure) and sufficiently
well-coupled (to the disk gas) magnetic field.
The degree of coupling between the magnetic field and the gas is crucial to the efficacy of the
MRI in transporting angular momentum. In fully ionized disks, such as those around black holes,
this coupling is very strong. However, the colder environments present in protoplanetary disks leadto low ionization levels; the result is substantially weakened gas-field coupling such that the MRI
cannot operate at full capacity throughout the disk. Early considerations of low ionization effects
on the MRI led to a model, put forth by Gammie (1996), in which large Ohmic diffusivity within
the very weakly ionized mid-plane quenches the MRI there, creating a so-called “dead zone”. At
larger heights from the mid-plane, however, the gas is ionized by cosmic-rays (Gammie, 1996),
X-rays (Igea and Glassgold, 1999), and FUV photons (Perez-Becker and Chiang, 2011a).
More recently, it has been realized that two other low-ionization, “non-ideal” magnetohydro-

88
dynamic (MHD) effects, namely, the Hall effect (which arises from ion-electron drift) and ambipolar
diffusion (AD; which arises from ion-neutral drift), are also important in protoplanetary disks (Wardle, 1999; Sano and Stone, 2002; Fromang et al., 2002; Bai and Goodman, 2009; Bai and Stone,2011). The relative strength of all three non-ideal MHD effects is dependent primarily on the density and ionization of the gas and thus the radius away from the central star (Kunz and Balbus,
2004; Desch, 2004). Broadly speaking, the Ohmic term dominates at high densities, the Hall term
dominates at intermediate densities, and the ambipolar term dominates at low densities; thus,
moving away from the star, the dominant non-ideal term transitions from Ohmic diffusion to theHall effect to ambipolar diffusion.
Considerable theoretical work has been carried out to study the role that both the Hall
effect and ambipolar diffusion play in determining the gas dynamics in protoplanetary disks. For
ambipolar diffusion, in particular, analytic work has shown that the relevant control parameter (atleast in the linear limit) for the growth of the MRI is the ratio of the rate of ion-neutral collisions
to the orbital frequency Ω, defined to be Am. If Am . 1, AD damps the MRI, although unstablemodes will always exists for a sufficiently weak field (Blaes and Balbus, 1994; Kunz and Balbus,
2004). This damping has been the topic of several numerical studies as well, seeking to understandthe full, non-linear evolution of the MRI under these conditions. Hawley and Stone (1998) studied
this problem treating the ions and neutrals as separate fluids coupled with a collisional drag term.
They establish limits on the relevance of AD, showing that the fluid behaves essentially the same as
the ideal limit for collision frequencies on the order of 100Ω, while for collision frequencies . 0.01Ω
the neutral fluid is quiescent, unable to be affected by the motion of the ions. Bai and Stone (2011)
study the non-linear evolution of the ambiploar MRI in an unstratified shearing box using a single
fluid method. In contrast to the two-fluid approach, in the single fluid MHD limit, the ion-neutral
coupling is strong (which is shown to be appropriate for protoplanetary disks Bai 2011a) the ion
density is a fixed fraction of the neutral density instead of obeying it’s own continuity equation,and the ion-density makes up a negligible fraction of the total mass-density. Bai and Stone (2011)
showed that for any given Am value, one can always find a magnetic field strength at which the89
MRI is present in its full nonlinear state. For smaller Am, weaker magnetic field is required, which
in turn leads to a lower level of turbulence.
With these results in hand, Bai (2011a,b) use complex chemical networks in combination
with ionization models for FUV, X-rays, and cosmic rays to calculate the ionization fraction andnon-ideal MHD parameters in the disk as a function of z. They then combine these calculations
with the α − βmin − Am relationship from their unstratified simulations to predict the turbulenceand stress in the disk as a function of z, and estimate the accretion rate. Modeling of this nature
implicitly assumes that each layer of the disk operates more or less independently, an assumptionthat remains untested. Indeed, several effects are relevant in suggesting that processes cannot beheight-independent, including vertical turbulent diffusion, explicit diffusion of magnetic fields due
to non-ideal effects, buoyancy, and advective transport due to large scale velocity modes that havebeen described analytically (Lubow and Pringle, 1993; Ogilvie, 1998) and have been seen in someshearing box experiments (Oishi and Mac Low, 2009; Gole et al., 2016). The combination of theseeffects may play a significant role in determining the steady state structure of the outer regions ofdisks and may predict a significantly different accretion rate.
The first such work to self-consistently model the AD-dominated outer disk with verticalgravity included was carried out by Simon et al. (2013b) and Simon et al. (2013a). These authors
carried out local and vertically stratified calculations, in which AD was quite strong near the midplane. This work included ionization of upper disk layers by FUV photons (Perez-Becker andChiang, 2011a) leading to a scenario similar to the Gammie (1996) model; active regions vertically
sandwich a region of weak or no turbulence. Furthermore, these works determined that a net
vertical field is required to predict high enough accretion rates to match the expected life-time ofdisks around T-Tauri stars (assuming that disk accretion is at least approximately in steady-state).
The theoretical work done to-date has established a new understanding of turbulence in theouter regions of protoplanetary disks. However, as alluded to above, the precise nature of thisturbulence, and in particular, the role of vertical transport of energy (in any form) has not yet beenaddressed in detail. It is a primary motivation of this chapter to better understand turbulence in

90
the AD-dominated outer regions of protoplanetary disks, and in particular to address whether or
not each vertical layer can be treated (largely) independently.
Additionally, beyond purely theoretical considerations, it has now become possible to put direct observational constraints on turbulence in the outer regions of disks with ALMA. Predictionsof the strength of MRI-driven turbulence in an AD-dominated disk (Simon et al., 2015a) have been
directly tested with observations of molecular line broadening (Flaherty et al., 2015, 2017, 2018) as
well as dust settling (Okuzumi et al., 2016; Pinte et al., 2016). These observations have (surprisingly) suggested that turbulence in these regions is quite weak. These considerations motivate amore precise understanding of the physics of the outer disk regions for two primary reasons. First,
the presence of weak turbulence has only been established for a few systems; if disks with strong
turbulence do exist, it behooves us to explore this turbulence in more detail. Second, if indeed all
disks show significantly weak turbulence, then we need to understand precisely why the MRI is
prevented from operating.
Thus, both theoretical and observational considerations strongly motivate us to deepen our
understanding of turbulence in the outer regions of protoplanetary disks. In this work, we carryout a series of numerical simulations to do just that.
This chapter is organized as follows: in §4.2 we describe the numerical algorithm used in our
simulations and explain the physics and initial conditions of the simulations, in §4.3 we present theresults from control simulations with no FUV layers, in section §4.4 we demonstrate the gas andmagnetic field behaviour in simulations with FUV layers and contrast them with results from §4.3,
in §4.5 we do a temporal correlation analysis to characterize the nature of the MRI and verticalmagnetic flux transport in this situation, and in §4.6 we summarize and discuss our primary results.

4.2

Numerical Algorithm, Parameters, and Initial Conditions
We carry out a series of numerical simulations with Athena, which is a second-order accurate

Godunov flux-conservative code that uses constrained transport (CT; Evans and Hawley, 1988a)
to enforce the ∇ · B = 0 constraint, the third-order in space piecewise parabolic method (PPM) of

91
Colella and Woodward (1984) for spatial reconstruction, and the HLLD Riemann solver (Miyoshi
and Kusano, 2005; Mignone, 2007) to calculate numerical fluxes. A full description of the Athenaalgorithm along with results showing the code’s performance on various test problems can be foundin Gardiner and Stone (2005), Gardiner and Stone (2008), and Stone et al. (2008).
We employ the shearing box approximation (Hawley et al., 1995) to simulate a small, corotating patch of an accretion disk. Taking the size of the shearing box to be small relative to the
distance from the central object, we define Cartesian coordinates (x, y, z) in terms of the cylindricalcoordinates (R, φ, z 0 ) such that x = R − R0 , y = R0 φ, and z = z 0 . The box co-rotates about thecentral object with angular velocity Ω, corresponding to the Keplerian angular velocity at R0 . The
equations of MHD, with ambipolar diffusion as the only non-ideal MHD term, in the shearing box
approximation are,∂ρ+ ∇ · (ρv) = 0,∂t∂(ρv)
+ ∇ · (ρvv − BB + ¯ĪP *) = 2qρΩ2 xî − ρΩ2 z k̂ − 2Ωk̂ × ρv,∂t(J × B) × B ∂B−∇× v×B−= 0.
∂tγc ρi ρ∂E
+ ∇ · (E + P *)v − B(B · v) = 0
∂t
J=∇×B
E=

ρ(v · v) B · B
P
+
+
γ−1
2
2

(4.1a)(4.1b)(4.1c)(4.1d)(4.1e)(4.1f)

Here ρ is the gas density, ρi the ion density, γc is the coefficient of momentum transfer in ion-neutral
collisions, B is the magnetic field, J is the current density, P * is the total pressure (related to themagnetic field and the gas pressure, P , via P * = P + B · B/2), ¯Ī is the identity tensor, v is thevelocity, E is the total energy, and q is the shear parameter defined by q = −d ln Ω/d ln R, whichis 3/2 in this case of a Keplerian disk.
We use an adiabatic equation of state with a constant background temperature profile. The
disk is heated and cooled towards this profile with a characteristic cooling/heating time τcool . This

92
adds the following term to the energy evolution equation:∂Eρ(T − T0 )=−∂tτcool (γ − 1)

(4.2)

where T ≡ P/ρ and T0 is the constant background temperature. We set τcool = 0.1Ω−1 .
We use the standard boundary conditions appropriate for the shearing box approximation.
At the azimuthal boundaries we exploit the symmetry of the disk and apply periodic boundary
conditions. At the radial boundaries, we apply shearing periodic boundaries; quantities that are
remapped from one radial boundary to the other are shifted in the azimuthal direction by the
distance the boundaries have sheared relative to each other at the given time (Hawley et al., 1995).
Furthermore, the energy and azimuthal velocity are appropriately adjusted to account for thedifference in angular momentum and orbital energy (see Hawley et al., 1995). Additionally, theseradial shearing-periodic boundaries are incorporated in such a way as to maintain conservationof conserved variables (see Stone and Gardiner, 2010). In the vertical direction, we employ themodified outflow boundary condition of Simon et al. (2013b), which are standard outflow boundaries
with the gas density extrapolated exponentially into the ghost zones; this method reduces theartificial build-up of toroidal magnetic flux near the vertical boundaries. Such outflow boundaries
can lead to significant mass loss, particularly in the case of strong vertical magnetic flux (Simon
et al., 2013a). In order to keep the box in a steady state, we re-normalize the mass of the box atevery time step. Finally, we employ Crank-Nicholson differencing to conserve epicyclic energy tomachine precision. A full discussion of these methods and Athena’s shearing box algorithm canbe found in Stone and Gardiner (2010).
Our standard choice of domain size is Lx × Ly × Lz = 2H × 4H × 8H, where H is theisothermal scale height, defined to beH=
where cs is calculated using T0 .1

√ cs
2 .
Ω

(4.3)

We also carry out simulations with 4H × 8H × 8H to test

convergence with domain size (see table 4.1). Our simulations without an FUV layer will only1

While our simulations are not isothermal, the cooling timescale is sufficiently short that departures away fromisothermality are minimal; see below

93
extend 5H in the vertical direction, as the low density regions at large |z| will impose a very strict
constraint on the time-step from the strong diffusion in those regions.
We initialize the box with net-magnetic-flux in the vertical direction such that the βz,0 ,defined as
βz,0 =

2Pmid
2 ,
Bz,0

(4.4)

is βz,0 = 104 . Note that a factor of 4π has been subsumed into the definition of the magnetic
pressure, following the standard definition of units in Athena (Stone et al., 2008). In the caseof a constant vertical field, the onset of the MRI creates a very strong channel mode that cancreate numerical problems (Miller and Stone, 2000). To avoid this, the initial field has a spatial
dependence in the x direction:Bz = Bz,0
1
1 + sin(kx ∗ x) ,
4

(4.5)

where kx = 2π/Lx , and Lx is the length of the box in the radial direction.
In order to seed the MRI, we apply random perturbations to the initial gas pressure andvelocity. Following the procedure used by Hawley et al. (1995), the perturbations are uniformly
distributed throughout the box and have zero mean. Pressure perturbations are a maximum of
2.5% of the local pressure, and velocity perturbations are a maximum of 5 × 10−3 cs .
We initialize the disk is in hydrostatic equilibrium. In an vertically isothermal disk, balancing
gravity against the gas pressure gradient yields a Gaussian density profile. While the disk is not
isothermal in this case, it is heated/cooled towards an isothermal profile. The heating/cooling timeis shorter than the dynamical time, meaning that temperature perturbations will be restored to theisothermal temperature on a faster time-scale than the disk’s density profile will respond on thescale of H. For this reason, we will consider H to be constant in our initialization and discussion.
We set ρ0 = 1 as the initial density at the mid-plane. A density floor of 10−4 is applied to keep
the Alfvén speed from getting too high (and thus the time step too small) in the region far fromthe mid-plane. The density floor also prevents β from getting too low, which can cause numerical√problems. We set the constant, background isothermal sound speed cs = 1/ 2, and choose the94
angular velocity to be Ω = 1, which gives H = 1.

4.2.1

Resistivity Profile

For the FUV-ionized layers of the disk we use a resistivity model informed by the considerations of Perez-Becker and Chiang (2011b) and adapted from other works (e.g., Simon et al., 2013a).
FUV radiation will penetrate some constant column density of the disk given by ΣFUV . We will callthe point at which this column density occurs zFUV . As the box evolves, turbulence and magnetic
pressure can be significant enough to change the density profile. As a result, we need to re-calculate
zFUV at each time-step. This is done by integrating the horizontally averaged density profile fromthe top of the box towards the mid-plane until ΣF U V is reached. The same calculation is made forthe bottom half of the box separately.
We will describe the strength of ambipolar diffusion in terms of the ambipolar Elsasser number, Am, defined as
Am =

γρi
.
Ω

(4.6)

Am in the active layer (|z| > zFUV ) is generally  1. We will use the prescription given in Bai andStone (2013a):7

Am = 3.37 × 10



f
10−5



ρ



ρmid

r
1AU

−5/4
.

(4.7)

We take the ionization fraction f to be 10−5 and place the box at r = 115 AU in a minimum masssolar nebula model (Weidenschilling, 1977; Hayashi, 1981).
At the mid-plane the ionization fraction is generally controlled by X-rays and cosmic rays,
which are able to penetrate farther into the disk. Bai (2011b) model this ionization using a complex
chemical reaction network and show that Am in the mid-plane region (|z| < zFUV ) is generally
constant as a function of space and is on the order of 1. Including small grains in this model
decreases Am and gives it a dependence on β. We will not include this more complicated behaviour,
but will consider a range of values for Ammid to account for this uncertainty. The ionization profile
will be smoothed over a thin region to prevent potential numerical problems, following the same

95
procedure as Simon et al. (2013b). Our choice of ΣFUV is 0.1g/cm2 .

4.2.2

Resolution Considerations

Our fiducial choice of resolution for these simulations is 30 grid cells per H. This choice isbased on the size-scales on which the MRI operates in the linear limit. This is typically characterized
by the critical wavelength λc and the most unstable wavelength λm . All size-scales larger than λc
are unstable, and λm corresponds to the mode withf the largest growth-rate. For ideal MHD,
−1/2

λm /H = 9.18β0, while λc is roughly half that (see Hawley et al., 1995). However, ambipolar

diffusion increases both of these wavelengths significantly when Am . 10 (Wardle, 1999; Bai andStone, 2011). We will consider a variety of Am parameters for the mid-plane region, spanning 10−2
to 102 , while Am will always follow equation 4.7 for |z| > zFUV . Figure 4.1 shows the most unstable
and critical wavelengths at the mid-plane and at z/H = 1.5 for our choice of β0 = 104 . The MRI
is reasonably well resolved for z > zFUV , with 9 grid cells per λm at z/H = 1.5. However, at the
mid-plane our resolution only yields 3 grid cells per λm for Am & 1 . While this is not ideal,
previous convergence studies have shown that non-linear state of the MRI is not greatly changed
even if, in the linear limit, it is under-resolved (Simon et al., 2013a). This is potentially due totwo effects. First, even in the linear limit there are larger unstable modes than λm that will bebetter resolved. Second, as the MRI evolves into the non-linear regime the field will generally getstronger, increasing the most unstable size-scale such that it is better resolved.
We should also note that at very high diffusivities, λc can exceed the size of the ambipolar
damping zone, which is ≈ 2H, for our fiducial choice of parameters. For example, at Am = 0.01,
λc ∼ 4H for β0 = 104 . If no MRI modes fit within the damping zone, then this is another way inwhich ambipolar diffusion can effectively quench the MRI, in addition to that described in moredetail below. However, it is also worth mentioning that this argument only applies to vertical MRI
modes and not toroidal or radial modes, which may still persist within the damping region. We
will also conduct simulations with no FUV layers, for which the vertical extent of the box is 5H.
In these simulations, our most diffusive case will be Am = 0.032, and at this diffusivity both λc

96
and λm fit within the vertical extent of the domain.

4.2.3

List of Simulations and Varied Parameters

Our main control parameter for this investigation will be Ammid . In order to rule out othereffects and help confirm our results, we will run a few simulations at higher resolution, and a few
with a larger domain in the horizontal directions. In addition to the simulations with the FUV
layer, we will perform simulations with constant Am across the domain as a comparison. In table
4.1 we list all of our simulations. Our labeling convention is as follows: The first letter is S or L
for small or large box size. The second letter is L or H for low or high resolution. The next letter
will be F or N for simulations with FUV or no FUV respectively. The number following Am is theAm parameter at the mid-plane.

4.3

Stratified Ambipolar Simulations with no FUV
In this section we will consider a series of vertically stratified shearing boxes with no FUV-

ionized layers (SLN Am0.032 through SLN Am10). These will serve as a baseline for comparison
with the FUV simulations.

4.3.1

Expectations

A simple way to consider the effect of resistive conditions on the MRI is to compare the
MRI growth time-scale (∼ Ω−1 ) with the diffusion time-scale of the magnetic field. If the diffusion
time-scale is shorter than the growth time-scale, field will diffuse faster than the MRI can build it
up and the MRI will not operate. In the case of Ohmic diffusion, this is relatively straightforward
because the Ohmic diffusion coefficient is not a function of the field strength. So it is possible
to say whether the MRI will operate or not based simply on the Elsasser number (which is itself
essentially just a ratio of these time-scales). However, ambipolar diffusion is more involved. The
ambipolar term in the induction equations expands into many terms (see equation 4.1c). Justfor the sake of a heuristic argument, let’s just consider the term that looks the most like simple

97

Figure 4.1 The critical and most unstable MRI wavelengths at the mid-plane and at z/H = 1.5) interms of our fiducial cell size δx = H/30. β0 = 104 at the mid-plane, which corresponds to β ∼ 103
at z/H = 1.5 due to the decrease in gas pressure.

98

Label
SLF Am0.01
SLF Am0.018
SLF Am0.032
SLN Am0.032
SLF Am0.056
SLF Am0.1
SLN Am0.1
LLF Am0.1
SHF Am0.1
SLF Am0.18
SLF Am0.32
SLN Am0.32
SLF Am0.56
SLF Am1
SLN Am1
SHF Am1
LLF Am1
SLF Am1.8
SLF Am3.2
SLN Am3.2
SLF Am5.6
SLF Am10
SLN Am10
SLF Am18
SLF Am32
SLF Am56
SLF Am100

Ammid
0.01
0.018
0.032
0.032
0.056
0.1
0.1
0.1
0.1
0.18
0.32
0.32
0.56
1.0
1.0
1.0
1.0
1.8
3.2
3.2
5.6
10
10
18
32
56
100

Table 4.1 Simulation Parameters

FUV
yes
yes
yes
no
yes
yes
no
yes
yes
yes
yes
no
yes
yes
no
yes
yes
yes
yes
no
yes
yes
no
yes
yes
yes
yes

H/∆x
30
30
30
30
30
30
30
30
60
30
30
30
30
30
30
60
30
30
30
30
30
30
30
30
30
30
30

(Lx, Ly, Lz)/H
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
4,8,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
4,8,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8
2,4,8

Orbits
100
50
100
75
200
400
300
100
25
400
400
300
400
400
300
200
400
400
400
300
400
100
300
100
100
100
100

99
diffusion, where the time derivative of a component is proportional to a double spatial derivativeof the same component:


∂B
∂t


=
AD

B2
· ∇2 B.
ΩρAm

(4.8)

This consideration is particularly valid in a shearing box for the By component, which is expected
to dominate over the other field components due to the shear in the disk. Therefore this term with
By2 as a pre-factor will dominate over any cross-terms. Notice that in this diffusion approximation,
the diffusion coefficient itself scales with B 2 , so the diffusion time-scale for length scale L does aswell:
L2 ΩρAm
=
tAD =
B2



L
H

2

β Am
,
Ω

(4.9)

with the second equivalency being true for vertically isothermal disk with constant sound speed,a reasonable approximation for our disk with an isothermal equilibrium temperature and a shortcooling/heating time. Considering a disk with a relatively weak seed-field, the MRI will be able tooperate initially while the field is weak because the diffusion is also weak. However, as the MRI
winds-up the field, a saturation field strength will be reached where the diffusion time-scale andgrowth time-scales are now equal and the MRI no longer operates. Setting tAD = Ω−1 in the above
equation yields the simple relation that the saturation β parameter will be ∝ Am−1 .
The fully non-linear version of this problem was considered by Bai and Stone (2011) in a
3D, unstratified, fully ambipolar diffusion dominated shearing box. They fit the results of theirsimulations to show that the minimum β parameter that permits the MRI is given by2 
2 1/250
8βmin (Am) =
+
+1
.
Am1.2
Am0.3

(4.10)The simple expectation is roughly recovered for the limit Am . 1, where βmin ∝ Am−1.2 .
4.3.2

Magnetic Field Strength

We characterize the magnetic field strength in terms of the β parameter corresponding to thetotal field:
β=2P
Bx2 + By2 + Bz2.

(4.11)

100
When space-time averages of β are calculated, the pressure and total field averages are taken first,
and then hβi is calculated from the ratio, i.e. hβi = h2P i/hBi. For the MRI operating in a shearing
box, we expect By to be the dominant field component as radial field will be sheared into toroidal
field by the shear velocity profile. In this case, B ' By and β ' βy .
In figure 4.2 we show the space-time averaged β parameter in the mid-plane region of ourstratified disk simulations. We find that the average field strength in the saturated state closely
follows the βmin − Am relationship described by Bai and Stone (2011) for their unstratified simulations.
We need to be cautious in using the average field strength, as it is possible that there exist
regions of very strong field separated from weak field regions within the mid-plane (as shown by
Simon et al. 2017a in similar calculations); such non-uniform structure can lead to an average β
value that is not representative of the typical field strength found at the mid-plane. It would thus
be possible to sustain some turbulence sourced from small patches of weak-field (β > βmin ) even
if on average β < βmin . We calculate the distribution of mid-plane field strength, integrated overtime and for |z| < 0.5H. Examples of these distributions are shown in figure 4.3. Especially forthe stronger diffusion cases, there are often regions where β < βmin , but the field strength does
not exceed this limit drastically and βmin is close to the center of each distribution. So while βmin
does not seem to be a strict limit on the field strength, it is generally gives a good description of
the saturated state field strength in the stratified case in addition to the unstratified case.

4.3.3

Turbulence and StressThroughout this chapter we will consider perturbations to various quantities, defined (for aquantity a) as
δa(x, y, z, t) ≡ a(x, y, z, t) − haixy (z, t).

(4.12)

where the xy subscript on the angled bracket denotes a horizontal average. For vector quantities,
the perturbation of each component is calculated and then the total perturbation is the sum of the
3 components in quadrature.

101

Figure 4.2 Temporally and spatially averaged β within 0.5H of the mid-plane. The black line shows
βmin (see Bai and Stone, 2011). βy generally follows βmin .

Figure 4.3 β distributions for simulations with Am = 0.1, 1, 10 in blue, green, and red respectively.
The dashed vertical lines corresponds to βmin for each simulation. In general, βmin gives a reasonably
good description of the typical field strength.

102
The first set of turbulence diagnostics that we calculare are the perturbed ρ and v. In
figure 4.4, we show that both of these perturbations scale as a power law with Am. The strengthof turbulent velocity perturbations (relative to the gas sound speed) in particular has been ofmuch interest recently as it provides a direct comparison to observations of turbulent line widths(Flaherty et al., 2015, 2017, 2018). Many of these observations probe optically thick gas at large |z|(usually |z| ∼ 2–3H) in the outer regions of protoplanetary disk. To date, these observations haveconstrained turbulence to be quite weak (δv/cs < 0.05) for two systems: HD163296 (Flaherty et al.,2015, 2017) and TW Hya (Flaherty et al., 2018). In an attempt to explain this weak turbulence,
Simon et al. (2017a) proposes a scenario in which the outer disk is shielded from ionizing sources
and thus strong ambipolar diffusion, analogous to what we are simulating here.
Given these considerations, we plot the the vertical profiles for the perturbed velocity in
Fig. 4.5, finding (in agreement with Simon et al. 2017a) that at large heights above the mid-plane,
δv/cs approaches values above the observational constraint. Even in our most diffusive simulation(Am = 0.032), turbulence at large |z| is above that of HD163296 and approaches the limit fromTW Hya.
While these simulations do not extend into the corona far enough to cover the entire observed
region, we expect that δv would continue to increase as z increases. More specifically, we speculate
that, in the Am = 0.032 simulation, δv would exceed or equal the limit for TW Hya if we were toconsider the region between z/H = 2 and 3. So even in the scenario where a disk has no ionizedlayers, an even lower ionization (i.e., ambipolar diffusion; Am < 0.032) or a weaker backgroundvertical field may be required to be consistent with observations.
Finally, we examine the Maxwell stress, and we calculate this stress in two different ways:
background-subtracted field perturbations, −hδBx δBy i, and the average field components, −hBx ihBy i,
evaluated for each constant-z plane. These perturbation and laminar contributions are plotted infigure 4.6. The stress scales as a power law with Am, with perturbations providing the dominantcomponent for Am < 1.

103

Figure 4.4 Hydrodynamic diagnostics of turbulence δρ and δv time and space averaged (within
0.5H of the mid-plane) for simulations SLN Am0.032 through SLN Am10. The level of turbulence
increases as the damping decreases, consistent with expectations.

Figure 4.5 The velocity perturbation profiles for our series of simulations with no FUV layer(constant Am). The profiles are shown in a gradient from blue to red, corresponding to
Am = 0.032, 0.1, 0.32, 1.0, 3.2, and 10.0. The observed limits on TW Hya and HD163296 (Flaherty et al., 2015, 2017, 2018) are shown in the dotted and dashed lines respectively. (Note that in
the case of TW Hya, this limit only applies to |z| ∼ 2–3H.) Even in a disk with no FUV-ionizedlayers, for β0 = 104 it would take a diffusivity of at least Am = 0.032, and likely higher, to explain
the observations of weak turbulence.

104

Figure 4.6 The magnetic stress calculated from perturbations (black) and from laminar fields (blue).
Note that at Am = 1, the points are directly on top of one another. The stress clearly scales with
Am. For Am < 1, the stress is dominated by small-scale perturbations.

105

4.4

Stratified Ambipolar Simulations with FUV
Having established a baseline of simulations, and accompanying diagnostics, with no FUV

ionization, in this section, we analyze simulations with FUV ionization (SLF Am0.01 throughSLF Am100) and consider the behaviour of the gas and magnetic field at the mid-plane as Ammid
is varied. We define the mid-plane, above-transition, and corona regions to be |z/H| = 0 to 0.5,
1.0 to 1.5, and 3 to 3.5 respectively. The simulations generally saturate in terms of turbulence andfield strength within about 10 orbits, so the time-averaged quantities shown in the plots in this andfollowing sections are taken starting at 10 orbits and extending to the end of the run. This gives
us a baseline of 40 orbits at minimum and over 100 orbits for most runs (see table 4.1).

4.4.1

Mid-plane field strength

We calculate the time and space averaged β in various regions of the box for our simulations
with an FUV ionized active layer. We find that the mid-plane β is less than the expected minimum
described in the previous section for simulations with Ammid < 1 (Figure 4.7).
As described briefly in §4.3, it is possible to have a situation where the average magnetic fieldexceeds the maximum that permits the AD-dominated MRI while still having patches of weaker
field where the MRI is active. This idea has been seen recently in similar simulations by Simon
et al. (2017a). Thus, simply calculating the average β parameter is not sufficient, and so we also
calculate the distribution of mid-plane β, integrated over time and for |z| < 0.5H. Examples ofthis distribution are shown in figure 4.8. We find that the fraction of cells for which β > βmin ison the order of a few percent for Am < 1. So while there may be small patches for which the MRI
is permitted, this likely does not dominate the dynamics of the mid-plane. We also see that the
distribution for the Am = 1 case is qualitatively different than the Am < 1 cases, and that the
Am < 1 cases are similar to one another.
These results have two implications:(1) The magnetic field at the mid-plane is not produced by the local MRI, as the field strength

106

Figure 4.7 The time-averaged β parameter plotted versus Am for SLF Am0.01 S0.1 throughSLF Am100 S0.1 for 3 spatial regions: the mid-plane, above the transition, and the corona. The
black line shows the minimum β for the mid-plane MRI (given in equation 4.10). Because the fieldstrength exceeds what would be expected from the local MRI for Am < 1, we expect that this field
is sourced elsewhere and transported to the mid-plane.

Figure 4.8 Distribution of β taken over the mid-plane region of SLF Am0.032 through SLF Am1.0,
with the mid-plane Am shown on each plot. The vertical dashed lines indicate the minimum betathat permits the MRI for each Am parameter. While the MRI might occasionally be active insmall patches, it is unlikely to dominate the dynamics of the mid-plane when Am < 1.

107
exceeds what the local, AD-dominated MRI would be able to maintain for the given Am
parameter. The field here must be transported from the MRI active regions above andbelow the mid-plane.
(2) In the majority of the mid-plane region (> 90%), the MRI cannot be active because thefield strength makes ambipolar diffusion strong enough to damp field perturbations fasterthan they can grow.

4.4.2

Hydrodynamic Turbulent Properties

We again calculate the time and space averaged perturbations of density and velocity; in figure
4.9 we see that hydrodynamic diagnostics of turbulence are independent of the magnetic diffusivityof the mid-plane for Am < 1. This agrees with the expectation from the previous section, as
the field is too strong for the MRI to be operating locally and generating density and velocity
perturbations. Any turbulence at the mid-plane has trickled down from the corona (possibily
related to the mechanisms seen in Fleming and Stone, 2003; Oishi and Mac Low, 2009; Gole et al.,2016, for Ohmic dead zones), and given that the properties of the corona are not varied in thesesimulations, one would expect to recover a constant perturbation amplitude at the mid-plane as afunction of Am for Am < 1. For Am > 1 we see a power law emerge, similar to the case with no
FUV layer where we see a power law across the entire range of Am tested. This is indicative that
the MRI is operating in this regime and gets stronger as Am increases, as expected. We should
also note that the properties of the corona are not effected by the diffusivity of the mid-plane.

4.4.3

Stress and Magnetic field perturbations

In Fig. 4.10, we calculate the average field perturbation as a function of Am for different
vertical regions. The mid-plane field perturbations scale roughly as a power law of Am with anindex of 1/3, even for Am < 1. This is in contrast to the hydrodynamic diagnostics of turbulence
discussed in the previous section that did not depend on Am in this regime. Considering a simple
diffusion equation for the magnetic field (e.g., equation 8), field perturbations will be diffused away

108

Figure 4.9 Velocity and density perturbations for SLF Am0.01 through SLF Am100. The midplane quantities do not scale with Am for Am < 1 indicating that these perturbations are notgenerated by local MRI. In addition, the mid-plane turbulence does not seem to affect the level of
turbulence in the corona.

109
in direct proportion to the magnetic diffusivity. Thus, we see that the field perturbations scalewith Am. Hydrodynamic perturbations, on the other hand, have no explicit diffusion, but canbe damped out by turbulent motions (e.g., through a turbulent cascade to the grid scale, where
numerical diffusion dominates); thus, these particular quantities do not scale with Am for Am < 1.

The expected α parameter in the ambipolar region determined by Bai and Stone (2011) isgiven byα = 1/(2hβi).

(4.13)

If we assume β will saturate to βmin , we can combine this relation with equation 4.10 to get an
expected relationship between Am and the stress for the AD-dominated MRI. Note that this is anupper limit on this stress, as β will not always reach βmin , depending on the field geometry (Bai andStone 2011 figure 15). We do not necessarily expect to recover this relation for Am < 1 because
the MRI is not active, but will use it as a point of comparison.
We again calculate the Maxwell stress from both laminar field and field perturbations, shownfor the mid-plane region in figure 4.11. Given that each individual field component perturbation
scales as Am1/3 , we expect and do see that the stress due to perturbations scales roughly as Am2/3 ,being the multiplication of two field components. We find that the laminar stress dominates forhighly diffusive mid-planes, equaling the perturbed stress around Am = 1. This is qualitatively
different from the results of our simulations with no FUV layers (figure 4.6), where we see that thestress from magnetic field perturbations remains dominant for all simulations with Am < 1.
For Am < 1, there are two effects that make the observed mid-plane stress different than theprediction in equation 4.13 :(1) As shown in figure 4.7, hBy i is stronger than the prediction based on βmin . This will
enhance the laminar magnetic stress defined to be -hBx ihBy i. Figure 4.11 shows that thisis the dominant effect on the total stress for Am . 0.1, where the laminar field is able toproduce a stress of α ∼ 10−3 .

110Figure 4.10 Magnetic field perturbations versus Am within different vertical regions. The mid-plane
perturbations scale roughly as Am1/3 , indicated by the black dashed line. Despite the mid-planescaling, the mid-plane does not influence the other regions.

Figure 4.11 The magnetic stress calculated from perturbations (black) and from laminar fields(blue). The red line indicates the predicted upper limit of stress from Bai and Stone (2011). The
black dashed line indicates the expected power law scaling of the stress from perturbations: Am2/3
for Am < 1 (based on the scaling seen in figure 4.10), and Am1.2 for 1 ∼ Am (based on equations
4.10 and 4.13). The mid-plane stress is dominated by the laminar component for very diffusive
mid-planes (Am < 1).

111(2) The MRI is not active and is not generating perturbations locally. Any perturbations will
just be due to trickle down turbulence from the corona. This causes the stress due toperturbations to be less than the prediction. This is the dominant effect on the total stressfor 0.1 . Am . 1, where the total stress is significantly less than the upper limit.
For 1 . Am, we expect β to scale roughly as Am−1.2 until Am >> 1, where the MRI will
saturate and is essentially in the ideal MHD limit (equation 4.10). Based on equation 4.13, thismeans the stress should scale as Am1.2 . We recover this scaling between Am = 1 and ∼ 10, at
which point the stress saturates, remaining constant as Am increases.
We should note that, for completeness, we performed tests at higher resolution and with alarger box size in the horizontal directions to ensure that these would not affect our results (seetable 4.1). In all of these tests we find that the behaviour is consistent with our fiducial choice ofresolution and box size.

4.5

Vertical transport of magnetic field in the FUV simulations
Figure 4.12 shows the space-time diagram of the horizontally averaged toroidal field By forrun SLF Am0.1. Upon visual inspection of the diagram, it appears that even for a AD damped
mid-plane, the MRI dynamo starts at the mid-plane and then is transported upwards via buoyancy.
If this is true, it is inconsistent with our claim in the previous section that the MRI-active corona
is able to transport flux to the mid-plane. In this section, we will perform a temporal correlation
analysis to further elucidate this mystery.
We compute the temporal correlation of the mean By for each x-y plane with that at the
mid-plane:c(z, t0 ) =

Z T
t=0

By x,y (z = 0, t) ∗ By x,y (z, t − t0 )dt,

(4.14)

After computing this correlation function, the time-shift (t0 ) at which c is maximized is found foreach z. In our definition, t0 < 0 at a given z corresponds to magnetic field at that z lagging thatof the mid-plane. Thus, in the example of buoyantly rising magnetic fields, t0 < 0 and |t0 | would

112
increase as the field rose further away from the mid-plane.
As shown in Fig. 4.13, while field well above the transition region (denoted by the verticaldashed lines) has t0 < 0, suggesting that field is buoyantly rising away from the transition region, themid-plane region has t0 > 0 suggesting that the mid-plane itself lags behind the transition region.
This result holds for all of our simulations with Am < 10. The dominant region for controlling the
MRI dynamo is right at the damping zone transition. The butterfly diagram still has its typicalappearance because the corona still lags behind both the transition region and the mid-plane.
These correlation results are somewhat surprising, as it is generally difficult to transport
magnetic flux towards the mid-plane of a disk. Usually the dominant effect in vertical transportof flux is magnetic buoyancy, which gives space-time diagrams like figure 4.12 their characteristic“butterfly” appearance. Considering a magnetic flux tube with field strength Bpert in a disk,
the magnetic flux provides additional pressure such that the gas pressure required to maintain2hydrostatic balance against the vertical component of the central’s stars gravity is reduced by Bpert(in our units). In the isothermal case, the gas density will also be reduced within the flux tube by
2 /c2 . The flux tube is then lighter than it’s surroundings and experiences a force away fromBperts

the mid-plane. This effect is known as magnetic buoyancy. The buoyant force Fb on the gas in theflux tube is equal to the weight of the displaced fluid:
Fb =

2
Bpert
g,
c2s

(4.15)

where g is the acceleration due to vertical gravity. Dividing by the gas density in the tube andusing g = Ω2 z yields a simple expression for the buoyancy acceleration ab :
ab =

2
Bpert
Ω2 z
2
Ω
z
=
,
ρc2s
βpert

(4.16)

where βpert is the β parameter corresponding to Bpert . In order to arrive at a buoyant time-scale
tb , we consider the time it would take for a perturbation of constant β to rise H at this rate ofacceleration. This yields−1tb = Ω



2H
βpert
z

1/2
.

(4.17)

113

Figure 4.12 Space-time diagram of sign(By )By2 /2P0 for SLF Am0.1. Here, P0 is the mid-planepressure. The horizontal dashed lines indicate zFUV .

Figure 4.13 Time shift of maximum correlation relative to the mid-plane as a function of z for the
SLF Am0.1 simulation. The dashed vertical lines indicate the center of the transition between FUV
ionization and no FUV ionization. We see that the mid-plane lags behind the transition region,
implying that diffusion towards the mid-plane is able to beat out the effect of buoyancy.

114
While the flux tube is accelerating upwards it will also be diffusing due to ambipolar diffusion.
In order to transport a significant portion of the flux to the mid-plane, this diffusive effect will haveto occur on roughly the same time-scale as buoyancy or shorter. The ambipolar diffusion time
under some simplifying assumptions is given in equation 4.9. Over length-scale H, the diffusion
time is tAD = 2βpert AmΩ−1 . Taking the ratio of these two time-scales yields:td=tb 1/2z
1/2βpert Am.
h(4.18)√
Note that h = H/ 2. In other words, diffusion will be the dominant transport mechanism if
−1/2. If we consider an MRI-generated perturbation of
Am < Amcrit , where Amcrit = hz βpert
β = 3 at z/H = 1, Amcrit ∼ 0.5. This is roughly consistent with the results of our simulations
with FUV layers, specifically figure 4.8, which shows that the mid-plane magnetic field strength
distribution changes in nature somewhere between Am = 0.32 and 1. This likely corresponds to atransition from diffusion dominated transport to buoyancy dominated transport.
If ΣFUV is decreased, perturbations will start at larger |z/H| and we expect the effect of thecorona on the mid-plane will lessen due to two effects:(1) The mid-plane is farther from the initial perturbation, and the time-scale for diffusivetransport scales as length squared.
(2) Buoyant acceleration is ∝ z/H, so perturbations will rise faster from their starting position(this is accounted for in the relations given above).

4.6

Conclusions
In order to further understand the properties of turbulence in the outer regions of protoplan-etary disks, we have run a series of stratified shearing box simulations with varying strength ofambipolar diffusion and with and without an FUV layer. Our main conclusions are as follows:(1) Our simulations without an FUV layer generally match the behaviour of the similar but
unstratified simulations of Bai and Stone (2011). The β parameter at the mid-plane has a

115
narrow distribution centered near βmin , and both hydrodynamic and magnetic diagnostics
of turbulence scale with Am. This is consistent with the AD-damped MRI operating.
(2) In our simulations with an FUV layer and Am < 1 within the mid-plane region, we havedemonstrated that the behaviour of the mid-plane, even qualitatively, does not match the
behaviour of unstratified or stratified simulations without an FUV layer. This is due tothe interaction with the MRI-active layers, which can transport magnetic flux to the midplane to a point such that the MRI is no longer permitted. There are several importantdifferences compared to the no-FUV-layer case:• The magnetic stress is dominated by laminar fields rather than small scale field perturbations.
• Hydrodynamic perturbations do not scale with Am at the mid-plane; they are set by
the fluid motions in the FUV layers rather than local MRI.
• For most of the box, β < βmin. The field strength distribution is wider than the
no-FUV case and does not vary greatly as a function of Am.
All of these points are consistent with the MRI being arrested in most of the mid-plane
region by the magnetic field. The dynamics of the mid-plane are primarily determined by
the FUV layers instead of the local MRI. While the MRI may operate in small patches of
weak field, this does not seem to be a dominant effect and occurs in a small fraction of thebox, if at all.
(3) For Am . 10, toroidal field at the mid-plane is temporally correlated with toroidal field in
the corona with a lag time on the order of an orbit. This correlation implies that diffusion(turbulent or ambipolar) of flux towards the mid-plane is able to overcome magnetic buoyancy, which is usually the dominant effect in the vertical transport of magnetic flux in a
disk.
(4) The properties of the MRI-active region are not affected by the diffusivity of the mid-plane116
region.
What do these results imply for studies of turbulence in protoplanetary disks? First andforemost, if an FUV layer is present, there is a significant influence of the strongly ionized layers
on the highly diffusive mid-plane. Thus, one cannot in general treat each vertical layer of the disk
independently; analytic or semi-analytic models that attempt to do so will potentially be missing
physical effects (e.g., vertical diffusion of magnetic field) that will change estimates for the degreeof turbulence and angular momentum transport. In the absence of an FUV layer, this issue is lessclear however. While there will be no interaction with the mid-plane from an FUV layer for thesesimulations, buoyancy still likely plays an important role, which implies that one cannot treat eachlayer independently, even without an FUV layer.
Furthermore, our results (particularly the second and fourth conclusions above) are quite
relevant to recent observational (Flaherty et al., 2015, 2017, 2018) and theoretical (Simon et al.,2017a) results that suggest a picture in which turbulence in the outer regions of protoplanetarydisks is quite weak. If indeed, FUV photons are somehow prevented from reaching the outer disk,
as suggested by Simon et al. (2017a), then the turbulent velocity in the disk is largely set by thestrength of the ambipolar-diffusion-dominated MRI. Considering those simulations without FUV,
δv/cs ∼ 0.01 within the mid-plane region for the most diffusive case and rises up to δv/cs ∼ 0.05 at|z| ∼ 2.5H, as shown in Fig. 4.5. Given the general trend of increasing δv with |z|, it seems likelythat at larger |z|, the turbulent velocity will be marginally consistent with observations (whichshow δv/cs ∼ 0.05–0.1; Flaherty et al. 2017, 2018) at best, though potentially still larger thanobservational constraints. In either case, this result indicates that even for very diffusive disks(which can arise from small grains soaking up charges, or significant reduction of ionizing flux), tobe consistent with observations, one needs to reduce the background magnetic field strength (Baiand Stone, 2011; Simon et al., 2017a).
In summary, our results provide a more in-depth understanding of turbulence and angularmomentum transport in the outer regions of protoplanetary disks. Not only will this aid in future

117
modeling of these disks, but with observations of such systems continuing well into the next decade,such work is crucial to making progress in comparing theoretical expectations with observationaldata.

Chapter 5

Forming Planetesimals in a Disk with Driven Turbulence

5.1

Introduction
Perhaps the biggest problem in planet formation theory is the intermediate regime known

planetesimal formation, where particles must grow from cm size-scales to km size-scales. At micronsizes, simple coagulation is fairly efficient at growing grains to larger sizes (Dominik et al., 2007;
Dullemond and Dominik, 2005; Schmitt et al., 1997). Collisions between grains are frequent andthese collisions are more likely to stick grains together than break them apart. This process stops
being efficient around cm size-scales, at which point a particles become more likely to fragment upon
collision Zsom et al. (2010); Birnstiel et al. (2011). A new mechanism is needed to grow particles
to larger than cm-scales, and there is a strict time-constraint on this next growth mechanism dueto the radial drift of particles.
Radial drift is a consequence of the sub-Keplerian rotation of the gaseous component of the
disk. The gas has an outward-declining pressure gradient that gives some support against theradial force of gravity. This causes the centrifugal force required to balance gravity to be less thanthat of an object in orbit, leading to a slower rotation velocity for the gas than orbiting objects.
As a result, all solids in the disk experience a “headwind”, trying to move faster than the gas
they are embedded in. This headwind causes solids to loose angular momentum due to gas dragand drift inward towards the star, a phenomenon known as “radial drift” (Weidenschilling, 1977;
Takeuchi and Lin, 2002). The speed of the drift depends on the size of the particle considered, withthe highest speeds achieved for particles of roughly 1m in size. At this scale, particles will drift

119
inward to the central star in hundreds to thousands of years. This places a strict constraint on thetime-scale over which planetesimal formation process must occur.
A promising mechanism for rapid planetesimal formation is the streaming instability (SI).
This instability, first derived by Youdin and Goodman (2005), is able to concentrate solids that are
aerodynamically coupled to a fluid into streams and clumps. Considerable numerical work has been
done to confirm the SI (Youdin and Lithwick, 2007; Johansen and Youdin, 2007; Bai and Stone,
2010b) and quantify it’s properties (Johansen et al., 2009; Bai and Stone, 2010c; Carrera et al.,2015). While the properties of the SI depend on several physical parameters including the solid-togas ratio, the friction time for the particles, and the radial pressure gradient of the gas (more on thisin §5.2.3), studies generally find that the SI is able to concentrate particles into clumps with overdensities of 100-1000 times the background particle density. Very importantly, this concentration
happens on a timescale only an order of magnitude or so longer than the dynamical time Ω−1 ,which is faster than radial drift. If the density of the clumps formed is sufficiently high, self-gravity
will take over and collapse the clump down to a gravitationally bound object roughly the size of aplanetesimal. This combination of the streaming instability and self-gravity (SI/SG) has become
the leading candidate for the mechanism of planetesimal formation.
A focus of recent work has been to measure the planetesimal mass spectrum resulting fromdN
the SI/SG mechanism. This mass spectrum is often represented as a simple power law: dM∝ M −p ,where the power law index p is generally found to be in a range between 1.4 and 1.7. This power law
index is, somewhat surprisingly, largely invariant with respect to the friction time of the particles(Simon et al., 2017b), the radial pressure gradient (Abod et al., 2018), and the strength of selfgravity (Simon et al., 2016). These results pose the question: can any effects significantly change
the power law index of this spectrum?
In this chapter, we will explore the effect of a turbulent environment on planetesimal formation via the SI/SG mechanism. Disks are generally thought to be a turbulent environment, withviscosity arising from this turbulence playing the dominant role in angular momentum transport
within the disk (Shakura and Sunyaev, 1973). This study will not deal with the physical origin of

120
the turbulence, instead choosing to explicitly drive the fluid at an appropriate energy spectrum.
Our investigation has 2 primary goals: 1. Determine the threshold at which the turbulent velocity
dispersion is too high for the SI/SG mechanism be effective. 2. Determine the effect of turbulenceon the mass-spectrum of planetesimals.
This chapter is organized as follows: In §5.2 we describe the numerical algorithm, setup, andparameter choices for our simulations. In §5.3 we describe the tools and statistical methods wewill use to analyze the planetesimals. §5.4 we measure the effect of Kolmogorov turbulence on
the ability of the SI to form gravitationally-bound clumps. In §5.5, we propose and implement
a new paradigm for measuring planetesimal masses in these types of simulations, and use it tomeasure the planetesimal mass spectrum. Finally, in §5.6 we summarize our results and discuss
their implications.

5.2

Methods, Setup, and Testing
In this section we will describe the methods used to perform our coupled particle/gas shearing-box simulations, define and justify our choices of parameters, and test our scheme for driving
turbulence.

5.2.1

Fluid Algorithm and Parameters

In general, our setup will be very similar to that of recent papers (Simon et al., 2016; Abodet al., 2018; Bai and Stone, 2010b). We will evolve the equations of hydrodynamics (HD) using
Athena. Athena is a second-order accurate Godunov flux-conservative code that uses thirdorder in space piecewise parabolic method (PPM) of Colella and Woodward (1984) for spatialreconstruction and the HLLD Riemann solver (Miyoshi and Kusano, 2005; Mignone, 2007) to
calculate numerical fluxes. A full description of the Athena algorithm along with results showing
the code’s performance on various test problems can be found in Gardiner and Stone (2005),
Gardiner and Stone (2008), and Stone et al. (2008).
We employ the shearing box approximation to simulate a small, co-rotating patch of an121
accretion disk (Hawley et al., 1995). Taking the size of the shearing box to be small relative to the
distance from the central object, we define Cartesian coordinates (x, y, z) in terms of the cylindricalcoordinates (R, φ, z 0 ) such that x = R − R0 , y = R0 φ, and z = z 0 . The box co-rotates about thecentral object with angular velocity Ω, corresponding to the Keplerian angular velocity at R0 .
The equations governing the gas in these coupled gas-particle simulations in the shearing box
approximation are∂ρg+ ∇ · (ρg vg ) = 0,
∂t

 

ρp vp − vg∂(ρvg )22¯.
+ ∇ · (ρvg vg + ĪP ) = ρg 2qΩ xî − Ω z k̂ − 2Ωk̂ × vg +
∂t
ρg
tstop

(5.1a)(5.1b)

Here ρg is the gas density, P is the pressure, ¯Ī is the identity tensor, vg is gas the velocity, the shear
parameter defined by q = −d ln Ω/d ln R, which is 3/2 in this case of a Keplerian disk. The lastterm in the second equation is the particle feedback term, where ρp is the particle density, vp is theaverage particle velocity averaged over a fluid element (i.e. a grid cell), and tstop is the stopping
time for the particles. We will discuss the equations and parameters governing the evolution of
the particles in the next section. Finally, to close the system of fluid equations, we assume an
isothermal equation of state with a constant sound speed cs :P = ρg c2s .

(5.2)

We use a standard set of boundary conditions appropriate for the shearing box approximation.
At the azimuthal boundaries we exploit the symmetry of the disk and apply periodic boundary
conditions. At the radial boundaries we apply shearing periodic boundaries, such that quantities
that are remapped from one radial boundary to the other are shifted in the azimuthal direction by
the distance the boundaries have sheared relative to each other at the given time (Hawley et al.,1995). In the vertical direction, we employ the modified outflow boundary condition of Simon
et al. (2013b), which are the standard outflow boundaries but with the gas density extrapolated
exponentially into the ghost zones. Such outflow boundaries can lead to significant mass loss. To
compensate for this, the density at every point in the box is multiplied by a factor at each time-stepsuch that the total mass of the box remains constant. In general, these boundary conditions can

122
break conservation, so methods described in Stone and Gardiner (2010) are used in order to keep
the correct variables conserved. Furthermore, we employ Crank-Nicholson differencing to conserve
epicyclic energy to machine precision. A full discussion of these methods and Athena’s shearing
box algorithm can be found in Stone and Gardiner (2010).
We run all simulations using a domain size of Lx × Ly × Lz = 0.2h × 0.2h × 0.2h, where h is
the scale height, defined here as
h=

cs
.
Ω

(5.3)

We assume that the disk is initially in hydrostatic equilibrium. Balancing the vertical componentof gravity from the star with the gas pressure gives a Gaussian density profile:ρ(z) = ρ0 exp

−z 2
2h2


.

(5.4)

We set ρ0 = 1 as the initial density at the mid-plane. We set the constant isothermal sound speed
cs = 1, and choose the angular velocity to be Ω = 1, which gives h = 1.

5.2.2

Particle Algorithm

We co-evolve particles along side the gas using Athena’s hybrid gas-particle scheme as
described in Bai and Stone (2010a). The evolution of the each particle (indexed by i) is given by
0 )∂(vp,i∂t00= 2(vp,i,y− ηvkep )Ωî − (2 − q)vp,i,xΩĵ − Ω2 z k̂ −0 − v0
vp,i
g

tstop

+ Fsg .

(5.5)

where the prime indicates a velocity calculated in a frame where the background shear has been
subtracted, vkep is the Keplerian velocity, η is the factor by which the background gas velocity
differs from the Keplerian velocity, tstop is the stopping time, and Fsg is the force due to selfgravity between particles.
The self-gravity term is calculated in two steps. First, the gravitational potential due toparticle self-gravity (ψp ) is calculated by solving Poisson’s equation:∇2 ψp = 4πGρp .

(5.6)

123
The Poisson equation is solved following the methods described in the appendix of Koyama andOstriker (2009). The force on any given particle can then simply be written as
Fsg = −∇ψp .

(5.7)

It is important to note that this method for self-gravity prohibits it from acting on distances smaller
than a grid cell. In other words, despite the fact that particles can have any position and are not
restricted to the grid, a clump of particles that collapses due to it’s own self-gravity will not collapse
smaller than about the size of a grid cell. This limitation plays a very important role in interpreting
the results of our study and will be discussed further in the following sections.

5.2.3

Particle Parameters

The behaviour of the streaming instability (SI) and gravitational collapse of particle clumps
is uniquely defined by 4 dimensionless quantities (Youdin and Goodman, 2005; Simon et al., 2016):
the dimensionless stopping time τs , the solid-to-gas ratio Z, the radial pressure gradient parametere In this section we will define each ofΠ, and the dimensionless particle-self-gravity parameter G.
these quantities and briefly justify our choice for the parameter. As we wish to be able to compare
our results to previous studies, we will generally adopt similar parameters to these other recent
work.
The dimensionless stopping time is defined asτs = tstop Ω(5.8)

and represents the time-scale for the velocity of a moving object to be significantly affected byaerodynamic interactions with the gas. Massive objects like planetesimals are not well coupled
to the gas and will have long stopping times (τs  1). Small objects like dust particles arenearly perfectly coupled to the gas and have very short stopping times (τs  1). We choose anintermediate value τs = 0.3 for all particles in our simulations. While in principle it is possible tohave particles with a distribution of stopping times, using the same value for all particles simplifies
the interpretation of our results.

124
The solid-to-gas ratio (or metallicity) is defined as
Z=

Σp
,
Σg

(5.9)

where Σp and Σg correspond to the particle and gas surface densities respectively. We take Z = 0.02,
about twice the solar value. Depending on the value of the other parameters, a low metallicity canprohibit the growth of the SI completely (Bai and Stone, 2010a). We choose a value that permits
the SI in test simulations while still being well within reason.
Protoplanetary disks have an outward-decreasing pressure profile. This pressure gradientprovides a small, radially-inward force on the fluid, but does not act on particles. This additionalradial force on fluid elements causes the fluid of the disk to rotate at sub-Keplerian velocities. As
a result, particles in a disk experience a “headwind”. This effect is usually parameterized by theparameter Π as follows:
Π=η

vkep
.
cs

(5.10)

We take Π = 0.05, a fairly intermediate value relative to previous studies which should be in areasonable range based on most disk models (Abod et al., 2018; Bai and Stone, 2010b).
e denotes the relative strength of particle-self-gravity as compared to tidal shear:
Finally, G
e = 4πGρ0 .
GΩ2(5.11)

e corresponds to changing the gas density
The interpretation here is that changing the value of G
(and particle density along with it, assuming constant Z), or changing Ω (ie. changing the radiuse to be 0.1.
within the disk). We take G
These parameters lead to a natural choice for a mass scale: the “Gravitational Mass”, given
by

√
MG =

2 9/2 3 e
π Z Gρ0 h3 .
2

(5.12)

Physically, this corresponds to the mass enclosed by a circle with surface density Σp and a diameter
equal to the critical Toomre wavelength for gravitational collapse. Unless otherwise stated, allmasses throughout the paper will be in units of MG .

125
5.2.4

Stirring Turbulence: Spectrum of Velocity Perturbations

In order to simplify this problem as much as possible, we will be agnostic to the sourceof the turbulence in our simulations and will directly perturb the velocity field of the fluid toprovide turbulence. In order to calculate the velocity perturbations to be applied, we must consider
the desired power spectrum of the fluid’s velocity field. In traditional ”Kolmogorov” turbulence(Kolmogorov, 1941), the power spectrum of the resulting velocity field scales as a power law withthe spatial wavenumber k as follows:
E1 (k) ∝ k −p , p = 5/3
where k is the magnitude of the wavenumber (i.e. |k| =(5.13)q
kx2 + ky2 + kz2 in 3D Cartesian coordinates)

and E1 (k) is an linear energy density in k-space such that E1 (k)dk has units of energy.
Our goal is to set up a turbulent stirring scheme such that we can control the power law index
p of the power spectrum. In reality, we can only control the spectrum of the velocity perturbations;the resulting power spectrum of the fluid will be some convolution of the seed perturbations andthe turbulent cascade. We will implement a scheme similar to the one used by Dubinski et al.
(1995). Our plan will be to calculate the amplitude of each mode in 3D k-space (kx , ky , kz ), assign
each mode a random phase between 0 and 2π, then take the inverse transform of these amplitudes
to find the velocity perturbations to be applied to the fluid. The amplitudes in k-space will begiven byA(kx , ky , kz ) = (rand)(|k|)−b ,(5.14)

where A is the amplitude of the mode and ”rand” is a random number drawn from a normal
distribution. The power law index b can be related to p via dimensional arguments. Starting withE1 (k) ∝ k −p

(5.15)and multiplying both sides by dk yields an expression with units of energy:
E1 (k)dk ∝ k −p dk.

(5.16)

126
Considering we are assigning velocity amplitudes in 3D k-space, we need to convert this to be avolume density (E3 ) instead of a linear density:
E1 (k)dk ∝ E3 (kx , ky , kz )dkx dky dkz ∝ E3 (k)(4πk 2 )dk.

(5.17)

By inspection, the power law index for E1 must be 2 greater than the power law index for E3 :E3 (k) ∝ k −p−2 .

(5.18)

This is very close to the quantity we want to set the amplitude of the perturbations, but we aresetting them in velocity space, not energy space:p
E3 (k) ∝ A(kx , ky , kz ).

(5.19)

where A is the amplitude of a velocity mode in 3D k-space, which is by definition a volume-density
in k-space. Finally, combining the two previous proportionalities yields
A ∝ k −b , b =

p
+ 1.
2

(5.20)

So if we want E1 (k) ∝ k −5/3 , then we need to set the velocity modes with amplitude A ∝ k −11/6 .
We assign these amplitudes in a 3D lattice for all (kx , ky , kz ) such that |k| is between klow
and khigh . We define klow to be the wavenumber such that exactly one mode fits in the box andkhigh to correspond to the Nyquist wavenumber for the given resolution: khigh = (Nx /2)klow .
After these amplitudes are defined in k-space, the inverse Fourier transform is taken to
calculate velocity perturbations with the desired spectrum. The total energy of these perturbations
is then scaled such that the energy injection rate throughout the simulation will be a constant:dE
dt drive .
5.2.5

Stirring Turbulence: Magnitude of Velocity Perturbations

The energy injection rate mentioned in the previous section is our direct control parameter,
but we really want to be able to specify the amount of turbulence in terms of the α parameter

127(Shakura and Sunyaev, 1973). Throughout this chapter, we will very simply calculate this effective
α parameter as follows:
α=

X

2
h|δvi (x, y, z, t)|i ,(5.21)

i=x,y,zwhere
δvi (x, y, z, t) ≡ vi (x, y, z, t) − hvi ixy (z, t).

(5.22)The brackets indicate a spatial average, with the xy subscript denoting an average over the x-yplane. The over-bar indicates an average over time. Note that we are working in units such thatcs = 1, so the velocities in the equations above are already non-dimensionalized. If this was not
the case, all velocities would need to be normalized by cs .
In order to drive our simulations at given effective α parameter, we need to measure therelationship between the energy injection rate (the parameter we control directly) and the observed
δv of the box. In general, this energy injection rate is in balance with energy lost due to numericaldissipation. In this simple case the energy injection rate - δv scaling would match the naive guess2
based on an energy argument: dE
dt drive ∝ (δv) . However, we also expect to see outflows throughthe top and bottom boundaries of the box (Bai and Stone, 2013a; Ogilvie, 2012; Li et al., 2018),
which add another sink for kinetic energy. To complicate further, we observe that the strength ofthe outflow scales with δv. Thus, runs with more turbulence have proportionally more energy loss,
introducing a slight non-linearity into this relation.
We adopt a rather brute-force, yet effective, method to understand the energy injection rate
- δv scaling over the region of parameter space we wish to study. We perform a series of tests at low
resolution in which energy injection rates are varied by factors of 10 over 5 orders of magnitude.
Once each simulation reaches a saturated state, the average δv is calculated. The resulting relationis shown in figure 5.1. As demonstrated in the figure, we create a relation between the energy
injection rate and the resulting velocity perturbation by linearly interpolating between these 5 datapoints in log-space. This interpolation tells us the appropriate dE
dt drive parameter in order toachieve the desired δv (ie. α) in the saturated state of our high resolution simulations.

128

Figure 5.1 The relationship between the driving amplitude dE/dt and the resulting average velocity
perturbation δv. Each black point is generated by averaging over a low-resolution simulation.
The grey lines show an linear interpolation between the black points in log-space. The red lines
demonstrate how the appropriate driving rate is found to produce a simulation with δv = 10−2 , ie.
α = 10−4 .

129
5.2.6

Characteristics of Driven Turbulence

Before considering simulations with particles, we will characterize the properties of our boxof driven turbulence to ensure that our driving algorithm meets several requirements:• The turbulence is isotropic. The space-time-averaged perturbations of the 3 componentsof the velocity are the same.
• The turbulence is homogeneous. The time-averaged velocity perturbations are the same
throughout the box.
• The turbulence has the desired spectrum. The velocity perturbations calculated using themethods described in §5.2.4 must translate to the intended turbulent power spectrum inthe shearing box.
It is important that we do these tests in the absence of particles for several reasons. First, it
is not obvious that the arguments leading to the Kolmogorov power spectrum hold with particles
feeding back on the gas. In addition, particles tend to settle to the mid-plane, making the physics
of the gas inherently inhomogeneous. Finally, the streaming instability will itself drive turbulence
to some degree, and it is not necessarily true that this turbulence will have a Kolmogorov spectrum.
Figure 5.2 shows a vertical profile of each component of the velocity perturbation averaged
over the horizontal directions and from t = 5Ω−1 to 10Ω−1 . While we acknowledge that this is afairly short time-average, the box saturates with turbulence quite quickly with this driving schemeand the majority of our computing time was reserved for simulations with particles. The profile
for each component is relatively flat, confirming that our driving scheme produces homogeneous
turbulence. The 3 components of the velocity are also roughly equal throughout the box, confirming
isotropy, with the possible exception of δvy . We find it plausible that this slight anisotropy in they-direction (∼ 10%) is due to the shear velocity in the the box, which is also in y-direction andintroduces an inherently anisotropic effect.
We will now calculate the perturbed kinetic energy power spectrum of the fluid. We will

130

Figure 5.2 A demonstration that our stirring algorithm produces a box of homogeneous, isotropic
turbulence. Data is taken from simulations hr moderate np. Time averages are taken from 5Ω−1
to 10Ω−1 . The dashed horizontal line shows the expected δv given the strength of the driving.

131
follow the methods described in Simon et al. (2009). The Fourier transform of a quantity f (x) isrepresented by fe(k) and is defined as usual:Z Z Z
fe(k) =

f (x)e−ik·x d3 x.

(5.23)The kinetic energy density power spectrum is then calculated as follows:
1
1 X √^ 2
2
PS ρ|δv| (k) =| ρδvi (k)|2
2

(5.24)

i=x,y,zIn order to arrive at a spectrum that represents a linear density in k-space, we integrate the powerover spherical shells of constant |k|:




Z
1
1
2
2
PS ρ|δv| (k) = PS ρ|δv| (k) 4πk 2 dk.
2
2(5.25)

Here, the convention is that k (not bold-ed) is the magnitude of the vector k. We will plot an“compensated” version of this power spectrum (Lemaster and Stone, 2009), in which the power ismultiplied by k p , where p is the expected power law index of the spectrum (p = 5/3 in our case).
When plotted in this way the spectra will appear flat if following the power law and any deviations
from this power law become more obvious. After compensation, we normalize the spectra such thatthe average is 1.
The power spectra for simulations without particles can be seen in figure 5.3. We see that the
spectra for 3 different driving levels are roughly the same, with slight deviations at large scales. The
compensated spectra show the typical signatures of turbulence: flatter than the expected powerlaw at low k (appearing upward sloped here), following the expected power law at intermediate k,
and steeper than the expected power law at high k as we approach the dissipation scale. In general,the mapping of the spectrum of the artificially driven velocity perturbations to the actual observed
perturbations (which are simultaneously undergoing the turbulent cascade) is a very interesting
problem. While this may merit further study and characterization in the future, we are satisfied
that these spectra are Kolmogorov-like enough for our purposes in this study.

132

Figure 5.3 A demonstration that our stirring algorithm produces the roughly same spectrum for 3
different levels of driving, and that the spectrum roughly follows a p = −5/3 power law over theinertial range. The power plotted here is normalized such that the mean is zero. The power is also
compensated with a k ∝ −5/3 dependence, such that a −5/3 power law would appear perfectly
horizontal on this plot. Time averages are taken from 10Ω−1 to 50Ω−1 .

133
5.2.7

Varied Parameters and Table of Simulations Performed

We perform a series of simulations with two primary goals:• Determine the ability of the streaming instability to form gravitationally-bound particle
clumps in the presence of Kolmogorov (p = 5/3) turbulence of varying strength.
• In cases which form gravitationally-bound clumps, determine the resulting mass-spectrumof the clumps that form and compare this spectrum to the case with no driven turbulence.
In order to achieve these goals, we set up a suite of simulations that vary the strength of the driven
turbulence (αdrive ). Table 5.1 shows the simulations performed and their parameters. We starteach simulation with no particle self-gravity. During this first phase, the streaming instability willconcentrate particles into clumps (to the extent allowed by the turbulence). We find that 20Ω−1 isenough time to let the SI saturate, noting that the maximum particle density has reached a roughly
steady state by this point. After 20Ω−1 , self-gravity is turned on and any clumps that are denseenough will collapse due their own self-gravity. We will then measure the characteristics of theresulting clumps using the methods explained in the next section.

5.3

Analysis Techniques for Clumps and Particles
In this section, we will describe the tools that we will use to measure the characteristics ofclumps in our simulations. This includes identifying gravitationally-bound clumps and carefullymeasure the statistics of the gravitationally-bound clumps.

5.3.1

Finding Gravitationally-Bound Particle Clumps

As described in the introduction, modeling planet formation is a many-stage process, startingwith micron-sized grains and ending with Earth-sized objects. We study an intermediate regime
here, hoping to end up with a shearing box full of planetesimal-sized objects. The most usefulway to characterize the properties of these objects to pass along to the next stage of modeling is
to measure the mass-spectrum of the planetesimals that form. This spectrum can also provide a

134

Label
hr control
hr strong
hr moderate
hr weakhr strong np
hr moderate np
hr weak np
lr control
lr stronglr moderate
lr weaklr strong np
lr moderate np
lr weak np
Table 5.1 Simulation Parameters.

αdrive0.0
10−3
10−3.510−4
10−3
10−3.5
10−4
0.0
10−3
10−3.510−4
10−3
10−3.510−4resolution
5123
5123
5123
5123
5123
5123
5123
1283
1283
1283
1283
1283
1283
1283

Number of Particles
5123
5123
5123
5123
0
0
0
1283
1283
1283
1283
0
0
0

tmax
29
40
38
30
10
10
10
50
50
50
50
50
50
50

135
useful comparison to the mass spectrum of objects in the asteroid and kuiper belts, which havebeen generally found to have a steeper spectrum than is explainable the by the SI (Morbidelli et al.,2009; Fraser et al., 2010, 2014).
To identify gravitationally bound clumps, we use PLanetesimal ANalyzer (“PLAN” for short),a code developed by Rixin Li (Li, in prep.). We will briefly summarize the algorithm here. PLAN
uses the same approach as HOP (Eisenstein and Hut, 1998), a code used to identify dark matter
halos in simulations of galaxies and galaxy clusters. First, particles in dense regions are identified
using a Barnes-Hut tree algorithm (Barnes and Hut, 1986) similar to those used in n-body simulations, which is capable of grouping together physically related particles. For particles whose
immediate surroundings have a particle density greater than a certain threshold, the algorithm
looks for each particle’s densest neighbor. This process is repeated, creating a chain of particles
with each in a denser region than the last, until a true peak is reached. All particles that get
mapped to the same peak are now tentatively considered a single clump. Clump’s whose mass issmaller than the mass required to achieve a Hill radius of at least one grid cell are discarded. The
mass of each clump is then measured by drawing larger and larger spheres centered on the peak.
For each sphere, the interior particle density is calculated. If this average density exceeds the Roche
density, given by
ρr
9
= ,
e
ρ0
G

(5.26)

the clump will be bound by it’s own self-gravity Johansen et al. (2012). The mass contained in thendefined to be the particle mass contained within the largest possible sphere such that the internal
particle density exceeds the Roche density. We note that the Roche density in our simulations
is about 90 ρ0 . If turbulence prevents clumps from achieving this density, we expect that thisalgorithm will find no clumps at all.

5.3.2

Statistics of the Planetesimal Mass Spectrum

Even with a rigorous algorithm to identify clumps, there is no obvious or universal way torepresent the statistics of these clumps. The mass-spectrum of planetesimals has been represented

136

Figure 5.4 Mid-plane slices of the particle density at the onset of the streaming instability in asimulation with no driven turbulence (hr control).

137
in many different ways by different authors, usually taking the form of some variation on a powerlaw. Our plan will be to try several fits and select the best model using an information criteria
analysis, following the general method of Li, in prep.
To start, we define the differential mass-histogram of the planetesimals following the convention of Simon et al. (2016):
dN
2=.
dMp i Mp,i+1 − Mp,i−1(5.27)

where Mp is the planetesimal mass, indexed by i. We also define the cumulative histogram N (> M )to be the number of planetesimals with mass greater than M.
First, we will consider a simple power-law fit. While this is usually not a particularly good
fit to the spectrum, it has been used by authors in the past and will allow for direct comparison
to these studies (Abod et al., 2018; Simon et al., 2016; Schäfer et al., 2017; Simon et al., 2017b).
It is also very simple to calculate and gives a useful first-order look at the data. In this case, thedifferential mass distribution represented as:
dN= C1 Mp−p
dMp(5.28)

where C1 is a normalization factor and p is the power law index of interest. Note that earlier inthe paper we used “p”, as the power-law index for turbulence, but from this point forward it willrefer to this power law index unless otherwise noted. We will use a maximum likelihood estimator
derived by (Clauset et al., 2009) to determine the power law index and associated uncertainty:
X 

Mp,i
p−1
p=1+n
ln
, σp = √ .
(5.29)
Mp,minni

As seen in figure 5.5, the mass-spectra tend to have a downward curvature as M increases,
motivating power law variations that are tapered, truncated, or piece-wise. In this study we willconsider 6 fits to the spectrum (including the simple power law). These functions are usually
expressed as cumulative probability P (> M ), which can be trivially converted into a cumulative
histogram (N (> M )) by multiplying by the number of planetesimals. Here, we will give thefunctional form and free parameters for all the fits, but will refer the reader to the original works
for the physical motivations and subtleties of each function (also see Li, in prep. for a review).

138

Figure 5.5 Typical differential and cumulative mass-histograms of the planetesimals. Data shownis from hr control at t = 27Ω. The black dashed lines show a simply power law fit to the data.

139
The Simple Power Law (SPL) has only one parameter: α. Note that this is the same as the simple
power law fit to the differential histogram with α = p − 1.
P (> M ) =

M
Mmin

−α
(5.30)The Simply Tapered Power Law (STPL, Abod et al. (2018)) has 2 parameters: α and Mexp .
P (> M ) =

M
Mmin

−α



M − Mmin
exp −
Mexp


(5.31)The Variably Tapered Power Law (VTPL, Schäfer et al. (2017)) has 3 parameters: α, β, and Mexp .
P (> M ) =

M
Mmin

−α


exp −

β
M β − Mmin



β
Mexp

(5.32)The Broken Cumulative Power Law (BCPL, Li, in prep.) has 3 parameters: α1 , α2 , and Mbr .

−α1


M

 MminM < MbrP (> M ) =(5.33)


−α2

M

M >M
α −α2

Mbr1

−α

Mmin1

br

The Truncated Power Law (TPL, Schäfer et al. (2017)) has 2 parameters: α and Mtr .
p(M ) =αM(M/Mmin )−α1 − (Mtr /Mmin )−αThe Broken Power Law (BPL, Li, in prep.) has 3 parameters: α1 , α2 , and Mbr .
 
−α1


C0
M






 M MminM < Mbr1
1
1
Mbr −α1 −1p(M ) =, C0 =
+
−

α1
α2 α1
Mmin
 C0 (M/Mmin )−α2

M
M
>
M
α
−α
br
(M /M
) 1 2
br

(5.34)(5.35)min

We fit all of these models using a maximum likelihood method based on the likelihood functions given in Li, in prep. These functions allow us to determine the likelihood of a model withsome given set of parameters denoted as θ. The parameters of maximum likelihood are then be determined using a nested-grid-search algorithm. To summarize, this algorithm starts by calculating
the likelihood of the model for many parameter sets θ, arranged on large, widely-spaced grid. This
first iteration finds the general vicinity of the maximum likelihood model in this parameter space.
The search is then refined to a smaller grid with tighter spacing, centered around the maximum

140
likelihood parameters from the first search. This process is repeated 5 times, yielding a maximum
likelihood parameter vector θ in which each parameter is optimized with an accuracy better thana part in 1000.
We use a non-parametric bootstrapping method to determine the uncertainties in the parameters (Efron and Tibshirani, 1994a; Burnham and Anderson, 2002a). This involves re-sampling thesample of clump masses from a given simulation 1000 times and then using the nested-grid-search
to determine the maximum likelihood parameters for each re-sampling. This generates a 1000 point
distribution for each parameter. The 1σ confidence interval is then is then defined to be the region
between the 16th percentile and 84th percentile of this distribution for each parameter.
We evaluate the goodness of fit of the various models using Bayesian and Akaike information
criteria (Kass and Raftery, 1995; Akaike, 1974). These criteria account for the number of fit
parameters in each model, giving an advantage to simpler models. They are defined as follows:
BIC = K ln N − 2 ln L,(5.36)AIC = 2K − 2 ln L,(5.37)

where K is the number of parameters in the model, N is the number of data points (clump masses inthis case), and L is the likelihood of the maximum likelihood model. These values do not bear muchmeaning in isolation, as the likelihood values are affected by arbitrary constant and the sample size.
To define a metric that has physical meaning, we measure the differential between the model withthe minimum BIC/AIC and all other models, defined as∆BIC = BIC − BICmin ,(5.38)∆AIC = AIC − AICmin .

(5.39)

With this definition, the best model will have ∆BIC and/or ∆AIC = 0. In most cases we findthat both criteria identify the same best model. The other models will have some positive value
for these criteria, with higher values corresponding to worse models. While these criteria aresomewhat abstract, as a general rule we consider there to be at least some evidence for models with

141
∆BIC (∆AIC ) between 0 and 10, and essentially no evidence for models with values greater than 10.
For a more nuanced interpretation of these criteria, we refer the reader to Burnham and Anderson(2002a).

5.3.3

When to Measure the Clumps

A issue that is not discussed enough in the literature is the somewhat arbitrary choice ofat what point in time to measure the statistics of the planetesimals. In general, two contrasting
factors must be weighed:(1) The behaviour that our simulations should accurately capture is the clumping of the streaming instability and subsequent collapse or non-collapse of each of those clumps in isolation.
This is the mass-spectrum that we want to measure. However, not all SI seeded clumps
collapse instantly when self-gravity is turned on. Even with no driven turbulence, the firstclumps become gravitationally bound after a dynamical time or two, and the number ofclumps continues to increase for roughly 10 dynamical times. Giving the clumps more time
to form increases the number of clumps present and makes our statistics more robust.
(2) After some time, the number of clumps will saturate and stay roughly steady while the total
mass of all clumps will continue to increase as planetesimals accrete and merge. A majordifficulty here is that this phase of evolution is not entirely physical. We are restricted by
the grid-scale: in general, clumps will not collapse further than approximately the gridscale because there is no sub-grid self gravity. Even at a respectable resolution of 5123 in a
box that is 0.2h on a side, the grid-scale is ∼ 4 ∗ 10−4 h. At a radius of 10AU in a disk withan aspect ratio of h/r = 0.05, this corresponds to a distance of ∼ 3 ∗ 104 km. This scale isabout 5 Earth radii, and is on the order of 1000 times larger than a typical planetesimal(∼ 1 to 100 km). Due to this non-physical scale, we do not expect our simulations tocorrectly simulate the evolution of planetesimals after their formation.
The selection of the time to measure the clumps’ statistics involves weighing these two factors:

142
we must allow enough clumps to form to do good statistics, but we really don’t want to capture
the post-formation behaviour of the clumps. As a result, we will try to select a time around or
slightly before the number of clumps saturates, but will keep in mind that these statistics are timedependent and will comment on their time-dependence throughout the analysis. Additionally, in§5.5 we propose a better way of measuring the initial masses.

5.4

Results I: Snapshots
In this section we will perform a snapshot-based analysis of the simulations, investigating theeffect that the turbulence has on the formation of clumps. Note that throughout the rest of thispaper, all times will be defined such that the time at which self-gravity is turned on is defined to
be 0.

5.4.1

The Turbulence Threshold for Planetesimal Formation

After running our high resolution simulations through the PLAN algorithm, we find thatsome amount of planetesimals (or “clumps”, these words are being used interchangeably in thiswork) form in every simulation except one. No gravitationally-bound clumps in the simulation
being driven at α = 10−3 after 20 Ω−1 with self-gravity on. While this is an incredibly simple
and qualitative result, it is probably the most important finding of this work. The simulation
being driven at α = 10−3.5 does form gravitationally bound clumps, meaning that the thresholdfor turbulence to disrupt planetesimal formation via the streaming-instability/self-gravity (SI/SG)
mechanism is somewhere between α = 10−3.5 and 10−3 . Equivalently, the threshold in terms ofδv/cs would be between 0.018 and 0.032.
This threshold is at a fairly modest value of α, making the implications very interesting. In
general, we expect a fully MRI active disk (ionized throughout) to easily exceed this threshold.
However, for most radii within a disk, we expect non-ideal MHD effects to damp the MRI at the
mid-plane, creating a dead zone. The mid-plane turbulence in dead zones will depend on several
factors including the column density of the disk and ionization state of the gas, but for most disks

143
we would expect a level of turbulence less than this threshold, permitting planetesimal formation.
We save a full discussion of the consequences and implications of this threshold for §5.6.

5.4.2

Mass Fraction in Planetesimals

To quantify the extent to which planetesimals form, a simple and useful metric is the total
amount of mass that ends up in planetesimals relative to the total mass of all particles in the simulations: Mplan /Mpar . We find that this ratio is significantly effected by the presence of turbulence.
Table 5.2 shows a summary of the results from our simulations, including the approximate massfraction. We should note that this metric suffers from the time-dependence problem explained in§5.3.3, making it difficult to decide when to measure it. For this section, we select the snapshots tomeasure by-eye, with these snapshot times being listed in table 5.2. Figures 5.6 and 5.7 show thetime-evolution of the planetesimal mass fraction and demonstrate the time at which we choose tomeasure.
In our control run, about 10% of the mass is contained in planetesimals at the fiducial time,compared to 6% and 2% at αdrive = 10−4 and 10−3.5 respectively. These smaller mass ratios are
despite waiting longer to make the measurement in the runs with turbulence. This result begs thequestion: does turbulence cause smaller clumps to form, or less clumps? The number of clumps
measured is also shown in table 5.2, and seems to scale fairly directly with the mass fraction. This
indicates that all of the simulations tend to form the same size clumps on average, with the majordifference being the number of clumps formed.

5.4.3

Differential Mass Spectra and Simple Power Law Fits

We will now take a first look at the mass spectra by fitting a simple power law to thedifferential mass histogram using the maximum likelihood estimator defined in §5.3. These fits tothe control and fiducial runs are shown in figures 5.6 and 5.7 respectively, along with the timeevolution of the power law index.
As demonstrated by the time-evolution of the power law indices, this is a difficult measure-

144

Figure 5.6 An analysis of our control run: hr control, which does not have any driven turbulence.
The top 3 plots from top to bottom show: the time evolution of the slope of the simple power law
fit to the differential histogram (with a 1σ error bar), the mass fraction of planetesimals, and thenumber of gravitationally bound clumps. The horizontal dashed lines roughly show the range ofpower law indices determined by previous studies. The dashed blue line shows the time at which
we choose to measure the spectrum, with that spectrum shown at the bottom in blue. Similarly,
the red line shows where the the spectrum shown in red is taken. The oval denotes the intermittent
population.

145

Figure 5.7 An analogous figure to the previous one for our fiducial run: hr moderate.

Label
hr control
hr weakhr moderate
hr strong

αdrive0.0
10−4
10−3.510−3

tfirst
1.6
2.0
5.0
never

tsnap
5.3
9.1
16.7
n/a

Nclumps
92
54
20
n/a

Mclumps /Mp
10%
6%
2%
0%

SPL slope1.462 ± 0.048
1.476 ± 0.065
1.586 ± 0.131
n/aTable 5.2 High level summary of clumping behaviour in high resolution simulations. tfirst denotes
the time of the formation of the first clump. tsnap shows the time at which we measure the values
for the next 3 columns.

146
ment to make, especially picking a frame to measure by-eye. Initially, while there are only a handfulof clumps in the box, the slope is quite steep. The slope comes down rather quickly and approaches
a somewhat constant value, although even this is quite sloppy. In figure 5.6, the slope seems toalmost oscillate between about 1.5 and 1.3. This behaviour is explained by the existence of a lowmass population of planetesimals, circled in the red histogram but absent from the blue one. This
population is intermittent, with clumps in this regime not being consistently identified by PLAN.
As a result, it is difficult to trust that these are bonafide clumps that should be included in theanalysis, so we try to select a frame in which this population is minimized. This is still somewhat
subjective though, and is part of the reason we will propose a different approach to the analysis inthe §5.5.
In this snapshot-based analysis, we find that the presence of turbulence does not change thepower law index of the spectrum of planetesimals that form due to the SI/SG mechanism. The
power law indices, shown in table 5.2, are all consistent with each other to 1σ. Furthermore, these
indices are also consistent with previous studies of streaming instability clumping that do not have
driven turbulence (Schäfer et al., 2017; Simon et al., 2016; Abod et al., 2018; Simon et al., 2017a),
who find p values in a range from 1.4-1.7.

5.4.4

More Sophisticated Fits to the Cumulative-Mass-Spectrum

We will now try a several different functional forms, defined in §5.3, to fit the cumulativemass-spectrum. The spectra and fits are shown in figure 5.8, with data again being taken from thesnapshot times in table 5.2. We use the information criteria analysis to determine the best models
and report these models, along with their parameters, in table 5.3.
The spectra in the simulations with no turbulence and weak turbulence favor the broken
power law model based on the information criteria, while the moderate turbulence simulation
favors more “straight” power laws like the TPL and SPL. However, if we compare like-models
across the 3 simulations (for example the STPL and BCPL), we see that the distributions are
generally consistent with each other. This is, in part, due to the relatively large errors on the147parameters resulting from the low sample size, in particular in the case of moderate turbulence.
So while the favored models may be different, this consistency in parameters leads us to believethat these distributions are not physically different from one another. Furthermore, these results
are consistent with previous studies that did not include forced turbulence. Abod et al. (2018) fits
the STPL to their distribution, and for our choice of the Π parameter finds that α = 0.29 ± 0.05
and Mexp = 0.35 ± 0.05, both of which are consistent with all of our simulations. These results inthe presence of turbulence provide even more evidence for a sense of “universality” (Simon et al.,2017a) to the mass-spectrum. While the number of clumps that form may change, the shape of thespectrum does not.

5.5

Results II: Clump Tracking
As pointed out in the previous 2 sections, there are problems with drawing conclusions looking

at only single snapshots (or “frames”) of these simulations. In this section, we propose a newparadigm for measuring the mass spectrum of the clumps that should, in principle, give the mostreliable results for the mass spectrum. We develop an algorithm to track clumps from one frame to
the next and create a mass-history for each clump. From this data, we can then look at every clumpthat ever existed in the simulation and measure it’s mass *at the time of formation*. Thus, we
can truly measure an initial-mass-spectrum that does not include any non-physical post-formation
accretion or merger physics.

5.5.1

Methods and Algorithm

Tracking clumps in this way is a non-trivial problem. PLAN can output lists of the particles
belonging to each clump it finds, but does not use a consistent clump labeling scheme from frame
to frame. Thus, at each frame we are provided a list of clumps, and a list of particles belonging
to each clump, but no innate way to say that a clump in one frame is the same clump seen in aprevious frame. We develop an algorithm to look through these particle lists and track clumps overtime, outlined as follows:

148

Label

hr control

hr weakhr moderate

Best Fit By...
BPL

Parameters
α1 = −0.48+0.29−0.20

α2 = 1.25+0.52−0.14

Mbr = 0.14+0.10−0.04(BCPL)

α1 = 0.283+0.083−0.063α2 = 1.18+0.38−0.19

Mbr = 0.100+0.027−0.052(STPL)

α = 0.15+0.12−0.06

Mexp = 0.261+0.087−0.095BPL

α1 = −0.80+0.38−0.25

α2 = 1.48+1.1−0.21

STPL (∆IC = 9.2)

α = 0.09+0.19−0.13

Mexp = 0.23+0.20−0.12(BCPL)

α1 = 0.21+0.15−0.03

α2 = 1.19+0.94
−0.05

TPL

α = 0.26+0.15−0.23

Mtr = 0.98+0.00−0.43

STPL (∆IC = 2.5)

α = 0.35+0.14−0.18

Mexp = 0.37+0.19−0.23

SPL (∆IC = 2.9)

α = 0.59+0.13−0.06

BCPL (∆IC = 6.4)

α1 = 0.28+0.15−0.03

α2 = 0.91+0.38−0.34

Mbr = 0.15+0.14−0.06
Mbr = 0.082+0.050−0.028

Mbr = 0.045+0.069−0.021

Table 5.3 Summary of fits to the cumulative mass spectrum of planetesimals based on single snapshots from our high resolution simulations. ∆IC denotes the average of the two information criteria
considered. We give the parameters for the STPL and BCPL for all runs regardless of ∆IC as away of comparing from run-to-run and because the interpretation of their parameters is intuitive.

149

Figure 5.8 Cumulative spectra and models for snapshots of hr control, hr weak, and hr moderatefrom top to bottom. Fit parameters to the best models can be found in 5.3. The square points
show the characteristic mass scale relevant to each fit. The bottom histogram shows the true values,
with each plot above offset by a factor of 10 for convenient visualization.

150(1) At the starting frame, read in the particle lists corresponding to each existing clump andgive each clump an internal ID that will be consistent from frame to frame.
(2) Read in the particle lists for each clump in the following frame. For each clump in thecurrent frame, count the number of overlapping particles of each clump in the new frame.
This presents several possible outcomes:• Case A: Only one new-frame clump has overlapping particles with the prior clump.
We identify these as the same clump and update the clump’s particle list with the
new-frame particles.
• Case B: Multiple new-frame clumps have overlapping particles with the prior clump,
indicating there has been a split of some sort. In this case, we identify the clump
with more overlapping particles as the same clump and update the clump’s particle
list with the new-frame particles. If the other clump with overlapping particles is
newly-formed this frame and more than 50% of it’s particles came from the parent
clump, it is flagged as a “splitter”. We will track this clump from this frame on, but
it will not be counted towards the statistical analysis.
• Case C: No new-frame clumps have overlapping particles with the prior clump. The
clump has “disappeared”. It could have been destroyed by turbulence or an interaction
with other clumps, or it could have been very marginally gravitationally bound in thefirst place such that PLAN would identify it as gravitationally bound in one frame,
but not the next. This happens relatively frequently. In this case, we do not match
this clump with any new-frame clumps and do not update it’s particle list. However,
we will continue to look for this clump in subsequent frames with the old particle list.
(3) At this point we have mapped clumps that existed in the previous frame to the new frame.
Now, we need to deal with new clumps (ie. clumps that are “unclaimed” by clumps thatexisted in the previous frame). We read in the particle lists for these new clumps and will
look for their matches in subsequent frames.

151(4) Continue to step through the frames by repeating steps 2 and 3.
This scheme is illustrated graphically in figure 5.9. We would like to emphasize that measuringinitial masses with this method is able to solve 3 problems stemming from the artificially large sizeof clumps in these simulations:(1) Accretion: After formation, clumps will accrete particles at a super-physical rate. Being
able to look only at the initial mass ignores this accretion.
(2) Mergers: Clumps will merge at an enhanced rate. If two clumps merge, a snapshot would
count them as one large clump. In our tracking method, we measure the initial formation
mass of the two separate clumps.
(3) Splits: If a clump was able to fully collapse, one would expect the rate of clumps splitting or
fragmenting to be quite low. However, we see it quite often with our artificially large clumps.
Our algorithm is able to ignore clumps that formed as the result of a split, measuring only
the initial mass of the parent clump.

5.5.2

Results from Clump Tracking

In this section we repeat a a very similar analysis to section §5.4, only this time using only
the initial masses of the planetesimals and ignoring clumps that formed from a split. With thisprocedure we not longer have to somewhat arbitrarily select a snapshot by-eye, making the analysis
more robust. Figure 5.10 shows the cumulative spectra along with the series of models we fit to it.
The model parameters are given in table 5.4.
First we will just consider the effect of this analysis on the distribution in the control simulation. Comparing to the snapshot analysis, we see that the masses of the clumps have generally
decreased by a factor of 3 or 4, with the average clump mass changing from 0.55 to 0.18MG . The
characteristic masses found by the models change by similar factors as well. This is quite significant, and implies that previous studies have over-predicted the masses of planetesimals formed by

152

Figure 5.9 An illustration of our particle tracking scheme. The example shown here is intended to
be demonstrative of our scheme and does not represent actual data from a simulation. Each circle
represents a gravitationally-bound clump. Red corresponds to a new clump, gray to a previously
found clump that was identified again, white to a clump that was not found (but will continue tobe looked for in future frames), and blue to a new clump that formed from a split. Clumps are
consistently labeled from frame to frame. Clumps that merge are given 2 labels and are continued
to be tracked as both clumps. In taking statistics, we use the masses of all of the red clumps.

153
the SI/SG mechanism due to super-physical accretion and merging. This is well quantified by α1of the BCPL model shown in the tables, which changes from 0.283 in the snapshot for the control
run to 0.065 with then new analysis. In fact, the low-mass slopes of the BCPL are all within a few
σ of zero. This effect is mostly due to the removal of clumps that form as the result of splits, which
are a non-physical process in these simulations and lead to an over-abundance of planetesimals at
the low-mass end of the spectrum.
We will now consider the original point of this study using this new analysis: does turbulencechange the mass-spectrum of planetesimals? Inspecting the parameters in table 5.4, we find thatthe clump tracking results generally confirm the null result found in the previous section. In thiscase, the clumps in each simulation favor a broken power law description, and the parameters of
the BPL are consistent across the 3 simulations with the exception of Mbr , which shows a fairly
unremarkable deviation for the moderate turbulence case in comparison to the control.
While the nature of this new analysis technique makes meaningful comparisons to previous
studies somewhat moot (as they did a different analysis), we would like to emphasize that it
is possible that this measurement of strictly the initial masses might change the conclusions ofprior work on the SI/SG mechanism. While prior studies have generally confirmed that the mass
spectrum is invariant with respect to most parameters that control the behaviour of the SI, it ispossible that their measurement of the spectrum is dominated by post-formation evolution of theclumps. This post formation evolution is not necessarily dependent on the SI parameters in thesame way that the initial masses are, and as a result these studies may have been missed relevant
effects that do influence the true initial mass spectrum.

5.6

Conclusions and Discussion
Our primary results from this chapter are as follows:(1) Turbulence at the level of α = 10−3 does not permit the formation of gravitationally boundparticle clumps via the SI/SG mechanism. Due to the velocity dispersion, the SI is unable

154

Figure 5.10 An analogous figure to figure 5.8, again for hr control, hr weak, and hr moderate from
top to bottom. In this figure, only the initial mass of each clump is used, as determined by ourclump tracking algorithm. The square points show the characteristic mass scale relevant to each
fit.

155

Label

hr control

hr weakhr moderate

Best Fit By...
BPL

Parameters
α1 = −1.54+0.17−0.23

α2 = 1.612+0.18−0.14

Mbr = 0.04204+0.0042−0.0044(BCPL)

α1 = 0.065+0.021−0.017α2 = 1.06+0.14−0.07

Mbr = 0.0189+0.0024−0.0025(STPL)

α = −0.0474+0.0080−0.0094

Mexp = 0.0463+0.0028−0.0043BPL

α1 = −1.72+0.54−0.38

α2 = 1.35+0.45−0.11

Mbr = 0.040+0.013−0.015(BCPL)

α1 = 0.0834+0.033−0.019α2 = 1.20+0.18−0.13

Mbr = 0.027+0.002−0.012(STPL)

α = −0.012+0.082−0.038Mexp = 0.078+0.033−0.041BPL

α1 = −1.70+0.49−0.47

α2 = 1.63+0.30−0.27

Mbr = 0.0377+0.002−0.027(BCPL)

α1 = 0.09+0.16−0.02

α2 = 1.25+0.46−0.07

Mbr = 0.0226+0.0051−0.0152(STPL)

α = −0.046+0.051−0.034Mexp = 0.047+0.006−0.034

Table 5.4 Summary of fits to the cumulative mass spectrum of planetesimals based on the initial
masses of each clump from our clump-tracking analysis. We give the parameters for the STPL
and BCPL for all runs regardless of ∆IC as a way of comparing from run-to-run and because theinterpretation of their parameters is intuitive. For all cases here, no fits other than the BPL had∆IC < 10 For reference, the SPL fits to the differential spectra for this method are 1.372 ± 0.023,
1.357 ± 0.035, and 1.380 ± 0.047 for the 3 runs shown from top to bottom, respectively.

156
to concentrate particles to a sufficient density for clumps to be bound by their own self
gravity. The threshold of turbulence above which planetesimals do not form is somewhere
between α = 10−3.5 and α = 10−3 , or equivalently δv = 0.018 and 0.032.
(2) Intermediate turbulence (α = 10−4 to α = 10−3.5 ) permits the formation of planetesimals,
but they form more slowly, are fewer in number, and contain less total mass than thosethat form in the control simulation.
(3) The mass spectrum of planetesimals is not greatly affected by the presence of turbulence. However, clumps in moderate turbulence show a preference for ”straight” powerlaws whereas clumps in weaker turbulence or no turbulence are better fit by broken power
laws.
(4) Due to the time-variance of the planetesimal characteristics and the non-physical post
evolution behaviour of clumps in these simulations, analyzing snapshots is not a very good
way of measuring the initial mass spectrum. We propose a method in which clumps canbe tracked, allowing us to measure only the true initial masses and ignore non-physicalprocesses. This method reveals that other studies have likely overestimated masses by afactor of ∼ 3, and that the slope of the mass spectrum at low masses is actually shallower
than previously determined.
As stated previously, the most striking result from this work is the turbulence threshold for
disrupting the SI/SG planetesimal formation mechanism. The α = 10−3 threshold is at a fairly“intermediate” value of α, making the implications very interesting. We will consider a few scenarios
here, assuming that the main driver of turbulence in protoplanetary disks is the MRI, and keeping
in mind that the SI/SG mechanism is occurring at the mid-plane of the disk. First, we considera region of a disk that is ionized throughout it’s vertical extent. This is expected to be the casevery close to the star, and will result in MRI activity throughout the whole column. In this case,the expected level of turbulence at the mid-plane is about α = 10−2 , or possibly greater if a strong

157
net-vertical-field exists (Simon et al., 2013a; Salvesen et al., 2016). Based on our threshold, theSI/SG mechanism cannot form planetesimals in a fully ionized region of the disk. This provides
evidence in favor of migration scenarios (Masset and Papaloizou, 2003; Beaugé and Nesvorný, 2012)
over in situ formation arguments (Batygin et al., 2016) for very close in planets like hot Jupiters.
As discussed in all other chapters of this thesis, non-ideal MHD effects are expected to greatlychange the behaviour of the MRI if the disk is not sufficiently ionized. This introduces the “deadzone” (DZ) paradigm: the disk will have MRI-driven turbulence at it’s radiatively ionized regions
at high |z|, and a DZ at the mid-plane due to non-ideal effects damping or stopping the MRI (seechapters 2 and 4). In general, we expect the DZ to be dominated by the Ohmic term at higher
densities (smaller radii) and the Ambipolar term at lower densities (larger radii). As seen in chapter
2, Ohmic DZs can have interesting fluid motions that are not necessarily well characterized by asingle parameter like α. We find that smaller disks tend to have more active mid-planes due totrickle-down-turbulence from the active regions, but still only reach roughly α ∼ 10−4 , or possibly
up to a few ×10−4 if the large scale-velocity patterns are included. As a result, we expect the
SI/SG mechanism to be able to form planetesimals fairly efficiently in an Ohmic DZ for disks ofmost masses, but may struggle slightly in the smallest disks.
In chapter 4, we considered fully ambipolar disks as well as disks with MRI-active upper layers
and an ambipolar damping region. As a reminder: Am is the parameter we use to characterize
the ionization level, with smaller Am values corresponding to less ionization and more damping.
Inspecting figures 4.4 and 4.5, we can see that, for a mid-plane of a fully ambipolar disk to be below
the δv threshold we have established here, Am must be about 0.1 or less. For a disk with FUV
ionized layers, figure 4.9 shows that δv/cs is about 0.04 or 0.05 at the mid-plane, just barely above
the threshold for planetesimal formation via SI/SG. Am turns out to be an interesting control
parameter for planet formation in this scenario, motivating further chemical and dust modeling(perhaps similar to Zhu et al. (2015)) to determine the appropriate Am value at the mid-plane ofthe outer disk.

Chapter 6

Summary and Looking Forward

This thesis has discussed several important topics relating to protoplanetary disks. Chapters
2-4 each take a close look at phenomena associated with non-ideal MHD effects in disks, while
chapter 5 places important constraints on the conditions of the disk that permit planetesimalformation. Here, we will briefly summarize the conclusions of each chapter of this thesis and thenconsider some of the most exciting prospects for the future of disk and planet formation theory.

6.1

Summary
In chapter 2, we perform a shearing box study investigating Ohmic dead zones and theinteraction between the active region and the dead zone. We find that the Reynolds stress at the
mid-plane of the disk scales strongly with the thickness of the dead zone (ie. the total mass of thedisk), with very thick dead zones having effectively negligible stresses and thin dead zones having
turbulent stresses on the order of α ∼ 10−4 . We also find that, while the fluid in the dead zone has
a considerable amount of kinetic energy, most of this energy is contained in non-turbulent motions
that do not transport angular momentum transport. We identify these large-scale motions as thelargest scale “r-modes” Lubow and Pringle (1993), found by considering the modes in a purelyhydrodynamic disk. While these modes may be amplified in a shearing box compared to a realdisk, we note that they would produce an interesting warping effect on settled particle layers.
In chapter 3, we construct a 1D model that couples a classic Shakura-Sunyaev-like disk modelfor the steady state of the disk with time evolution models for the net magnetic field and surface

159
density of the disk. We use this model to study the interaction between the stellar magnetic cycle
and the disk, and specifically the ability of the cycle to produce EXor-like events as proposed byArmitage (2016). We find that this mechanism is able to produce cyclic outbursts in disks with
sufficiently low accretion rates (∼ 10−8 M /yr or less ). These outbursts produce disk-luminosity
enhancements of roughly an order of magnitude lasting a quarter of the stellar magnetic cycle period
or less. While these results are somewhat promising at a “proof of concept” level, the accretion
rates required are quite low and more modeling is needed to account for the erratic behaviour ofsmall stars.
In chapter 4, we use shearing boxes to study turbulence in the outer regions of PPDs, where
ambipolar diffusion is the dominant non-ideal effect. Our simulations with FUV-ionized surface
layers demonstrate that one should use caution in constructing 1D models of disks that assume eachconstant-z-slice of the disk can be treated as local, obeying already established physics at a given
density and ionization fraction (Bai, 2011b). Vertical transport within the disk plays an importantrole. For mid-planes with Am < 1 the behaviour of the fluid is qualitatively different than in
unstratified simulations or simulations with no FUV layers: Magnetic stresses are dominated bylarge scale fields and the mid-plane has a magnetic field strength in excess of that producible by
the local MRI. This indicates the the MRI is completely arrested and is not even acting in it’sdamped capacity. Finally, we show that toroidal magnetic field from the FUV layers is able to betransported to the mid-plane via diffusive effects for Am < 10, despite magnetic buoyancy trying
to carry magnetic flux away from the mid-plane.
Finally, in chapter 5 we study the effect of turbulence on the formation of planetesimals via
a combination of the streaming instability and self gravity. While in general the mass spectrum ofplanetesimals does not change drastically based on the level of turbulence, we see that the total
amount of mass the goes into planetesimals does change drastically, with no planetesimals forming
at a turbulence level of α = 10−3 . We also point out several flaws in the way that these typesof simulations are typically analyzed, relying heavily on the statistics of single snapshots that are
selected by-eye. We propose and implement a better way of measuring the planetesimals that

160
involves tracking each gravitationally-bound-clump over time and measuring only it’s initial mass.
This generally reduces the masses found by a factor of ∼ 3. In addition, this method yields a flatter
slope in the cumulative mass spectrum at the lowest masses, finding that many of the planetesimals
formed in that regime were the result of a splitting event, which will occur at a super-physical ratein these simulations.

6.2

Looking Forward
This is an exciting time in protoplanetary disk and planet formation theory. For many

years, theorists in this field worked without ever seeing a spatially resolved image of a disk in any
capacity, and with only one example of a planetary system. Both of these things has changed in
the past 10 years. Thanks to advances in interferometry, spatially resolved images of disks now
available (ALMA Partnership et al., 2015), along with incredible measurements of disk properties
like turbulent velocity (Flaherty et al., 2015, 2017, 2018) and the distribution of various chemical
abundances throughout the disk (Schwarz et al., 2016; Huang et al., 2018; Cleeves et al., 2015).
In addition, thanks to monumental efforts by many groups all around the world, we have now
found and characterized hundreds of exoplanet systems, a huge increase in sample size compared to“one”, which was the case for most of the history of planet formation theory (for a review, see Winnand Fabrycky (2015)). These advances have sent theorists back to the drawing board on several
topics including determining the primary mechanism of angular momentum transport in disks (theweak turbulence problem) and understanding key aspects of planet formation and planetary systemevolution.
Thinking in this very broad context, I hope that this thesis has made meaningful contributions
to these areas in which much theory work is yet to be done. I would hope that chapters 2 and
4 of this work has contributed to the addressing the weak turbulence problem by pinning downsome of the finer points relating to non-ideal MHD effects in disks and their effect on the MRI andturbulence. In chapter 5, our work has for the first time studied planetesimal formation via thestreaming insatiability in a more realistic context with turbulence, placing important constraints

161
on the conditions of planetesimal formation. This result also re-emphasizes the importance ofunderstanding non-ideal MHD effects and the level of turbulence within dead zones.
In the world of disk simulations, it has long been impossible to perform simulations that areboth global and high enough resolution to resolve relevant small-scale physics like the MRI. Thesetypes of simulations are now entering the realm of possibility (eg. Zhu and Stone (2018); Hu et al.
(2019)), thanks to advances in both algorithms (Felker and Stone, 2018) and raw computational
power. Over the next several years, non-ideal global MHD simulations are likely to be incredibly
helpful in linking local and global disk physics and will hopefully give further insight into therelative roles of winds and MRI turbulence in transporting angular momentum.
While there are many aspects of planet formation that remain active areas of research, theplanetesimal formation process remains one of the most interesting and puzzling problems. We
are near a point of saturation for studies similar to chapter 5, due to the inherent limitations ofsimulations on a uniform grid and the implementation of the self-gravity solver. The next generationof simulations must have either sub-grid self gravity, adaptive mesh refinement (Felker and Stone,2018), or some other method that allows the full gravitational collapse of planetesimals. This is
not an easy problem, as self-gravity is inherently non-local, computationally expensive, and hardto parallelize efficiently. However, simulations of this type have become necessary to make further
progress on the topic of planetesimal formation.

Bibliography

C. P. Abod, J. B. Simon, R. Li, P. J. Armitage, A. N. Youdin, and K. A. Kretke. The Mass and Size
Distribution of Planetesimals Formed by the Streaming Instability. II. The Effect of the Radial
Gas Pressure Gradient. arXiv e-prints, October 2018.
H. Akaike. A New Look at the Statistical Model Identification. IEEE Transactions on AutomaticControl, 19:716–723, 1974.
V. Alexiades, G. Amiez, and Pierre-Alain Gremaud. Super-Time-Stepping acceleration of explicit
schemes for parabolic problems. Communications in Numerical Methods in Engineering, 12:
31–42, 1996.
ALMA Partnership, C. L. Brogan, L. M. Pérez, T. R. Hunter, W. R. F. Dent, A. S. Hales, R. E.
Hills, S. Corder, E. B. Fomalont, C. Vlahakis, Y. Asaki, D. Barkats, A. Hirota, J. A. Hodge,
C. M. V. Impellizzeri, R. Kneissl, E. Liuzzo, R. Lucas, N. Marcelino, S. Matsushita, K. Nakanishi, N. Phillips, A. M. S. Richards, I. Toledo, R. Aladro, D. Broguiere, J. R. Cortes, P. C.
Cortes, D. Espada, F. Galarza, D. Garcia-Appadoo, L. Guzman-Ramirez, E. M. Humphreys,
T. Jung, S. Kameno, R. A. Laing, S. Leon, G. Marconi, A. Mignano, B. Nikolic, L.-A. Nyman,
M. Radiszcz, A. Remijan, J. A. Rodón, T. Sawada, S. Takahashi, R. P. J. Tilanus, B. Vila Vilaro,
L. C. Watson, T. Wiklind, E. Akiyama, E. Chapillon, I. de Gregorio-Monsalvo, J. Di Francesco,
F. Gueth, A. Kawamura, C.-F. Lee, Q. Nguyen Luong, J. Mangum, V. Pietu, P. Sanhueza,
K. Saigo, S. Takakuwa, C. Ubach, T. van Kempen, A. Wootten, A. Castro-Carrizo, H. Francke,
J. Gallardo, J. Garcia, S. Gonzalez, T. Hill, T. Kaminski, Y. Kurono, H.-Y. Liu, C. Lopez,
F. Morales, K. Plarre, G. Schieven, L. Testi, L. Videla, E. Villard, P. Andreani, J. E. Hibbard,
and K. Tatematsu. The 2014 ALMA Long Baseline Campaign: First Results from High Angular Resolution Observations toward the HL Tau Region. Astrophysical Journal, Letters, 808:L3,
July 2015. doi: 10.1088/2041-8205/808/1/L3.
J.

J. Aly.
Electrodynamics of disk accretion
Astronomy and Astrophysics, 86:192–197, Jun 1980.

onto

magnetic

neutronstar.

P. J. Armitage. Lecture notes on the formation and early evolution of planetary systems. arXiv
Astrophysics e-prints, January 2007.
P. J. Armitage. Dynamics of Protoplanetary Disks. Annual Review of Astron and Astrophys, 49:
195–236, September 2011. doi: 10.1146/annurev-astro-081710-102521.
P. J. Armitage.
EXor Outbursts from Disk Amplification of Stellar Magnetic Cycles.
Astrophysical Journal, Letters, 833:L15, December 2016. doi: 10.3847/2041-8213/833/2/L15.

163
P. J. Armitage, M. Livio, and J. E. Pringle. Episodic accretion in magnetically layered protoplanetary discs. Monthly Notices of the RAS, 324:705–711, June 2001. doi: 10.1046/j.13658711.2001.04356.x.
Philip J. Armitage. Protoplanetary disk evolution, page 65 to 108. Cambridge University Press,2009. doi: 10.1017/CBO9780511802225.004.
M. Audard, P. Ábrahám, M. M. Dunham, J. D. Green, N. Grosso, K. Hamaguchi, J. H. Kastner, Á. Kóspál, G. Lodato, M. M. Romanova, S. L. Skinner, E. I. Vorobyov, and Z. Zhu.
Episodic Accretion in Young Stars. Protostars and Planets VI, pages 387–410, 2014. doi:
10.2458/azu uapress 9780816531240-ch017.
X.-N. Bai.
Magnetorotational-instability-driven Accretion in Protoplanetary
Astrophysical Journal, 739:50, September 2011a. doi: 10.1088/0004-637X/739/1/50.

Disks.

X.-N. Bai. The Role of Tiny Grains on the Accretion Process in Protoplanetary Disks.
Astrophysical Journal, 739:51, September 2011b. doi: 10.1088/0004-637X/739/1/51.
X.-N. Bai. Wind-driven Accretion in Protoplanetary Disks. II. Radial Dependence and Global
Picture. Astrophysical Journal, 772:96, August 2013. doi: 10.1088/0004-637X/772/2/96.
X.-N. Bai. Hall-effect-Controlled Gas Dynamics in Protoplanetary Disks. I. Wind Solutions at theInner Disk. Astrophysical Journal, 791:137, August 2014. doi: 10.1088/0004-637X/791/2/137.
X.-N. Bai. Hall Effect Controlled Gas Dynamics in Protoplanetary Disks. II. Full 3D Simulations toward the Outer Disk. Astrophysical Journal, 798:84, January 2015. doi: 10.1088/0004637X/798/2/84.
X.-N. Bai and J. Goodman.
Heat and Dust in Active Layers of Protostellar Disks.
Astrophysical Journal, 701:737–755, August 2009. doi: 10.1088/0004-637X/701/1/737.
X.-N. Bai and J. M. Stone. Particle-gas Dynamics with Athena: Method and Convergence. Astrophysical Journal, Supplement, 190:297–310, October 2010a. doi: 10.1088/00670049/190/2/297.
X.-N. Bai and J. M. Stone. Dynamics of Solids in the Midplane of Protoplanetary Disks: Implications for Planetesimal Formation. Astrophysical Journal, 722:1437–1459, October 2010b. doi:
10.1088/0004-637X/722/2/1437.
X.-N. Bai and J. M. Stone. Effect of Ambipolar Diffusion on the Nonlinear Evolution of Magnetorotational Instability in Weakly Ionized Disks. Astrophysical Journal, 736:144, August 2011.
doi: 10.1088/0004-637X/736/2/144.
X.-N. Bai and J. M. Stone. Local Study of Accretion Disks with a Strong Vertical Magnetic Field:
Magnetorotational Instability and Disk Outflow. Astrophysical Journal, 767:30, April 2013a. doi:
10.1088/0004-637X/767/1/30.
X.-N. Bai and J. M. Stone. Wind-driven Accretion in Protoplanetary Disks. I. Suppression of the Magnetorotational Instability and Launching of the Magnetocentrifugal Wind.
Astrophysical Journal, 769:76, May 2013b. doi: 10.1088/0004-637X/769/1/76.

164
X.-N. Bai and J. M. Stone. Magnetic Flux Concentration and Zonal Flows in Magnetorotational
Instability Turbulence. Astrophysical Journal, 796:31, November 2014. doi: 10.1088/0004637X/796/1/31.
Xue-Ning Bai and James M. Stone. The Effect of the Radial Pressure Gradient in Protoplanetary
Disks on Planetesimal Formation. Astrophysical Journal, 722(2):L220–L223, Oct 2010c. doi:
10.1088/2041-8205/722/2/L220.
Xue-Ning Bai, Jiani Ye, Jeremy Goodman, and Feng Yuan. Magneto-thermal Disk Winds
from Protoplanetary Disks. Astrophysical Journal, 818(2):152, Feb 2016. doi: 10.3847/0004637X/818/2/152.
S. A. Balbus and J. F. Hawley. A powerful local shear instability in weakly magnetized disks. I- Linear analysis. II - Nonlinear evolution. Astrophysical Journal, 376:214–233, July 1991. doi:
10.1086/170270.
S. A. Balbus and J. F. Hawley. Instability, turbulence, and enhanced transport in accretion disks.
Reviews of Modern Physics, 70:1–53, January 1998. doi: 10.1103/RevModPhys.70.1.
S. A. Balbus and J. C. B. Papaloizou.
On the Dynamical Foundations of α Disks.
Astrophysical Journal, 521:650–658, August 1999. doi: 10.1086/307594.
S. A. Balbus and C. Terquem. Linear Analysis of the Hall Effect in Protostellar Disks.
Astrophysical Journal, 552:235–247, May 2001. doi: 10.1086/320452.
J. Barnes and P. Hut. A hierarchical O(N log N) force-calculation algorithm. Nature, 324:446–449,
December 1986. doi: 10.1038/324446a0.
C. Baruteau, S. Fromang, R. P. Nelson, and F. Masset. Corotation torques experienced by planets
embedded in weakly magnetized turbulent discs. Astronomy and Astrophysics, 533:A84, September 2011. doi: 10.1051/0004-6361/201117227.
K. Batygin, P. H. Bodenheimer, and G. P. Laughlin. In Situ Formation and Dynamical Evolutionof Hot Jupiter Systems. Astrophysical Journal, 829:114, October 2016. doi: 10.3847/0004637X/829/2/114.
C. Beaugé and D. Nesvorný. Multiple-planet Scattering and the Origin of Hot Jupiters.
Astrophysical Journal, 751:119, June 2012. doi: 10.1088/0004-637X/751/2/119.
T. Birnstiel, C. W. Ormel, and C. P. Dullemond.
Dust size distributions incoagulation/fragmentation equilibrium:
numerical solutions and analytical fits.
Astronomy and Astrophysics, 525:A11, January 2011. doi: 10.1051/0004-6361/201015228.
O. M. Blaes and S. A. Balbus. Local shear instabilities in weakly ionized, weakly magnetized disks.
Astrophysical Journal, 421:163–177, January 1994. doi: 10.1086/173634.
R. D. Blandford and D. G. Payne. Hydromagnetic flows from accretion disks and the production ofradio jets. Monthly Notices of the RAS, 199:883–903, Jun 1982. doi: 10.1093/mnras/199.4.883.
Erika Böhm-Vitense. Chromospheric Activity in G and K Main-Sequence Stars, and What
It Tells Us about Stellar Dynamos. Astrophysical Journal, 657(1):486–493, Mar 2007. doi:
10.1086/510482.

165
A. Brandenburg, A. Nordlund, R. F. Stein, and U. Torkelsson. Dynamo-generated Turbulence andLarge-Scale Magnetic Fields in a Keplerian Shear Flow. Astrophysical Journal, 446:741, June
1995. doi: 10.1086/175831.
K.
P.
BurnhamandD.
R.
Anderson.
Modelselectionandmultimodel
inference:A
practical
information-theoretic
approach.
In
Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach Burnham and Anderson (2002b).
K. P. Burnham and D. R. Anderson. Model Selection and Multimodel Inference: A Practical
Information-Theoretic Approach. Springer-Verlag, 2002b.
Daniel Carrera, Anders Johansen, and Melvyn B. Davies. How to form planetesimals from mmsized chondrules and chondrule aggregates. Astronomy and Astrophysics, 579:A43, Jul 2015. doi:
10.1051/0004-6361/201425120.
Luca Casagrande, Chris Flynn, and Michael Bessell. M dwarfs: effective temperatures, radiiand metallicities. Monthly Notices of the RAS, 389(2):585–607, Sep 2008. doi: 10.1111/j.13652966.2008.13573.x.
E. I. Chiang and P. Goldreich. Spectral Energy Distributions of T Tauri Stars with Passive Circumstellar Disks. Astrophysical Journal, 490:368–376, November 1997. doi: 10.1086/304869.
A. Clauset, C. R. Shalizi, and M. E. J. Newman. Power-Law Distributions in Empirical Data.
SIAM Review, 51:661–703, January 2009. doi: 10.1137/070710111.
L. I. Cleeves, E. A. Bergin, C. Qi, F. C. Adams, and K. I. Öberg. Constraining the X-Ray
and Cosmic-Ray Ionization Chemistry of the TW Hya Protoplanetary Disk: Evidence for a Subinterstellar Cosmic-Ray Rate. Astrophysical Journal, 799:204, February 2015. doi: 10.1088/0004637X/799/2/204.
P. Colella and P. R. Woodward. The Piecewise Parabolic Method (PPM) for Gas-Dynamical
Simulations. Journal of Computational Physics, 54:174–201, September 1984. doi: 10.1016/00219991(84)90143-8.
S. W. Davis, J. M. Stone, and M. E. Pessah. Sustained Magnetorotational Turbulence in Local
Simulations of Stratified Disks with Zero Net Magnetic Flux. Astrophysical Journal, 713:52–65,
April 2010. doi: 10.1088/0004-637X/713/1/52.
S. J. Desch. Linear Analysis of the Magnetorotational Instability, Including Ambipolar Diffusion,
with Application to Protoplanetary Disks. Astrophysical Journal, 608:509–525, June 2004. doi:
10.1086/392527.
C. Dominik, J. Blum, J. N. Cuzzi, and G. Wurm. Growth of Dust as the Initial Step Toward Planet
Formation. Protostars and Planets V, pages 783–800, 2007.
J. Dubinski, R. Narayan, and T. G. Phillips.
Turbulence in Molecular Clouds.
Astrophysical Journal, 448:226, July 1995. doi: 10.1086/175954.
B. Dubrulle, G. Morfill, and M. Sterzik. The dust subdisk in the protoplanetary nebula. Icarus,
114:237–246, April 1995. doi: 10.1006/icar.1995.1058.

166
C. P. Dullemond and C. Dominik. Dust coagulation in protoplanetary disks: A rapid depletion of small grains. Astronomy and Astrophysics, 434:971–986, May 2005. doi: 10.1051/00046361:20042080.
B.

Efron and R. J. Tibshirani.
An introduction to
An Introduction to the Bootstrap Efron and Tibshirani (1994b).

thebootstrap.

In

B. Efron and R. J. Tibshirani. An Introduction to the Bootstrap. Chapman & Hall/CRC, 1994b.
D. J. Eisenstein and P. Hut. HOP: A New Group-Finding Algorithm for N-Body Simulations.
Astrophysical Journal, 498:137–142, May 1998. doi: 10.1086/305535.
C. R. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flows - A constrained transport
method. Astrophysical Journal, 332:659–677, September 1988a. doi: 10.1086/166684.
C. R. Evans and J. F. Hawley. Simulation of magnetohydrodynamic flows - A constrained transport
method. Astrophysical Journal, 332:659–677, September 1988b. doi: 10.1086/166684.
Kyle Gerard Felker and James M. Stone. A fourth-order accurate finite volume method for ideal
MHD via upwind constrained transport. Journal of Computational Physics, 375:1365–1400, Dec
2018. doi: 10.1016/j.jcp.2018.08.025.
C. Fendt, M. Camenzind, and S. Appl. On the collimation of stellar magnetospheres to jets. I.
Relativistic force-free 2D equilibrium. Astronomy and Astrophysics, 300:791, August 1995.
K. M. Flaherty, A. M. Hughes, K. A. Rosenfeld, S. M. Andrews, E. Chiang, J. B. Simon, S. Kerzner,and D. J. Wilner. Weak Turbulence in the HD 163296 Protoplanetary Disk Revealed by
ALMA CO Observations. Astrophysical Journal, 813:99, November 2015. doi: 10.1088/0004637X/813/2/99.
K. M. Flaherty, A. M. Hughes, S. C. Rose, J. B. Simon, C. Qi, S. M. Andrews, Á. Kóspál,
D. J. Wilner, E. Chiang, P. J. Armitage, and X.-n. Bai. A Three-dimensional View of Turbulence: Constraints on Turbulent Motions in the HD 163296 Protoplanetary Disk Using DCO+ .
Astrophysical Journal, 843:150, July 2017. doi: 10.3847/1538-4357/aa79f9.
K. M. Flaherty, A. M. Hughes, R. Teague, J. B. Simon, S. M. Andrews, and D. J. Wilner. Turbulence in the TW Hya Disk. Astrophysical Journal, 856:117, April 2018. doi: 10.3847/15384357/aab615.
T. Fleming and J. M. Stone. Local Magnetohydrodynamic Models of Layered Accretion Disks.
Astrophysical Journal, 585:908–920, March 2003. doi: 10.1086/345848.
Wesley C. Fraser, Michael E. Brown, and Megan E. Schwamb. The luminosity functionof the hot and cold Kuiper belt populations. Icarus, 210(2):944–955, Dec 2010. doi:
10.1016/j.icarus.2010.08.001.
Wesley C. Fraser, Michael E. Brown, Alessand ro Morbidelli, Alex Parker, and Konstantin Batygin.
The Absolute Magnitude Distribution of Kuiper Belt Objects. Astrophysical Journal, 782(2):100,
Feb 2014. doi: 10.1088/0004-637X/782/2/100.
S. Fromang, C. Terquem, and S. A. Balbus. The ionization fraction in α models of protoplanetary discs. Monthly Notices of the RAS, 329:18–28, January 2002. doi: 10.1046/j.13658711.2002.04940.x.

167
C. F. Gammie. Layered Accretion in T Tauri Disks. Astrophysical Journal, 457:355, January 1996.
doi: 10.1086/176735.
C. F. Gammie and K. Menou.
On the Origin of Episodic Accretion in Dwarf Novae.
Astrophysical Journal, Letters, 492:L75–L78, January 1998. doi: 10.1086/311091.
T. A. Gardiner and J. M. Stone. An unsplit Godunov method for ideal MHD via constrained transport. Journal of Computational Physics, 205:509–539, May 2005. doi: 10.1016/j.jcp.2004.11.016.
T. A. Gardiner and J. M. Stone. An unsplit Godunov method for ideal MHD via constrainedtransport in three dimensions. Journal of Computational Physics, 227:4123–4141, April 2008.
doi: 10.1016/j.jcp.2007.12.017.
P. Ghosh and F. K. Lamb. Accretion by rotating magnetic neutron stars. III. Accretion torques
and period changes in pulsating X-ray sources. Astrophysical Journal, 234:296–316, Nov 1979.
doi: 10.1086/157498.
D. Gole, J. B. Simon, S. H. Lubow, and P. J. Armitage. Turbulence, Transport, and Waves in Ohmic
Dead Zones. Astrophysical Journal, 826:18, July 2016. doi: 10.3847/0004-637X/826/1/18.
O. Gressel, R. P. Nelson, and N. J. Turner. On the dynamics of planetesimals embedded in turbulent
protoplanetary discs with dead zones. Monthly Notices of the RAS, 415:3291–3307, August 2011.
doi: 10.1111/j.1365-2966.2011.18944.x.
O. Gressel, R. P. Nelson, and N. J. Turner. Dead zones as safe havens for planetesimals: influence ofdisc mass and external magnetic field. Monthly Notices of the RAS, 422:1140–1159, May 2012.
doi: 10.1111/j.1365-2966.2012.20701.x.
X. Guan, C. F. Gammie, J. B. Simon, and B. M. Johnson. Locality of MHD Turbulencein Isothermal Disks. Astrophysical Journal, 694:1010–1018, April 2009. doi: 10.1088/0004637X/694/2/1010.
J. Guilet and G. I. Ogilvie. Global evolution of the magnetic field in a thin disc and its consequences for protoplanetary systems. Monthly Notices of the RAS, 441:852–868, June 2014. doi:
10.1093/mnras/stu532.
K. E. Haisch, Jr., E. A. Lada, and C. J. Lada. Disk Frequencies and Lifetimes in Young Clusters.
Astrophysical Journal, Letters, 553:L153–L156, June 2001. doi: 10.1086/320685.
L.

HartmannandS.
J.
Kenyon.
The
FU
Annual Review of Astron and Astrophys,
34:207–240,
1996.
nurev.astro.34.1.207.

Orionisdoi:

Phenomenon.
10.1146/an-

J. F. Hawley and J. M. Stone. Nonlinear Evolution of the Magnetorotational Instability in IonNeutral Disks. Astrophysical Journal, 501:758–771, July 1998. doi: 10.1086/305849.
J. F. Hawley, C. F. Gammie, and S. A. Balbus. Local Three-dimensional Magnetohydrodynamic Simulations of Accretion Disks. Astrophysical Journal, 440:742, February 1995. doi:
10.1086/175311.
C. Hayashi. Structure of the Solar Nebula, Growth and Decay of Magnetic Fields and Effects ofMagnetic and Turbulent Viscosities on the Nebula. Progress of Theoretical Physics Supplement,
70:35–53, 1981. doi: 10.1143/PTPS.70.35.

168
G. H. Herbig. Eruptive phenomena in early stellar evolution. Astrophysical Journal, 217:693–715,
November 1977. doi: 10.1086/155615.
G. H. Herbig. FU Orionis eruptions. In B. Reipurth, editor, European Southern Observatory
Conference and Workshop Proceedings, volume 33 of European Southern Observatory Conference
and Workshop Proceedings, pages 233–246, September 1989.
G. H. Herbig. History and Spectroscopy of EXor Candidates. Astronomical Journal, 135:637–648,
February 2008. doi: 10.1088/0004-6256/135/2/637.
P. F. Hopkins. A New Public Release of the GIZMO Code. arXiv e-prints, December 2017.
X. Hu, Z. Zhu, S. Okuzumi, X.-N. Bai, L. Wang, K. Tomida, and J. M. Stone. Non-ideal MHD
simulation of HL Tau disk: formation of rings. arXiv e-prints, April 2019.
J. Huang, S. M. Andrews, L. I. Cleeves, K. I. Öberg, D. J. Wilner, X. Bai, T. Birnstiel, J. Carpenter,A. M. Hughes, A. Isella, L. M. Pérez, L. Ricci, and Z. Zhu. CO and Dust Properties in the TW
Hya Disk from High-resolution ALMA Observations. Astrophysical Journal, 852:122, January
2018. doi: 10.3847/1538-4357/aaa1e7.
S. Ida, T. Guillot, and A. Morbidelli. Accretion and Destruction of Planetesimals in Turbulent
Disks. Astrophysical Journal, 686:1292–1301, October 2008. doi: 10.1086/591903.
J. Igea and A. E. Glassgold.
X-Ray Ionization of the Disks of Young Stellar Objects.
Astrophysical Journal, 518:848–858, June 1999. doi: 10.1086/307302.
M. Ilgner and R. P. Nelson. On the ionisation fraction in protoplanetary disks. I. Comparing
different reaction networks. Astronomy and Astrophysics, 445:205–222, January 2006a. doi:
10.1051/0004-6361:20053678.
M. Ilgner and R. P. Nelson. On the ionisation fraction in protoplanetary disks. II. The effect of
turbulent mixing on gas-phase chemistry. Astronomy and Astrophysics, 445:223–232, January
2006b. doi: 10.1051/0004-6361:20053867.
M. Ilgner and R. P. Nelson. On the ionisation fraction in protoplanetary disks. III. The effect of
X-ray flares on gas-phase chemistry. Astronomy and Astrophysics, 455:731–740, August 2006c.
doi: 10.1051/0004-6361:20065308.
A. F. Illarionov and D. A. Kompaneets. A Spin-Down Mechanism for Accreting Neutron Stars.
Monthly Notices of the RAS, 247:219, Nov 1990.
A. Johansen and A. Youdin. Protoplanetary Disk Turbulence Driven by the Streaming Instability: Nonlinear Saturation and Particle Concentration. Astrophysical Journal, 662:627–641, June
2007. doi: 10.1086/516730.
A. Johansen, A. Youdin, and H. Klahr. Zonal Flows and Long-lived Axisymmetric Pressure
Bumps in Magnetorotational Turbulence. Astrophysical Journal, 697:1269–1289, June 2009. doi:
10.1088/0004-637X/697/2/1269.
A. Johansen, A. N. Youdin, and Y. Lithwick. Adding particle collisions to the formation of asteroids
and Kuiper belt objects via streaming instabilities. Astronomy and Astrophysics, 537:A125,
January 2012. doi: 10.1051/0004-6361/201117701.

169
Christopher M. Johns-Krull. Measuring T Tauri star magnetic fields. In Klaus G. Strassmeier,Alexander G. Kosovichev, and John E. Beckman, editors, Cosmic Magnetic Fields: From
Planets, to Stars and Galaxies, volume 259 of IAU Symposium, pages 345–356, Apr 2009. doi:
10.1017/S1743921309030713.
A. H. Joy. T Tauri Variable Stars.
10.1086/144749.

Astrophysical Journal, 102:168, September 1945.

doi:

Robert E. Kass and Adrian E. Raftery.
Bayes factors.
Journal of the American
Statistical Association, 90(430):773–795, 1995. doi: 10.1080/01621459.1995.10476572. URL
https://amstat.tandfonline.com/doi/abs/10.1080/01621459.1995.10476572.
A. Koenigl. Disk accretion onto magnetic T Tauri stars. Astrophysical Journal, Letters, 370:
L39–L43, March 1991. doi: 10.1086/185972.
A. Kolmogorov. The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large
Reynolds’ Numbers. Akademiia Nauk SSSR Doklady, 30:301–305, 1941.
H. Koyama and E. C. Ostriker. Gas Properties and Implications for Galactic Star Formation in
Numerical Models of the Turbulent, Multiphase Interstellar Medium. Astrophysical Journal, 693:
1316–1345, March 2009. doi: 10.1088/0004-637X/693/2/1316.
M. W. Kunz. On the linear stability of weakly ionized, magnetized planar shear flows.
Monthly Notices of the RAS, 385:1494–1510, April 2008. doi: 10.1111/j.1365-2966.2008.12928.x.
M. W. Kunz and S. A. Balbus.
Ambipolar diffusion in the magnetorotational instability.
Monthly Notices of the RAS, 348:355–360, February 2004.
doi: 10.1111/j.13652966.2004.07383.x.
C. J. Lada and B. A. Wilking. The nature of the embedded population in the Rho Ophiuchi darkcloud - Mid-infrared observations. Astrophysical Journal, 287:610–621, December 1984. doi:
10.1086/162719.
M. Lambrechts and A. Johansen. Rapid growth of gas-giant cores by pebble accretion.
Astronomy and Astrophysics, 544:A32, August 2012. doi: 10.1051/0004-6361/201219127.
M. N. Lemaster and J. M. Stone. Power Spectra and Density Probability Distribution Functionsin Supersonic Hydrodynamic and MHD Turbulence. In Revista Mexicana de Astronomia y
Astrofisica Conference Series, volume 36 of Revista Mexicana de Astronomia y Astrofisica, vol. 27,
pages CD243–CD251, August 2009.
G. Lesur, M. W. Kunz, and S. Fromang. Thanatology in protoplanetary discs. The combined
influence of Ohmic, Hall, and ambipolar diffusion on dead zones. Astronomy and Astrophysics,
566:A56, June 2014. doi: 10.1051/0004-6361/201423660.
R. Li, A. N. Youdin, and J. B. Simon. On the Numerical Robustness of the Streaming Instability:
Particle Concentration and Gas Dynamics in Protoplanetary Disks. Astrophysical Journal, 862:
14, July 2018. doi: 10.3847/1538-4357/aaca99.
J. J. Lissauer and G. R. Stewart. Growth of planets from planetesimals. In E. H. Levy and J. I.
Lunine, editors, Protostars and Planets III, pages 1061–1088, 1993.

170
R. V. E. Lovelace, M. M. Romanova, and G. S. Bisnovatyi-Kogan. Spin-up/spin-down of magnetizedstars with accretion discs and outflows. Monthly Notices of the RAS, 275:244–254, July 1995.
doi: 10.1093/mnras/275.2.244.
S. H. Lubow and S. Ida. Planet Migration. In Seager (2010), pages 347–371.
S. H. Lubow and J. E. Pringle. Wave propagation in accretion disks - Axisymmetric case.
Astrophysical Journal, 409:360–371, May 1993. doi: 10.1086/172669.
S. H. Lubow, J. C. B. Papaloizou, and J. E. Pringle. Magnetic field dragging in accretion discs.
Monthly Notices of the RAS, 267:235–240, March 1994. doi: 10.1093/mnras/267.2.235.
D. Lynden-Bell and C. Boily. Self-Similar Solutions up to Flashpoint in Highly Wound Magnetostatics. Monthly Notices of the RAS, 267:146, March 1994. doi: 10.1093/mnras/267.1.146.
R. G. Martin and S. H. Lubow.
The Gravo-magneto Limit Cycle in Accretion Disks.
Astrophysical Journal, Letters, 740:L6, October 2011. doi: 10.1088/2041-8205/740/1/L6.
R. G. Martin and S. H. Lubow. The gravo-magneto disc instability with a viscous dead zone.
Monthly Notices of the RAS, 437:682–689, January 2014. doi: 10.1093/mnras/stt1917.
F. S. Masset and J. C. B. Papaloizou. Runaway Migration and the Formation of Hot Jupiters.
Astrophysical Journal, 588:494–508, May 2003. doi: 10.1086/373892.
K. Menou and E. Quataert. Ionization, Magnetorotational, and Gravitational Instabilities in Thin
Accretion Disks Around Supermassive Black Holes. Astrophysical Journal, 552:204–208, May
2001. doi: 10.1086/320466.
A. Mignone. A simple and accurate Riemann solver for isothermal MHD. Journal of Computational
Physics, 225:1427–1441, August 2007. doi: 10.1016/j.jcp.2007.01.033.
K. A. Miller and J. M. Stone. The Formation and Structure of a Strongly Magnetized Corona
above a Weakly Magnetized Accretion Disk. Astrophysical Journal, 534:398–419, May 2000. doi:
10.1086/308736.
T. Miyoshi and K. Kusano. A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics. Journal of Computational Physics, 208:315–344, September 2005. doi:
10.1016/j.jcp.2005.02.017.
Alessandro Morbidelli, William F. Bottke, David Nesvorný, and Harold F. Levison. Asteroids wereborn big. Icarus, 204(2):558–573, Dec 2009. doi: 10.1016/j.icarus.2009.07.011.
J. Muzerolle, P. D’Alessio, N. Calvet, and L. Hartmann. Magnetospheres and Disk Accretion in
Herbig Ae/Be Stars. Astrophysical Journal, 617:406–417, December 2004. doi: 10.1086/425260.
R. P. Nelson and J. C. B. Papaloizou. The interaction of giant planets with a disc with MHD
turbulence - IV. Migration rates of embedded protoplanets. Monthly Notices of the RAS, 350:
849–864, May 2004. doi: 10.1111/j.1365-2966.2004.07406.x.
G. I. Ogilvie. Waves and instabilities in a differentially rotating disc containing a poloidal
magnetic field. Monthly Notices of the RAS, 297:291–314, June 1998. doi: 10.1046/j.13658711.1998.01507.x.

171
G. I. Ogilvie.
Jet launching from accretion discs in the local approximation.
Monthly Notices of the RAS, 423:1318–1324, June 2012. doi: 10.1111/j.1365-2966.2012.20958.x.
J. S. Oishi and M.-M. Mac Low. On Hydrodynamic Motions in Dead Zones. Astrophysical Journal,
704:1239–1250, October 2009. doi: 10.1088/0004-637X/704/2/1239.
S. Okuzumi and S. Hirose. Modeling Magnetorotational Turbulence in Protoplanetary Disks with
Dead Zones. Astrophysical Journal, 742:65, December 2011. doi: 10.1088/0004-637X/742/2/65.
S. Okuzumi and C. W. Ormel. The Fate of Planetesimals in Turbulent Disks with Dead Zones.
I. The Turbulent Stirring Recipe. Astrophysical Journal, 771:43, July 2013. doi: 10.1088/0004637X/771/1/43.
S. Okuzumi, T. Takeuchi, and T. Muto. Radial Transport of Large-scale Magnetic Fields inAccretion Disks. I. Steady Solutions and an Upper Limit on the Vertical Field Strength.
Astrophysical Journal, 785:127, April 2014. doi: 10.1088/0004-637X/785/2/127.
S. Okuzumi, M. Momose, S.-i. Sirono, H. Kobayashi, and H. Tanaka. Sintering-induced Dust Ring
Formation in Protoplanetary Disks: Application to the HL Tau Disk. Astrophysical Journal,
821:82, April 2016. doi: 10.3847/0004-637X/821/2/82.
C. W. Ormel and J. N. Cuzzi. Closed-form expressions for particle relative velocities inducedby turbulence. Astronomy and Astrophysics, 466:413–420, May 2007. doi: 10.1051/00046361:20066899.
E. C. Ostriker and F. H. Shu. Magnetocentrifugally Driven Flows from Young Stars and Disks.
IV. The Accretion Funnel and Dead Zone. Astrophysical Journal, 447:813, July 1995. doi:
10.1086/175920.
S.-J. Paardekooper, C. Baruteau, and W. Kley. A torque formula for non-isothermal Type I
planetary migration - II. Effects of diffusion. Monthly Notices of the RAS, 410:293–303, January
2011. doi: 10.1111/j.1365-2966.2010.17442.x.
E. S. Parsamian and R. Mujica. Spectral Observations of Subfuors (Exors). II. V1143 Ori.
Astrophysics, 47:433–442, October 2004. doi: 10.1023/B:ASYS.0000049780.78084.40.
D. Perez-Becker and E. Chiang. Surface Layer Accretion in Conventional and Transitional
Disks Driven by Far-ultraviolet Ionization. Astrophysical Journal, 735:8, July 2011a. doi:
10.1088/0004-637X/735/1/8.
D. Perez-Becker and E. Chiang. Surface Layer Accretion in Conventional and Transitional
Disks Driven by Far-ultraviolet Ionization. Astrophysical Journal, 735:8, July 2011b. doi:
10.1088/0004-637X/735/1/8.
D. Perez-Becker and E. Chiang. Surface Layer Accretion in Conventional and Transitional
Disks Driven by Far-ultraviolet Ionization. Astrophysical Journal, 735:8, July 2011c. doi:
10.1088/0004-637X/735/1/8.
C. Pinte, W. R. F. Dent, F. Ménard, A. Hales, T. Hill, P. Cortes, and I. de Gregorio-Monsalvo.
Dust and Gas in the Disk of HL Tauri: Surface Density, Dust Settling, and Dust-to-gas Ratio.
Astrophysical Journal, 816:25, January 2016. doi: 10.3847/0004-637X/816/1/25.

172
G. Salvesen, J. B. Simon, P. J. Armitage, and M. C. Begelman. Accretion disc dynamo
activity in local simulations spanning weak-to-strong net vertical magnetic flux regimes.
Monthly Notices of the RAS, 457:857–874, March 2016. doi: 10.1093/mnras/stw029.
T. Sano and S. M. Miyama. Magnetorotational Instability in Protoplanetary Disks. I. On the
Global Stability of Weakly Ionized Disks with Ohmic Dissipation. Astrophysical Journal, 515:
776–786, April 1999. doi: 10.1086/307063.
T. Sano and J. M. Stone. The Effect of the Hall Term on the Nonlinear Evolution of the Magnetorotational Instability. I. Local Axisymmetric Simulations. Astrophysical Journal, 570:314–328,
May 2002. doi: 10.1086/339504.
U. Schäfer, C.-C. Yang, and A. Johansen. Initial mass function of planetesimals formed by thestreaming instability. Astronomy and Astrophysics, 597:A69, January 2017. doi: 10.1051/00046361/201629561.
W. Schmitt, T. Henning, and R. Mucha. Dust evolution in protoplanetary accretion disks.
Astronomy and Astrophysics, 325:569–584, September 1997.
K. R. Schwarz, E. A. Bergin, L. I. Cleeves, G. A. Blake, K. Zhang, K. I. Öberg, E. F. van Dishoeck,
and C. Qi. The Radial Distribution of H2 and CO in TW Hya as Revealed by Resolved ALMA
Observations of CO Isotopologues. Astrophysical Journal, 823:91, June 2016. doi: 10.3847/0004637X/823/2/91.
S. Seager, editor. Exoplanets. University of Arizona Press, 2010.
D. Semenov, T. Henning, C. Helling, M. Ilgner, and E. Sedlmayr. Rosseland and Planck mean
opacities for protoplanetary discs. Astronomy and Astrophysics, 410:611–621, November 2003.
doi: 10.1051/0004-6361:20031279.
N. I. Shakura and R. A. Sunyaev. Black holes in binary systems. Observational appearance.
Astronomy and Astrophysics, 24:337–355, 1973.
J. B. Simon, J. F. Hawley, and K. Beckwith. Simulations of Magnetorotational Turbulence witha Higher-Order Godunov Scheme. Astrophysical Journal, 690:974–997, January 2009. doi:
10.1088/0004-637X/690/1/974.
J. B. Simon, K. Beckwith, and P. J. Armitage. Emergent mesoscale phenomena in magnetized accretion disc turbulence. Monthly Notices of the RAS, 422:2685–2700, May 2012. doi:
10.1111/j.1365-2966.2012.20835.x.
J. B. Simon, X.-N. Bai, P. J. Armitage, J. M. Stone, and K. Beckwith. Turbulence in the Outer
Regions of Protoplanetary Disks. II. Strong Accretion Driven by a Vertical Magnetic Field.
Astrophysical Journal, 775:73, September 2013a. doi: 10.1088/0004-637X/775/1/73.
J. B. Simon, X.-N. Bai, J. M. Stone, P. J. Armitage, and K. Beckwith. Turbulence in theOuter Regions of Protoplanetary Disks. I. Weak Accretion with No Vertical Magnetic Flux.
Astrophysical Journal, 764:66, February 2013b. doi: 10.1088/0004-637X/764/1/66.
J. B. Simon, A. M. Hughes, K. M. Flaherty, X.-N. Bai, and P. J. Armitage. Signatures of MRI-driven
Turbulence in Protoplanetary Disks: Predictions for ALMA Observations. Astrophysical Journal,
808:180, August 2015a. doi: 10.1088/0004-637X/808/2/180.

173
J. B. Simon, G. Lesur, M. W. Kunz, and P. J. Armitage. Magnetically driven accretion in protoplanetary discs. Monthly Notices of the RAS, 454:1117–1131, November 2015b. doi: 10.1093/mnras/stv2070.
J. B. Simon, P. J. Armitage, R. Li, and A. N. Youdin. The Mass and Size Distribution of Planetesimals Formed by the Streaming Instability. I. The Role of Self-gravity. Astrophysical Journal,
822:55, May 2016. doi: 10.3847/0004-637X/822/1/55.
J. B. Simon, X.-N. Bai, K. M. Flaherty, and A. M. Hughes. A New Model for Weak Turbulence inProtoplanetary Disks. ArXiv e-prints, November 2017a.
J. B. Simon, X.-N. Bai, K. M. Flaherty, and A. M. Hughes. Origin of Weak Turbulence in theOuter Regions of Protoplanetary Disks. Astrophysical Journal, 865:10, September 2018. doi:
10.3847/1538-4357/aad86d.
Jacob B. Simon, Philip J. Armitage, Andrew N. Youdin, and Rixin Li. Evidence for Universality
in the Initial Planetesimal Mass Function. Astrophysical Journal, 847(2):L12, Oct 2017b. doi:
10.3847/2041-8213/aa8c79.
S. T. Stewart and Z. M. Leinhardt. Velocity-Dependent Catastrophic Disruption Criteria for Planetesimals. Astrophysical Journal, Letters, 691:L133–L137, February 2009. doi: 10.1088/0004637X/691/2/L133.
J. M. Stone and T. A. Gardiner. Implementation of the Shearing Box Approximation in Athena.
Astrophysical Journal, Supplement, 189:142–155, July 2010. doi: 10.1088/0067-0049/189/1/142.
J. M. Stone, J. F. Hawley, C. F. Gammie, and S. A. Balbus. Three-dimensional Magnetohydrodynamical Simulations of Vertically Stratified Accretion Disks. Astrophysical Journal, 463:656,
June 1996. doi: 10.1086/177280.
J. M. Stone, T. A. Gardiner, P. Teuben, J. F. Hawley, and J. B. Simon. Athena: A New Codefor Astrophysical MHD. Astrophysical Journal, Supplement, 178:137–177, September 2008. doi:
10.1086/588755.
T. Takeuchi and D. N. C. Lin.
Radial Flow of Dust Particles in Accretion Disks.
Astrophysical Journal, 581:1344–1355, December 2002. doi: 10.1086/344437.
C. E. J. M. L. J. Terquem. New Composite Models of Partially Ionized Protoplanetary Disks.
Astrophysical Journal, 689:532–538, December 2008. doi: 10.1086/592597.
N. J. Turner and T. Sano.
Dead Zone Accretion Flows in Protostellar Disks.
Astrophysical Journal, Letters, 679:L131, June 2008. doi: 10.1086/589540.
N. J. Turner, T. Sano, and N. Dziourkevitch. Turbulent Mixing and the Dead Zone in Protostellar
Disks. Astrophysical Journal, 659:729–737, April 2007. doi: 10.1086/512007.
N. J. Turner, S. Fromang, C. Gammie, H. Klahr, G. Lesur, M. Wardle, and X.-N. Bai. Transportand Accretion in Planet-Forming Disks. Protostars and Planets VI, pages 411–432, 2014. doi:
10.2458/azu uapress 9780816531240-ch018.
T. Umebayashi. The Densities of Charged Particles in Very Dense Interstellar Clouds. Progress ofTheoretical Physics, 69:480–502, February 1983. doi: 10.1143/PTP.69.480.

174
A. A. van Ballegooijen. Energy release in stellar magnetospheres. Space Science Reviews, 68(1-4):
299–307, May 1994. doi: 10.1007/BF00749156.
A. A. Vidotto, M. Jardine, J. Morin, J. F. Donati, P. Lang, and A. J. B. Russell. Effects of M
dwarf magnetic fields on potentially habitable planets. Astronomy and Astrophysics, 557:A67,
Sep 2013. doi: 10.1051/0004-6361/201321504.
M. Wardle. The Balbus-Hawley instability in weakly ionized discs. Monthly Notices of the RAS,
307:849–856, August 1999. doi: 10.1046/j.1365-8711.1999.02670.x.
M. Wardle. Magnetic fields in protoplanetary disks. Astronomy and Astrophysics, Supplement,
311:35–45, October 2007. doi: 10.1007/s10509-007-9575-8.
S. J. Weidenschilling.
The distribution of mass in the planetary system and solar
nebula.
Astronomy and Astrophysics, Supplement, 51:153–158, September 1977.
doi:
10.1007/BF00642464.
S. J. Weidenschilling. Dust to planetesimals - Settling and coagulation in the solar nebula. Icarus,
44:172–189, October 1980. doi: 10.1016/0019-1035(80)90064-0.
J. N. Winn and D. C. Fabrycky. The Occurrence and Architecture of Exoplanetary Systems.
Annual Review of Astron and Astrophys, 53:409–447, August 2015. doi: 10.1146/annurev-astro082214-122246.
G. Wurm, G. Paraskov, and O. Krauss. Growth of planetesimals by impacts at 25 m/s. Icarus,
178:253–263, November 2005. doi: 10.1016/j.icarus.2005.04.002.
R. Xu and X.-N. Bai. On the Grain-modified Magnetic Diffusivities in Protoplanetary Disks.
Astrophysical Journal, 819:68, March 2016. doi: 10.3847/0004-637X/819/1/68.
A. N. Youdin and J. Goodman.
Streaming Instabilities in Protoplanetary Disks.
Astrophysical Journal, 620:459–469, February 2005. doi: 10.1086/426895.
A. N. Youdin and Y. Lithwick. Particle stirring in turbulent gas disks: Including orbital oscillations.
Icarus, 192:588–604, December 2007. doi: 10.1016/j.icarus.2007.07.012.
Z. Zhu and J. M. Stone. Global Evolution of an Accretion Disk with a Net Vertical Field: Coronal
Accretion, Flux Transport, and Disk Winds. Astrophysical Journal, 857:34, April 2018. doi:
10.3847/1538-4357/aaafc9.
Z. Zhu, L. Hartmann, N. Calvet, J. Hernandez, J. Muzerolle, and A.-K. Tannirkulam. The Hot Inner
Disk of FU Orionis. Astrophysical Journal, 669:483–492, November 2007. doi: 10.1086/521345.
Z. Zhu, L. Hartmann, N. Calvet, J. Hernandez, A.-K. Tannirkulam, and P. D’Alessio. LongWavelength Excesses of FU Orionis Objects: Flared Outer Disks or Infalling Envelopes?
Astrophysical Journal, 684:1281–1290, September 2008. doi: 10.1086/590241.
Z. Zhu, L. Hartmann, and C. Gammie. Nonsteady Accretion in Protostars. Astrophysical Journal,
694:1045–1055, April 2009. doi: 10.1088/0004-637X/694/2/1045.
Z. Zhu, J. M. Stone, and X.-N. Bai. Dust Transport in MRI Turbulent Disks: Ideal and Non-Ideal
MHD With Ambipolar Diffusion. Astrophysical Journal, 801:81, March 2015. doi: 10.1088/0004637X/801/2/81.

175
A. Zsom, C. W. Ormel, C. Güttler, J. Blum, and C. P. Dullemond. The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals? II. Introducing the bouncing barrier.
Astronomy and Astrophysics, 513:A57, April 2010. doi: 10.1051/0004-6361/200912976.