Master Thesis

The Eﬀects of Electron Heatingon the Magnetorotational
Instability
in Protoplanetary Disks

Shoji Mori

Okuzumi Laboratory
Department of Earth and Planetary Science
Tokyo Institute of Technology
October 28, 2016

Abstract
To elucidate mysterious planet formation, understanding protoplanetary disks is essential. Especially, turbulence in the disk greatly aﬀects dynamics of solid particles whose
size is from sub-micron to kilometers. Magnetorotaional instability (MRI) is thought
to be the most plausible mechanism generating the vigorous turbulence. However,
although many studies of MRI have addressed the turbulence strength, the strength
remains debatable mainly because of the uncertainty of MRI in the low ionization fraction. Thus, the role of MRI in planet formation is still a puzzle.
So far we have investigated the eﬀect of electron heating on protoplanetary disks. The
electron heating takes place when the electric field induced by the magnetic turbulenceheats up electrons. The heated electrons frequently collide with and stick to dust grains,which in turn decreases the ionization fraction. MRI is stabilized under the suﬃciently
low ionization faction, which indicates high resistivity. Our previous work showedthat the electron heating might quench the magnetic turbulence in extensive regions.
However, how much the electron heating suppresses the turbulence is unknown.
To answer the question, we first numerically investigate the eﬀect of electron heatingon the nonlinear evolution of MRI. We perform magnetohydrodynamical simulations
including increase of resistivity by electron heating. We introduce a simple analytic
resistivity model that enables us to treat the increasing resistivity.
Our simulations confirm the electron heating suppresses magnetic turbulence. We finda clear relation between magnetic turbulence strength and its current density. This
relation means that a lower current leads to lower turbulence strength.. We find thatwhen turbulence completely dies away, laminar accretion flow is caused by orderedmagnetic field. By solving the dispersion relation, we find an analytical expression of
the laminar state that exactly gives the physical quantities. Based on the simulation
results and the scaling relation between the accretion stress and current density, we
obtain a formula that successfully predicts the accretion stress in the presence of electron
heating once current density is given. In a protoplanetary disk, the current density canbe estimated by using the current-dependent resistivity and the saturated resistivity.
Thus, we can predict the accretion stress in extensive regions where the electron heating
occur. This helps us to construct realistic planet formation theory.
The simulations presented in this thesis are based on the simple resistivity model. In
future work, we will employ the more realistic resistivity model.

Contents
1 Introduction

5

2 Background

11

2.1

2.2

Magnetohydrodynamics (MHD) and Magnetorotaional Instability (MRI)

11

2.1.1

Ideal MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.1.2

Linear Stability Analysis of MRI . . . . . . . . . . . . . . . . .

17

2.1.3

Nonideal MHD Eﬀects . . . . . . . . . . . . . . . . . . . . . . .

23

Nonlinear Ohm’s Law . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

2.2.1

Outline of Calculation . . . . . . . . . . . . . . . . . . . . . . .

24

2.2.2

Plasma Heating . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.2.3

Charge Reactions . . . . . . . . . . . . . . . . . . . . . . . . . .

27

2.2.4

Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

3 Electron Heating Zone : Mori & Okuzumi (2016)
3.1

3.2

3.3

31

Disk and Ionization Models . . . . . . . . . . . . . . . . . . . . . . . .

32

3.1.1

Disk Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

3.1.2

Simplified Plasma Heating Model . . . . . . . . . . . . . . . . .

34

3.1.3

Ionization Balance and Accuracy of Simplified Approach . . . .

37

Active, Dead, and E-heating Zone . . . . . . . . . . . . . . . . . . . . .

40

3.2.1

Conditions for MRI Growth . . . . . . . . . . . . . . . . . . . .

40

3.2.2

Zoning Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

Location of the E-heating Zone . . . . . . . . . . . . . . . . . . . . . .

43

3.3.1

Fiducial Disk Model . . . . . . . . . . . . . . . . . . . . . . . .

45

3.3.2

Dependence on the Magnetic Field Strength . . . . . . . . . . .

47

3.3.3

Dependence on the Grain Size and Dust-to-Gas Mass Ratio

47

3

. .

4

CONTENTS
3.3.4

Dependence on the Disk Mass . . . . . . . . . . . . . . . . . . .

49

3.3.5

Ion Heating . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

3.4

Saturation of Turbulence in the E-heating Zone . . . . . . . . . . . . .

51

3.5

Charge Barrier against Dust Growth in the E-heating Zone . . . . . . .

54

3.6

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

3.6.1

Dust Diﬀusion. . . . . . . . . . . . . . . . . . . . . . . . . . .

60

3.6.2

Eﬀects of Grain Size Distribution and Porosity . . . . . . . . . .

61

3.6.3

Hall Eﬀect and Ambipolar Diﬀusion . . . . . . . . . . . . . . . .

62

4 Suppression of MRI by the Electron Heating

65

4.1

Numerical Setup and Procedure . . . . . . . . . . . . . . . . . . . . . .

66

4.2

Results of Shearing-Box MHD Simulations . . . . . . . . . . . . . . . .

72

4.3

Predictive Formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

4.3.1

Scaling Relation . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

4.3.2

Analytical Solution of Laminar State . . . . . . . . . . . . . . .

83

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88

4.4.1

Vertical Diﬀusion of Dust . . . . . . . . . . . . . . . . . . . . .

88

4.4.2

Weak Disk Turbulence Suggested by Recent Observations . . . .

89

4.4

5 Summary and Conclusion

93

5.1

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

5.2

Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

94

A Appendix
A.1 Time Evolution of Simulation Results . . . . . . . . . . . . . . . . . . .

105
105

Chapter 1
Introduction
Many exoplanets have been discovered since Mayor & Queloz have discovered the first
exoplanet. At October 28, 2016, 2049 exoplanets including candidates have been found(http://exoplanet.eu). As exoplanets had been discovered, the diversity of exoplanets was revealed. The first exoplanet is a typical jupiter-like planet. However, thediscovery was surprised because the planet is much closer to the central star than researchers had expected, and thus such planets are called a hot jupiter. In addition, hot
Neptune, a neptune-like planet orbiting close to the central star, and a super earth,which consists mainly of rock but has too larger mass than the Earth, were discovered.
Although many exoplanets have been discovered, how the diversity of exoplanets originates is vague. In order to understand the diversity, the process of planet formation
have to be revealed. Planets forms in protoplanetary disks consisting of gas and dust.
The outline of dust coagulation model of classical planet formation theory is the following. First of all interstellar molecular cloud collapses with own self-gravity in thestar formation region, and then a protester and a disk surrounding the star are formed.
In the disk, the dust particles grow via collisions and sticking. The micron sized dust
particles grow up into planetesimals of which size is approximately 1 km. After the
disk disperses owing to accretion into the central star, a system of the protoplanet become unstable, and then their orbit cross. Via the impacts between protoplanets, they
become rocky planets. When mass of protoplanets excesses a critical mass, gas in the5

6

CHAPTER 1. INTRODUCTION

disk accrete into the protoplanet, and then the gas planet is formed.
However, the classical scenario has serious diﬃculties in the formation of planetesimals
from dust particles. One of the dificulties is the “meter size barrier”. When dust
particles aerodynamically coupling with the gas grows into meter size rocks, the dust
particles feel a headwind and drift radially inward because of the gas drag. The infall
time scale is the order of 102 years (Adachi et al., 1976; Weidenschilling, 1977). Thus,
the most solid particles fall into the central star before they grow up into planetesimals.
Moreover, the turbulence in the disk significantly hinders also the growth of meter size
rocks. The disk turbulence excites the relative velocity of dust particles. Although theincrease of collisional frequency leads to the rapid growth under the perfect sticking,but when the collisional speed is higher than several cm/s, the silicate solid particles
bounce or disrupt (Blum & Wurm, 2008; Zsom et al., 2010).
The other classical planet formation scenario using gravitational instability is recognized
as a way to avoid the meter size barrier (Safronov, 1969; Goldreich & Ward, 1973). If
the disk turbulence absent, the dust particles settles on the midplane and forms a dust
layer. When dust density suﬃciently increases, the self-gravitational instability of thedust layer occurs, and eventually planetesimals forms. However, dust layers do not formin the turbulent because the disk turbulence stir up the dust layer (Weidenschilling &Cuzzi, 1993). Thus, from the dust coagulation model, in order to build a successfulplanet formation model, the disk turbulence should be more specified.
Magnetorotational instability (MRI; Balbus & Hawley, 1991), which is an instability
between magnetic fields and ionized gas, is thought to be a most plausible mechanism ofgenerating the disk turbulence. If the MRI grows, the vigorous magnetic turbulence isgenerated with complex magnetic fields (e.g., Hawley et al., 1995). The magnetic fields
lead to Maxwell stress, an accretion stress which is caused by the magnetic tension.
Maxwell stress of vigorous MRI turbulence provides an eﬀective disk viscosity thatallows the disk lifetime consistent with disk observations (Hawley et al., 1995; Fromang& Nelson, 2006; Simon et al., 2009; Flock et al., 2011). In the fully developed MRI
turbulence, the kinetic energy of the turbulence is also enough high to prevent dust
settling (Carballido et al., 2005) and disrupt particles by collisions (Carballido et al.,2010).

7
Since MRI turbulence greatly aﬀects both the disk evolution and planet formation, the
MRI has been studied long time, but the role of MRI in protoplanetary disks is still not
obvious. The awkward problem is that MRI closely links to the ionization fraction of
the disk. Since thermal ionization is relevant only close to the central star (Umebayashi,1983), the dominant part of the disks is ionized only by high-energy sources such as
stellar X-rays (Glassgold et al., 1997) and galactic cosmic rays (Umebayashi & Nakano,1981). Deep inside the disks, the ionization fraction is significantly low because theseionizing radiations are attenuated and because recombination proceeds fast. The low
ionization fraction gives rise to fast Ohmic dissipation that stabilizes the MRI (Sano& Miyama, 1999). Such a region is called thedead zone(Gammie, 1996; Sanoet al., 2000). The MRI is also suppressed by ambipolar diﬀusion near the surface of thedisks (Desch, 2004; Bai & Stone, 2011; Dzyurkevich et al., 2013). The Hall eﬀect caneither stabilize or destabilize the MRI depending on the orientation of the magnetic
field relative to the disk rotation axis (Wardle, 1999; Wardle & Salmeron, 2012; Bai,
2014). Thus, activity of MRI in the disks is closely linked to ionization fraction and its
distribution.
All the previous studies of MRI had neglected the ionization process caused by MRI
itself until Inutsuka & Sano (2005) investigated the possibility of ionization by electricfield in MRI turbulence with a simple estimation. They found heating of electron byelectric fields in MRI turbulence, which is called the “electron heating” and focusedin our work. The heating mechanism is the following. The vigorous MRI turbulence
generates strong electric fields associated with the growth of magnetic fields. Plasma
particles are accelerated by the strong electric fields and are scattered isotropically bycollisions with neutral gas particles, leading to increase of their thermal velocity. In
particular, electrons are more easily heated compared to ions because light particles areeasily scattered. Therefore, the suﬃciently developed electric fields of MRI turbulence
increase electron temperature in a weakly ionized gas. The heated electrons can become
suﬃciently high temperature for collisional ionization, which is known as a phenomenonof electron discharge. If this process works, MRI turbulence sustains (Muranushi et al.,2012).
However, the estimation does not consider balance between ionization and recombi-

8

CHAPTER 1. INTRODUCTION

nation. Okuzumi & Inutsuka (2015, henceforth OI15) investigated ionization balancevarying electric field strength. They found reduction of ionization fraction by heatedelectrons sticking to dust grains before the collisional ionization occurs. The heatedelectrons frequently collide with and stick to dust grains. As a result, the electron
heating decreases the ionization fraction. They also suggested that a region where MRI
is suppressed enlarges in the dusty disk.
Our previous work (Mori & Okuzumi, 2016, henceforce MO16) investigated where theelectron heating takes place and, moreover, might suppress MRI. We found the electron
heating occurs in large regions of protoplanetary disks, which is called “e-heating zone.”
Especially, because the region locates outside the dead zones, electron heating would
eﬀectively enlarge dead zones. In addition, we estimated the strength of magnetic
turbulence in e-heating zones with a simple scaling relation between the Maxwell stressand current density. In e-heating zone, the current density was suggested to be muchless than that in fully developed MRI turbulence. Thus, the previous concluded thatelectron heating would suppress the MRI turbulence. However, the existence of thescaling relation is not clear because many previous studies have not focused the relationof current density to the Maxwell stress. Moreover, the estimation has assumed acurrent density where MRI is saturated. Thus, actually, the possibility of suppressing
MRI by electron heating has been debatable.
In this work, we first investigate the eﬀect of electron heating on the nonlinear evolutionof MRI using local three-dimensional MHD simulations. The goal in this work is toshow the clear relation between the Maxwell stress and current density and, moreover,to demonstrate the potential of electron heating to greatly suppress MRI . We introduce
a resistivity model in which resistivity increases with electric field strength in order tomimic the realistic resistivity of electron heating. If the stress-current relation exists, itappears by varying parameters of current density which begin to cause electron heating
to show the relation.
We denote the following contents with the following order. In Section 2, we present
some fundamental background. The mechanism of MRI and electron heating is included
there. In Section 3, we present our previous study investigating e-heating zone. The
study also describes the implication for dust dynamics. In Section 4, we experimentally9
investigate the eﬀect of electron heating on the magnetic turbulence by performing
MHD simulations including the increase of resistivity. Finally, in Section 5, we present
the summary and conclusion. We also present the future work. In this work, we do notuse consistent notation occasionally. Thus, we define every variables in each chapter
for clarity.

Chapter 2
Background

2.1

Magnetohydrodynamics (MHD) and Magnetorotaional Instability (MRI)

In this section, we briefly denote the fundamentals of magnetohydrodynamics (MHD)
and magnetorotational instability (MRI), with referring to Chen (1977) and Johansen(2009). Magnetorotational instability (MRI) is an instability which takes place in differentially rotating ionized disk. MRI is the most plausible mechanism of driving massaccretion by turbulent viscosity in protoplanetary disks. MRI, which is first discovered in the context of general accretion disks by Velikhov (1959), is rediscovered in thecontext of protoplanetary disks by Balbus & Hawley (1991).

2.1.1

Ideal MHD

MHD deals with a motion of ionized gas holding a magnetic field, by solving the equations of both the fluid and field. In magnetic fields, neutral gas particles receives a force
via collisions with ionized particles which receives the Lorentz force from the fields. Conversely, magnetic fields is aﬀected by the motion of ionized gas as contribution of thecurrent density. An electric conductivity, which is an important physical quantity in
MHD, express the ease with which electric current flows thorough a medium. MHD
based on a limit that the electric conductivity reaches to infinity is called “ideal MHD”.
11

12

CHAPTER 2. BACKGROUND

In this subsection, we describe in the ideal MHD for first understanding of MRI. In thispager, we consistently adopt the Gaussian-cgs units and write magnetic flux density as
magnetic fields.

Fundamental Equations

In protoplanetary disks, contribution of molecular viscosity is much less than of induction and thereby is negligible 1 . It allows us to use the Euler equation, which neglects
the viscosity term in the Navier-Stokes equation, to express the motion of equation of
fluid. In MHD, the fluid receives a force from an electromagnetic field in addition to
gravity and pressure-gradient force. The electromagnetic force in an electric field E
and a magnetic field B per unit volume is generally written as
1
FE.M. = ρe E + J × B,cwhere c is light speed, ρe ≡∑α qα nα is electric charge density, and J(2.3)≡ qα nα uα iscurrent density. Under the first MHD approximation, which holds in an electricallyneutral state , the first term is neglected. Thus, we have a equation of motion of MHD
fluid,∂u
1
GM∗1+ (u · ∇)u = − ∇P − 2 eR + (J × B)
∂t
ρ
r
ρc

(2.4)

where u is the fluid velocity, P is the gas pressure, eR is the radial unit vector of
cylindrical coordinate system.
1

as

A ratio of the advection term to the diﬀusion term of molecular viscosity is approximately written
Lv(v · ∇)v≈= Re,νmol ∆v
νmol

(2.1)

where νmol , v, and L are , respectively, the molecular viscosity, fluid velocity, and scale length. The
value is called Reynolds number Re. To estimate the typical value, we take νmol = cs /nn σnn , L = H,
and v = cs , and then
cs nn σnn
Re == 1012
Ω

(

T
280 K

)−1/2 () ( r )−3/2
nn1015 cm−3
1 AU

(2.2)

where we take σnn = 10−15 cm2 . That’s why the typical value in protoplanetary disks is much larger
than unity and thereby viscosity term can be neglected.

2.1. MAGNETOHYDRODYNAMICS (MHD) AND MAGNETOROTAIONAL INSTABILITY (MRI)1
The equation of continuity is written as∂ρ+ (u · ∇)ρ = −ρ ∇ · u.
∂t

(2.5)

Assuming the isothermal change, we express the equation of state as,
P = c2s ρ.

(2.6)

In the MHD equations, the evolution equation of a magnetic field, the induction equation, should be solved. The Maxwell-Faraday equation is∂B= −c∇ × E,
∂t

(2.7)

which express the time evolution of magnetic fields. The relation between an electricfield and current density, which is called “the Ohm’s law”, in a comoving frame (whichmoves with fluid) is written as
J ′ = σc E ′(2.8)

where we express an electric conductivity as σc and values in the comoving frame withthe superscript “ ’ ” only in this section. Under the Lorentz transformation in u ≪ c,
E ′ and J ′ are related with E and J as
1
E ′ = E + u × B,
c

(2.9)

J′ = J.

(2.10)

Using Equations (2.7), (2.8), (2.9), and (2.10), we obtain∂B= −c∇ ×
∂t

(

)
J
1
− u×B .
σc c

(2.11)

Under the ideal MHD approximation which means σc → ∞, we obtain the induction
equation in ideal MHD,∂B
= ∇ × (u × B).
∂t(2.12)

Maxwell equations except for the Faraday’s equation include Gauss’ equation, equation

14

CHAPTER 2. BACKGROUNDof conservation of magnetic flux, and Maxwell-Ampére equation, respectively,∇ · E = 4πρe ,(2.13)∇ · B = 0,∂E
4πJ += c∇ × B.
∂t(2.14)(2.15)

Under the second MHD approximation, which holds when variation of electric fields
is small, the displacement current term is dropped, and then we obtain the Ampére’sequation,
4πJ = c∇ × B(2.16)

Usually, when the fluid velocity is much less than the speed of light 2 , the approximation
is justified.

The force in Euler’s equation 2.4 is transformed, using Ampére ’s equation (2.16),
into
1
1
J ×B =((∇ × B) × B) ,
c
4π
( 2)
B
1
=
(B · ∇) B − ∇,4π
8π(2.19)(2.20)

where the first term is called the magnetic tension while the second is called the magnetic
pressure. The importance of the magnetic force 3 represents a non dimensional number,plasma beta which is defined by the ratio of gas pressure to magnetic pressure as
β=

2

8πP
.
B2

(2.21)Let us derive the approximation condition from comparing two terms. The condition is approximately written asEB≪c .
(2.17)
T
L
Using approximated Maxwell-Faraday Equation (2.7), B/T = cE/L, the equation can be translated
into
L2≪ c2 .
(2.18)
T2
That’s why when the representative-velocity scale, which is fluid velocity in many cases of MHD, ismuch less than light speed, we can neglect the term.
3
Strictly speaking, the value expresses the unimportance of the magnetic force.

2.1. MAGNETOHYDRODYNAMICS (MHD) AND MAGNETOROTAIONAL INSTABILITY (MRI)1
The transformation leads to the equation of motion of MHD,
∂u
1
GM∗
1
1
+ (u · ∇)u = − ∇P − 2 eR +(B · ∇) B − ∇
∂t
ρ
r
4πρ
ρ

(

B2
8π

)
(2.22)

In summary, the system of equation of MHD consists of Equations (2.4), (2.5), (2.6),(2.12), and (2.16). This system does not explicitly include the electric fields E because
the infinite current streams because of the infinite conductivity and then the fields
disappeares.

Ohm’s Law

Here we denote the electric conductivity σc in Ohm’s law (2.8). In this work, we derive
the Ohm’s law in the comoving frame, but note that it is not necessarily and some books
or papers is written in the rest frame. The reason that we adopt the comoving frame
is to simply derive the nonlinear Ohm’s law which we will denote in Section 2.2.
In the comoving frame, the current density is expressed from the definition,J′ ≡∑qα nα ⟨vα ⟩′ ,(2.23)

αwhere ⟨vα ⟩′ is the mean velocity of α-particles. ⟨vα ⟩′ is derived by the kinetics asfollows. When the charged particles having mass mα and charge qα receive the electricforce in the electric field E’ and then accelerate into ⟨vα ⟩′ during mean collisional time∆tα , using the approximated equation of motion mα ⟨vα ⟩′ /∆tα = qα E ′ ,the velocity isexpressed as⟨vα ⟩′ =qα E ′ ∆tα.
mα(2.24)

Substitution for ⟨vα ⟩′ into Equation (2.23) leads to(J′ =∑ q 2 nα ∆tα
α

α

mα

)
E ′ = σc E ′ .

(2.25)The equation in the parenthesis is the electric conductivity σc . As is evident from theequation, σc is proportional to nα ∆tα (= nα /(nn ⟨σαn vαn ⟩)). Therefore, we can say that

16

CHAPTER 2. BACKGROUND

the ideal MHD approximation means the limit of the ionization fraction nα /nn .
σc is commonly treated as constant, but, strictly speaking, that is just assumption.
Okuzumi & Inutsuka (2015) get rid of the assumption and focus that the ionization
balance can be changed when electric fields are strong. The more detail is described in
Section 2.2.

Alfvén Velocity

Here we denote the Alfvén velocity which is the propagation speed of a magnetic field
displacement. Let us consider the propagation of the displacement B1 in the uniform
magnetic field B0 (= B0 ez ) with zero initial velocity (u0 = 0). Moreover, for simplicity,
we assume constant density, which leads to no pressure gradient(∇P = 0) and incompressible fluid(∇ · u = 0), and B1 to be perpendicular to B0 , B1 = B1x ex + B1y ey . The
equation of motion can be written as
∂u1
1
= −
(B0 × (∇ × B1 )) ,
∂t
4πρ
1
=
(B0 ∇z B1 ) ,4πρ

(2.26)

which have been derived from Equation (2.19). The induction equation is writtenas∂B1= ∇ × (u1 × B0 ),∂t= (B0 · ∇)u1 − B0 (∇ · u1 ),= B0 ∇z u1 .

(2.27)We then have the wave equation of B1 toward the z-direction,()∂ 2 B1
1= B0 ∇z
B0 ∇z B1 ,∂t2
4πρB02 ∂ 2=B1 .
4πρ ∂z 2(2.28)

2.1. MAGNETOHYDRODYNAMICS (MHD) AND MAGNETOROTAIONAL INSTABILITY (MRI)1
This equation expresses that the magnetic displacement is propagated with a velocity
B0vA = √,
4πρ

(2.29)

which is called “Alfvén velocity.”

2.1.2

Linear Stability Analysis of MRI

MRI appears by the linear stability analysis of ideal MHD in diﬀerentially rotating
disks in magnetic fields.
Let us consider the Kepler disk in the presence of an initially uniform magnetic field inz direction (B0 = (0, 0, B0 )). The angular velocity Ω at a distance r from the central
star is expressed as√Ω(r) =

GM∗.
r3(2.30)In the cylindrical coordinate (r, ϕ, z), the equation of motion is written as∂ur u2ϕ
1 ∂P
GM∗
1
−
+ (u · ∇)ur = −− 2 er + (J × B)r ,
∂t
r
ρ ∂r
r
ρc
∂uϕ ur uϕ1 ∂P1++ (u · ∇)uϕ = −+ (J × B)ϕ ,
∂t
r
ρr ∂ϕ
ρc
∂uz
1 ∂P
GM∗
1
+ (u · ∇)uz = −− 2 ez + (J × B)z .
∂t
ρ ∂z
r
ρc

(2.31a)(2.31b)(2.31c)Local Shearing Sheet Approximation

Some proper approximations help to simply derive the linear stability analysis of MRI.
Here let us adopt the local coordinate (x, y, z), where x-axis is toward r-axis and yis perpendicular to x. The coordinate rotates with the angular velocity Ω0 = Ω(r0 ),where r0 is the center of the coordinate and much larger than x and y, r0 ≪ x andr0 ≪ y. x and y are, respectively, expressed as(x, y) = (r − r0 , r0 (ϕ − Ω0 t)),(2.32)

18

CHAPTER 2. BACKGROUND

Local Sharing Boxy!

Central Star

x!
r0!
Ω0!

Figure 2.1: Schematic diagram expressing the local sharing coordinate. The origin
is located at a distance r0 from the central star and rotates with an angular velocityΩ(r0 ). In diﬀerentially rotating disks, an observer in the coordinate observes a shear
flow.

and the velocity in the local coordinate, (ux , uy ), is related to the one in the cylindricalcoordinate, (ur , uϕ ), as(ux , uy ) = (ur , r0 (Ω(r) − Ω0 )) .

(2.33)

where Ω(r) = ∂ϕ
. To know the background velocity field (unperturbed velocity) in the∂tlocal coordinate, we consider the steady flow and then have the velocity field,()
3(ux,steady , uy,steady ) = 0, − Ω0 x ,2(2.34)

where we use
√

(
)−3/2
x
GM∗
− Ω(r0 )1+Ω(r0 + x) − Ω(r0 ) =
r03
r0
(
)
3x
≈ Ω0 1 −− Ω(r0 )2 r0
3 x≈ − Ω0 .
2 r0(2.35)

2.1. MAGNETOHYDRODYNAMICS (MHD) AND MAGNETOROTAIONAL INSTABILITY (MRI)1
The tidal force (a centrifugal force − −a gravitational force) in Equation (2.31a) can be
written as
u2ϕ GM∗− 2 = (r0 + x)
r
r

(

∂ϕ
∂t

)2

GM∗
−
(r0 + x)2
(
)2
(
)−2
uy
x
GM∗
1+
= (r0 + x)
+ Ω0 − 2
r0
r0
r
(
)
( 0
)
x
uy Ω0
x
2
2
≈ r0 (1 + ) Ω0 + 2− r0 Ω0 1 − 2
r0
r0
r0≈ 3Ω20 x + 2uy Ω0 ,

while

∂uϕ
∂t

(2.36)

can be written as
( (
))
∂uϕ
∂
uy
=
r
+ Ω0
∂t
∂t
r0
ux uy
r ∂uy
=
+ ux Ω0 +
r0
r0 ∂t
∂uy
≈ ux Ω0 +.
∂t(2.37)

ur uϕ /r can be written as
ur uϕ
= ux
r

(

uy
+ Ω0
r0

)

≈ ux Ω0 .

(2.38)For the simplicity of the derivation, we here use the shearing sheet approximation which
neglects the vertical velocity and the equation of motion with respect to z direction.
In addition, we impose the incompressibility, ρ = ρ0 (constant), and isothermality, cs =constant, and then also have P = constant. Thus, we obtain the equations to be solved,∂ux
1+ (u · ∇)ux = 2Ω0 uy + 3Ω20 x+ (J × B)x ,
∂t
ρc
1
∂uy
+ (u · ∇)uy = −2Ω0 ux+ (J × B)y ,∂tρc∂B+ (u · ∇)B = (B · ∇)u,∂t(2.39a)(2.39b)(2.39c)20

CHAPTER 2. BACKGROUND

where we use ∇ · u = 0.

To investigate the growth of given perturbations, we perform a linear stability analysis.
We express values as a sum of the unperturbed value and perturbed value, and thentransform the above equations into the ones of the first-order vales. We consider thatthe unperturbed magnetic field is vertically uniform and the perturbed magnetic field
is only perpendicular to z direction. The values are expressed asux = ux1 ,(2.40a)3
uy = − Ω0 x + uy1 ,2(2.40b)Bx = Bx1 ,(2.40c)By = By1 ,(2.40d)Bz = B0 ,(2.40e)J=

1∇ × B1 = J 1 .
4πc(2.40f)Using these equations and neglecting the second-order values, we obtain∂ux1
1= 2Ω0 uy1 +(J1 × B0 )x∂t
ρ0 c∂uy1 31− Ω0 ux1 = −2Ω0 ux1 +(J1 × B0 )y
∂t
2
ρ0 c
∂Bx1
∂ux1
= B0
∂t
∂z
(
)
∂By1
∂uy1
3
= B0
+ Bx1 − Ω0 .
∂t∂z
2

(2.41a)(2.41b)(2.41c)(2.41d)The Lorentz force can be written as))
1 (( c1(J ′ × B0 ) =∇ × B1 × B0
ρ0 c
ρ0 c
4π
1
=
((B0 · ∇) B1 − ∇ (B0 · B1 ))
4πρ0
1
∂B1
=
B0
.
4πρ0
∂z

(2.42)

2.1. MAGNETOHYDRODYNAMICS (MHD) AND MAGNETOROTAIONAL INSTABILITY (MRI)2
Thus, we obtain∂ux1
B0 ∂Bx1= 2Ω0 uy1 +∂t
4πρ0 ∂z1
B0 ∂By1∂uy1= − Ω0 ux1 +∂t
2
4πρ0 ∂z
∂Bx1
∂ux1
= B0
∂t
∂z
∂By1
∂uy1 3= B0− Ω0 Bx1 .
∂t∂z
2(2.43a)(2.43b)(2.43c)(2.43d)Let us investigate an unstable mode in z direction, kz , by the Fourier analysis. The any
given perturbation can be expressed by a superposition of plane waves. We assume an
arbitrary perturbed value q1 to beq1 ∝ exp [i(kz z − ωt)].

(2.44)

Thus, we obtain
B0
Bx1
4πρ0
1
B0
−iωuy1 = − Ω0 ux1 + ikz
By1
2
4πρ0

−iωux1 = 2Ω0 uy1 + ikz(2.45a)(2.45b)−iωBx1 = ikz B0 ux1(2.45c)3−iωBy1 = ikz B0 uy1 − Ω0 Bx1 .
2(2.45d)Thus, by arranging this equations, we obtain the dispersion relation with respect tokz ,ω 4 − 2ω 2 (2vA2 kz2 + Ω20 ) + vA2 kz2 (vA2 kz2 − 3Ω20 ) = 0.

(2.46)

Equation (2.46) can be easily solved. Finally, we obtain the solution of the dispersion
relation,2ω±
vA2 kz2 1
=
+ ±
Ω20
Ω20
2

√

v2 k2 1
4 A 2z +Ω04(2.47)In Figure 2.2, we plot ω 2 as a function of kz . Since ω is given in Equation (2.44),
the perturbation is unstable when ω 2 (kz ) < 0 , while the perturbation is stable when

22

CHAPTER 2. BACKGROUND
1.5

0.5

2

ω /Ω0

2

1

0

-0.5
-1
0

0.5

1
1.5
kzvA/Ω02

2.5

Figure 2.2: Dispersion relation of MRI in ideal MHD. The frequency in the verticalaxis is normalized by Ω0 , and the perturbation wavenumber in the horizontal axis is22
normalized by Ω0 /vA . Blue line shows ω−/Ω20 , and red line ω+/Ω20 . When ω 2 < 0, the2perturbation is unstable, while when ω > 0, the perturbation is stable and propagate
as a MHD wave.
ω 2 (kz ) > 0. Thus, the unstable wave number is
0 < kz <
where

√

√
3Ω0 /vA ,(2.48)

3Ω0 /vA is the largest wavenumber in the wavenumber which trigger MRI. The

shortest unstable wavelength λcrit which corresponds to the largest unstable wavenumber isvA
λ = 2π √.
3Ω0(2.49)The unstable wavelength in which MRI growth is fastest is the one of min(ω 2 ), andthereby we obtain the most unstable growth rate,√
λ = 2π

16 vA
15 Ω0

(2.50)

when the wavelength is the most unstable wavelength,
3
Im(ω)max = Ω0 .
4

(2.51)The timescale of MRI growth is ∼ 1/Ω0 . This means that the instability grows during an
orbital period, and the growth rate is relatively fast in global eﬀects of protoplanetary2.1. MAGNETOHYDRODYNAMICS (MHD) AND MAGNETOROTAIONAL INSTABILITY (MRI)2
disks. Moreover, in ideal MHD, the growth rate does note depend on the verticalmagnetic field strength. In other words, MRI would appear even if the initial magnetic
fields is very weak. That’s why it has been thought that MRI should play an importantrole in protoplanetary disks.

2.1.3

Nonideal MHD Eﬀects

Ideal MHD approximation assumes infinite conductivity, whereas “nonideal MHD” havefinite conductivity.
In nonideal MHD, using Equations (2.16) and (2.14), the induction equation is expressedas∂Bc2= ∇ × (u × B) −∇ × (∇ × B)
∂t
4πσc
c2
= ∇ × (u × B) −(∇(∇ · B) − ∇2 B)
4πσc= ∇ × (u × B) + η∇2 B,(2.52)

where η is magnetic resistivity which is related to electric conductivity as
η=

c2
.
4πσc

(2.53)The right-hand side first term in Equation (2.52) , which is the same as in ideal MHD,
expresses the amplification of magnetic fields by gas motion coupling with the fields.
The second term expresses the diﬀusion of magnetic fields. When the first term is
negligible, the equation is expressed as a diﬀusion equation. Therefore, when η iseﬀective, the magnetic fields grow slowly or decrease. Thus the dissipation of magnetic
field by high resistivity or low conductivity is called the ohmic dissipation.
We can obtain a criterion which expresses whether the magnetic field is induced or
diﬀused, by taking a ratio of the induction term in Equation (2.52) to the diﬀusion
term. If we consider partial derivatives to be the most unstable wavenumber Ω/vA andcharacteristic velocity to be vA , we obtain the criterion, which is called the Elsasser

24

CHAPTER 2. BACKGROUND

number Λ,Λ=vA2|∇ × (u × B)|
∼
.
η∇2 B
ηΩ

(2.54)

When Λ ≫ 1, the ohmic dissipation is negligible and thereby we can assume ideal
MHD. In contrast, when Λ ≲ 1, the ohmic dissipation is not negligible and thereby weshould consider nonideal MHD.

2.2

Nonlinear Ohm’s Law

In this subsection, we summarize the nonlinear Ohm’s law in OI15. The nonlinear
Ohm’s law is a relation between electric field and current density including heatingof electron by strong electric field. So far, previous studies assumes that electric conductivity (electric resistivity) does not depend on the electric field strength. However,actually, MRI turbulence generates strong electric field after besides magnetic field
growth. Thus, the possibility that the electric field heat up electrons was suggested inInutsuka & Sano (2005). After that, in weakly ionized gas including dust, Okuzumi &Inutsuka (2015) found that an increase of electron energy leads to increase of electron’s
dust adsorption rate. Since the reduction rate of electron abundance is increased, theelectron abundance is decreased in a gas-phase. Thus, the conductivity depend on theelectric field because electric conductivity depends on the ionization fraction. For MRI
to grow, the ionization fraction should be suﬃciently high (see Section 2.1.3). Thus,
because MRI might be suppressed if this eﬀect suﬃciently work, this eﬀect will be
important.

2.2.1

Outline of Calculation

Here, we present the outline of derivation of the nonlinear Ohm’s law. The keys ofthe nonlinear Ohm’s law are heating of electron by electric field and the the ionization
balance. In electric fields, electron velocity distribution is changed from the thermalvelocity distribution, Maxwell distribution. The change of the velocity distribution
leads to change mean drift velocity ⟨vα ⟩ and mean energy ⟨ϵα ⟩ of plasma particles
in turn. Because the mean plasma energy is closely related to ionization fraction,

2.2. NONLINEAR OHM’S LAW

25

Figure 2.3: Sketch of the derivation procedure and the charge reaction model presented
in OI15. The model gives the electric current density J as a function of the electricfield strength E. The plasma velocity distribution functions fe and fi take into accountheating by the electric field, and therefore the resulting J is nonlinear in E.
ionization balance is aﬀected and therefore number density of plasma particle nα is∑changed. Current density is defined as J = α=e,i qα nα ⟨vα ⟩, where q is charge and α
express a value of plasma particle. Thus, current density including dependence of the
conductivity on electric fields can be calculated, and then the relation between electricfields and current density is obtained.
The order of the derivation is as follows:
1. Calculate electron and ion velocity distribution function fe and fi .
2. From the velocity distribution, calculate electron and ion mean velocity, ⟨ve ⟩ and⟨vi ⟩, and mean energy ,⟨ϵe ⟩ and ⟨ϵi ⟩.
3. Calculate electron and ion number density, ne and ni , from ionization balance.
4. Obtain the J–E relation from J =∑α qα nα ⟨vα ⟩.

We show the schematic diagram of the derivation procedure in Figure 2.3. In theirpaper, they focus also a collisional ionization of electrons by strong electric field besides
dust adsorption. However, because our work focus the nonlinearity by dust adsorption ofheated electron, we does not treat the collisional ionization here. In addition, because
the collisional ionization occurs in E ∼ 103 Ecrit , the neglecting it is reasonable inE ≲ 102 Ecrit . Actually, as we see Section 3.2 of our work, we confirm the electric fields
does not increase so much in a large region of the disk.

26

2.2.2

CHAPTER 2. BACKGROUND

Plasma Heating

Here we briefly present a derivation of heating of plasma particles by electric field(plasma heating) in OI15.
A mean velocity ⟨vα ⟩ and energy ⟨ϵα ⟩ of charged particles having a velocity distributionfunction fα is expressed as∫⟨vα ⟩ =⟨ϵα ⟩ =∫

vα fα (E, vα ) d3 vα ,(2.55)

ϵα fα (E, vα ) d3 vα .

(2.56)

where ϵα = mα vα2 /2. Using ⟨ϵα ⟩, the temperature of charged particles is defined as⟨ϵα ⟩ ≡ 3kB Tα /2 when the random (thermal) energy dominates the kinetic energy.
The kinetic energy of electrons exceeds the one of neutral particles when electric fieldstrength exceeds the threshold Ecrit , (Landau & Lifshitz, 1953),√
Ecrit =

6me kB T,
mn ele(2.57)

where me and mi are, respectively, the mass of electron and neutral particles, kB is theBoltzmann constant, T is the neutral gas temperature, e is the elementary charge, andle is the mean free path of electron.
The velocity distribution function of electrons depending on electric field strength isderived by Davydov (1935). OI15 used the distribution although here the specificformula is skipped for simplicity. Figure 2.4 shows the E dependent energy distributionof electrons in the case of E = 0 and 100Ecrit . When E ≈ 0, the energy distribution
is approximately the Maxwell distribution. From E = 0 to 100Ecrit , the top of theenergy distribution which corresponds to ⟨ϵe ⟩ shifted to the right. Thus, electric field
increases the mean kinetic energy. The velocity distribution function of ions fi areapproximately expressed as a Maxwell distribution that is shifted by ion mean drift
energy mi ⟨vi ⟩2 /2 (Hershey, 1939). The drifted velocity can be derived from the kinetics(e.g., see Appendix in OI15).
Figure 2.5 shows the dependence of electron and ion mean energy on electric fieldexpressions ca
(17)
27

101
E$0

0

10

$

$

E $ 100Ecrit

10!1
10!2
10!3
10!4
10!5 !1
10

100

101
102
Electron Energy Εe !kB T

103

101
Mean Energies "Εe #, "Εi # $eV%

Energy Distribution 4Πv3e fe0 "E, ve #temperature equal to that of the neutrals,
!
"3/2
!
"
me
ϵe
(M)
fe0 =
exp −
.
2πk
T
k
T
B
B
2.2. NONLINEAR OHM’S LAW

100

10"1

10"2

10"3
10"2

Figure 3. Mean kand ⟨ϵi ⟩ (Equation
Figure 2. Electron energy distribution as a function of the electron energy ϵe
The neutral gas tem
for
E
=
0
(dashedcurve)
andE=
100E(solidcurve).
The
ticks
indicate
Figure 2.4: Energy distribution of electronscritdepending on electric field strength invertical lines mark
ϵ
=
⟨ϵ
⟩,
where
⟨ϵ
⟩
is
the
mean
electronenergy(Equation(20)).
ethe casee of Ee = 0 (dashedline) and 100Ecrit (solid line). The vertical ticks on the
distributions indicate ⟨ϵe ⟩.

strength. When E ≲ Ecrit , the mean kinetic energy is constant and equal to the neutral
gas temperature. When E > Ecrit , the electron mean kinetic energy increases withincrease of E. Likewise, when E is larger than Ecrit,i , an electric field strength whichstarts to heat ions, the ion kinetic energy increases with increase of E. Thus, strong
electric field heat up plasma particles.
The critical Ecrit and Ecrit,i are diﬀerent by ∼ 300. Actually, mean free time of electrons
and ions are the same order, because polarization of ion increases the cross section
although electrons move faster than ions.

2.2.3

Charge Reactions

Ionization fraction can be obtained by solving the rate equation including ionization
and reduction reactions. OI15 consider ionization, recombination, and dust adsorption.
In addition, OI15 adopted an analytical method for solving the rate equation presented
in Okuzumi et al. (2009). They assume the equilibrium state and the ions represented

ution with thetotic forms in the weak and strong field limits. The asymptotic(M)
expressions can be derived by direct substitution of fe0
and

28
(17)

.

CHAPTER 2. BACKGROUND

T ! 100 K
100

"Εe #

10"1

10"2

10"3
10"2

3

10

E!Ecrit, i"Εi #E!Ecrit$ 100Ecrit

Mean Energies "Εe #, "Εi # $eV%

101

10"1
100
101
102
Electric Field Strength E!Ecrit103

Figure 3. Mean kinetic energies of electrons and ions, ⟨ϵe ⟩ (Equation (20))
2.5:of mean plasma energy on electric field strength. They tookand ⟨ϵDependence
electron energy ϵFigurei ⟩ (Equation (26)), as a function of the normalized field strength E/E crit .
e
here
gas
temperature
tobe
100 K.TVertical
axistoshows
kineticof electrons
The neutral gas temperature
is assumed
be 100mean
K. Thesolid energy
and dashed
The ticks indicate
verticalmark
E=
Ecrit and Eaxis= Eshows(27)),field
respectively.
n (20)).
(blue) andionslines(red),
andhorizontal
the electricstrength notarized
crit,i (Equation
Ecrit . Ecrit and Ecrit,i shows the electric fields where electron and ion heating starts up,
respectively.

by a single dominant species. The rate equations to be solved are
0 = ζnn − Krec ni ne − Kdi nd ni ,(2.58a)0 = ζnn − Krec ni ne − Kde nd ne ,(2.58b)0 = Znd + ni − ne ,(2.58c)where nx (x=n(neutral particles), e(electrons), i(ions), d(dust particles)) is number
density of x particle species, ζ is the ionization rate, Kdα is the dust adsorption rate ofcharged particles, Krec is the recombination rate, and Z is the number of charges per
a dust grain. Equation (2.58a) and (2.58b) can transform into a quadratic equation.
(eq)Thus, the analytical solution nα

is given by

ζnn
=
n(eq)
α
Kdα ndZ (eq) is given by
Z=

(

√
1
+
2

1
Krec ζnn+
4 Kdi Kde n2d)
1 ( (eq)(eq).
ne − nind)−1.

(2.59)(2.60)

2.2. NONLINEAR OHM’S LAW

29

Dust adsorption rate Kdα is the ensemble averaged volume that passes through collisional cross section in the mean free time. Kdα is given by∫
Kdα = Sασdα (vdα )vdα fα (vα )d3 vα(2.61)

where the cross section between a dust grain and a charged particle, σdα (vdα , Z), is
given by Rutherford scattering (Spitzer, 1941) and Sα is the absorption probability.
OI15 assumes Sα = 1.

2.2.4

ExamplesFigure 2.6 shows the ionization fraction depending on an electric field. In low electricfields, electrons and ions are in equilibrium. This means that gas-phase recombination
process controls ionization fraction. While electric field strength exceeds Ecrit , theelectron abundance decreases with increase in electric field. At this time, ions andcharges on dust grains are in equilibrium. The ion abundance is controlled by dust
adsorption rate, whereas the electron abundance is controlled by collisional frequency
to dust grains, i.e. electron temperature. While the electron fast sticks to dust grains,electron abundance is diﬀerent from ion abundance.
The nonlinear Ohm’s law can be derived by the definition of current density, using⟨vα ⟩ (E) and nα (E) . In low electric fields, J is linear in E, but the linear relation isbroken around E = Ecrit . Then, at E > Ecrit , current density decreases because thedust adsorption rate of electron, Kde , steeply increases.
As stated above, OI15 showed the nonlinear relation between current density and electric field strength by considering electron heating. The nonlinearity comes from thedependence of electric resistivity on electric field strength. The resistivity increased byelectron heating might suppresses MRI. This eﬀect has never considered in protoplanetary disk so far. The suppression of MRI impacts on both the disk evolution and dust
growth. In the next chapter, we investigate zones where the electron heating stabilizes
MRI.

10!10

10!18 !11
10!7
10!6
10
10!10
10!9
Electric Field Strength E CHAPTER!esu cm!2 " 2. BACKGROUNDElectric Field Stren

10!10

Model C!w

10!10

Model D

10!12

Nonlinear Ohm’s Law!Zxd
xe
Model A

E"Ecrit

!2
10!16
10

ption#"o adsor

xi

10!12

10!14
10!1

10!8

E"Ecrit, iAbundances
xe , xi , !Zxd
!2 !1
Current Density J !esu cm
s "

30

10!9

!2 !1
Current Density J !esu cm
s "Abundances
xe , xi , !Zxd

10!18 !11
10

10!14
10!1
!2
10
10!16

Model B

10!3
10!3
10!18 !11
10!18 !11
!10
!9
!8
!7
!6
10
10
10
10
10
10
10
10!10
10!9
10!4
10!4
Electric Field Strength E !esu cm!2 "Electric Field Stren
10x!5
10!5
Figure 6. Ion abundance
grain charge abundance −Zi (dotted curve), electron abundance xe (dashed curve), andfield strength
for modes A,of
B,plasmaC, and D.
The dot-dashedcurve shows
gas-phasechargeabundance in the case
Figure2.6: EDependence
abundanceon electricfield thestrength.
The
verticalanddashedvertical
lines
mark
E
=
E
and
E
=
E
,
respectively.
!6
!6
crit
crit,iaxis shows the abundances,and the horizontal axis shows electric field strength.
The
10
10
red dashed line shows the ion abundance xi = ni /nn , the blue dashed line the electron
the Coulomb
thelinegrain
electrons
regu- of charges
resultis dust!7 between
abundance
xe =repulsion
ne /nn10, thegreen solid
the andabundanceof number
10!7 on
!11
!10
!9
!8
!7
!6
lates
the
electron-ion
flux
ratio
measured
on
the
grain
surface
10
10
10
10
10
10
10!10
10!9n ≈1
grains −Zxd = Znd /nn . Ecrit and Ecrit,i shows the electric fields where electron and ioneto be unity.
heating
start, respectively. Electric Field Strength E !esu cm!2 "Electric Field Strengt

Current Density J !esu cm!2 s !1 "E"EcritCurrent Density J !esu cm!2 s !1 "We look at how the ionization balance changes by adding
In this high field limit (Edust grains. Model10B!2has the dust-to-gas mass ratio fdg that10!4Equation (38))(seeis 100 times higher than model A. We see that the gas-phase
C to decrease at high field
Model D
charge abundances 10xi!3andModelxe start
10!5
ne ≈
−8
−2
strengths (E ! 10 esu cm ≈ 10Ecrit ). This is an indication that plasma adsorption
takes over gas-phase recombina!610!4
10
The enhanced electrontion at high E. As mentioned above, the heating of electrons
the !7charge balance in gas
suppresses gas-phase !5recombination since Krec is a decreasing
10
10
creases with E while xe dfunction of ⟨ϵe ⟩. In contrast, electron heating accelerates√thetend to dominate over fre
adsorption of electrons
onto grains since Kde ∝ ⟨ve ⟩ ∝ ϵe .
!8
Je
10!6
In10this limit, the charge
As a result of the two eﬀects, the adsorption generally overlished between ions and c
Ji strengths. As
whelm gas-phase recombination
at high field!7
10!9electrons being a min10
freelong as ne ≫ |Z|nd , the ion number density also decreases
model B, this is only appr
with E to maintain the
charge neutrality ni ≈ ne . Obvi!10!8
10 !11takes !10
E10!
10−8 esu cm−2 ≈ 10
ously, plasma adsorption
over gas-phase
recombination
10!11 10!10 10!9 10!8
10
10
10!9
10!8
10!7
10!6
occurs at much lower E (at lower E when nd (or fdg ) is higher. This explains why
we!2
Electric
Strengt
Electric
Field
Strength
E
!esu
cm
"
A decreasein ζField
gives
a si
observe the decrease in xe and xi in model B but not in model2
becauseA ∝ A,nd /ζ.
modFigure 10. Current density J (solid curve) as a function of the electric field strength
E for models
B, C,Inand
D.
A.
contributions
fromelectrons
and
ions,
J
and
J
,
respectively.
FormodelsA
andB,
thecurves
for
J
and
J
are
in
100,
|Z|x
dominates
over
e
i
e
d
Figure2.7:
Nonlinear
Ohm’s
law
of
OI15.
The
blue
dashed
line
shows
the
electron
For later convenience, we here derive the approximate exE = Ecrit .
current
density,dashednumber
line the ion
current
solid line the sumof thesepressionof theredelectrondensityne density,
for the blackcase where
7 Okuzumi (2009) called the"√current
density.
The0.35
model
assumes
thatgas-phasethe dust-to-gas
ratio fdg = 0.01, the ionizationelectronadsorption
dominates
overrecombination.
10
electronandion–dust
plasmaproportional
to E
becausene ∝ 1/ TKrec= (Equation(52)),density of neutral gas ≈ 5 × 10−5 sta−17 −1rate= 10bes , neutral
temperature
100
K, 1,
numberThisζ canderivedfromgasEquation(47) withR≫
and thencase and xi ≈ |Z|xd in the latter fc−0.69
Krec ∝ n⟨ϵe =⟩ 1012(Equation⟨ϵe ⟩ grains∝ E (Equation
particlecm−3 , the (34)),
radius andof dusta = 1µm,(22)).
and the material
physicsdensity(e.g., Shukla" & Mamunn#0.85−3J ≈ Je is proportional to E
Therefore,atE≫
Ecrit , whichby the dimensionless Havnes
ρ• = 2g cm .
ζ par×is not significantly diﬀerent from the linear Ohm’s law J ∝ E
10−17 sseen in E ≪ Ecrit .
In model B, the behavior of J is the same as that in model A
Hence, Je is independent ofas long as gas-phase recombination dominates. A new, intertion factor Fde . In model B
esting behavior appears when plasma adsorption onto grains
E because the grain surface
takes over (E ! 10−8 esu cm−2 ). Again, this behavior canChapter 3

Electron Heating Zone : Mori &Okuzumi (2016)

In this chapter, we present our previous study, MO16. Inutsuka & Sano (2005) foundthe electron heating as a new ionization source in magnetic turbulence, i.e. collisional
ionization by energetic electrons. After that, OI15 has found the electron heating
would change ionization balance by electrons sticking to dust before the collisional
ionization. However, these studies did not clarify the locations and the impacts on
disk evolution. We here investigate the region where electron heating would aﬀect MRI
turbulence.

In Section 3.1, we present the disk model, simplified plasma heating model, and ionization balance. In Section 3.2, we present some conditions for MRI growth and somecriteria for mapping of turbulent state in a disk. We also briefly summarize the turbulent state and calculation steps. In Section 3.3, we show where the electron heatingaﬀects MRI turbulence. We also consider cases with various parameters. In Section 3.4,
we estimate how the electron heating suppresses MRI turbulence. In Section 3.5, we
discuss the eﬀect of heated electrons on the electric repulsion and the collisional growthof dust grains.
31

32

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

3.1

Disk and Ionization Models

3.1.1

Disk Model

We consider a gas disk around a solar-mass star. We assume that the surface densityof the disk gas obeys a power law
Σ(r) = 1.7 × 103 fΣ

( r )−3/2
g cm−2 ,
1 AU

(3.1)

where r is the distance from the central star, and fΣ is a dimensionless parameter. The
choice of fΣ = 1 corresponds to the minimum-mass solar nebula (MMSN) model ofHayashi (1981), which we take as the fiducial model.
We assume that the disk is optically thin and give the temperature profile as (Hayashi,
1981)
T (r) = 280( r )−1/2
K,
1 AU

(3.2)

where the central star is assumed to have the solar luminosity.
The sound speed is given by cs =√kT /mn , where mn is the mass of a neutral gas par-ticle, and k is the Boltzmann constant. Assuming mn = 2.34 amu and using Equation(3.2), we have( r )−1/4
cs (r) = 1.0 × 10
cm s−1 .
1 AU
5

(3.3)

We assume that the gas disk is hydrostatic in the vertical direction and give the vertical
distribution of the gas density as()
z2
ρ(r, z) = ρc (r) exp − 2 ,
2H

(3.4)

where ρc is the mid-plane density and H ≡ cs /Ω is the gas scale height with Ω =2.0 × 10−7 (r/1 AU)−3/2 s−1 being the orbital frequency (note that a solar-mass star is√∫∞assumed). Using the relation Σ = −∞ ρdz = 2πHρc , we haveρc (r) = 1.4 × 10−9 fΣ( r )−11/4
g cm−3 .
1 AU

(3.5)

3.1. DISK AND IONIZATION MODELS

33

Thus, the number density of gas particles nn = ρ/mn is given as()( r )−11/4
z2
nn (r, z) = 3.5 × 10 fΣ ×
exp − 2 cm−3 .
1 AU
2H
14

(3.6)

As we will describe in Section 3.2.1, the criteria for MRI depends on the magnetic
field strength in the disk. Following Sano et al. (2000), we consider a net (large-scale)vertical field Bz0 threading the disk and specify its strength with the plasma beta at
2
the midplane, βc ≡ 8πρc c2s /Bz0
. If we use Equations (3.3) and (3.5), the net verticalfield strength can be expressed as(
1/2
Bz0 (r) = 0.59fΣ

βc
1000

)−1/2 (

r )−13/8
G.
1 AU

(3.7)

For simplicity, we will assume that βc is constant in the radial direction.
The charge reaction model adopted in this study takes into account the eﬀects of graincharging on the ionization balance. For simplicity, we assume that dust grains are well
mixed in the gas so that the dust-to-gas mass ratio fdg is a global constant. We also
assume that the grains are spherical and single-sized with radius a (taken as a freeparameter) and internal density ρ• (fixed to be 3 g cm−3 ). From these assumptions,the number density of dust grains nd is given by 3fdg ρ/(4πa3 ρ• ), which is expressed
as
)−1 (
)−3
ρ•
a
nd (r, z) = 1.1 × 10 fΣ
3 g cm−3
0.1 µm()( r )−11/4
z2
×exp − 2 cm−3 .
1 AU
2H
(

3

fdg
0.01

)(

(3.8)

The disk is assumed to be ionized by galactic cosmic rays, stellar X-rays, and radionuclides. The ionization rate can be expressed asζ = ζCR + ζXR + ζRN ,(3.9)

where ζCR , ζXR , and ζRN stand for the contributions from cosmic rays, X-rays, andradioactive decay, respectively. The cosmic ray distribution is expressed as (Umebayashi

34

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)& Nakano, 2009)

ζCR


(
)[
(
)3/4 ]−4/3χχζCR,0 
exp −
=
1+
2 
χCR
χCR

(
)[
(
)3/4 ]−4/3 Σ−χΣ−χ
+ exp −
1+
,

χCR
χCR

(3.10)

where ζCR,0 = 1.0×10−17 s−1 is the characteristic ionization rate of cosmic rays, χ(r, z) =∫∞ρ(r, z ′ )dz ′ is the vertical gas column density above height z, and χCR = 96 g cm−2 iszthe attenuation depth of ionizing cosmic rays. The ionization rate of X-rays is expressedas (Bai & Goodman, 2009)ζXR =( r )−2.2
LX
1029 erg s−1 1 AU
{
[
( (
( (
)0.4 ))0.4 )]χΣ−χ
× ζXR,1 exp −+ exp −
χXR,1
χXR,1
( (
[
( (
)0.65 )
)0.65 )]}Σ−χχ
+ exp −,+ζXR,2 exp −
χXR,2
χXR,2

(3.11)

where χXR,1 and χXR,2 are taken to be 6 × 10−3 g cm−2 and 3 g cm−2 respectively,
ζXR,1 and ζXR,2 are taken to be 6 × 10−12 s−1 and 1 × 10−15 s−1 respectively. We takeLx = 2 × 1030 erg s−1 in accordance with the median X-ray luminosity of solar-mass
young stars (Wolk et al., 2005). The ionization rate of the radionuclide is expressed as(Umebayashi & Nakano, 2009)(
ζRN = 7.6 × 10

3.1.2

−19

fdg
0.01

)
s−1 .

(3.12)

Simplified Plasma Heating Model

As we will describe in Section 3.2.1, the criterion for MRI depends on the ionization
fraction in the disk. We employ a simple ionization model proposed by OI15 to calculate
the ionization fraction taking into account plasma heating by a strong electric field.
The model determines the ionization fraction of the gas at each location of a disk
from the balance between ionization by external high-energy sources (e.g., cosmic rays

3.1. DISK AND IONIZATION MODELS

35

and X-rays), recombination in the gas phase, and adsorption of ionized gas particles
onto dust grains. The rates of recombination and adsorption generally depend on thetemperatures of ions and electrons, Ti and Te . Previous ionization models assumed
that Ti and Te are equal to the neutral gas temperature T . By contrast, the model of
OI15 determines Ti and Te as a function of the electric field strength E. For simplicity,
positive ions are represented by the single species HCO+ , which is good as a first-orderapproximation when heavy molecular ions that recombine through dissociation reactions
dominate (Umebayashi & Nakano, 1990; Dzyurkevich et al., 2013). We do not consider
negative ions. Although production of negative ions is rare in cool protoplanetary disks,electrons heated to ≳ 3 eV can produce negative hydrogen ions H− via dissociativeelectron attachment H2 + e− → H− + H (Wadehra, 1984). However, H− would beinstantly destroyed by CO, the most abundant molecule after H2 , via the reaction
H− + CO → HCO + e− (Ferguson, 1973). For this reason, we may safely neglect the
dissociative electron attachment during electron heating.
In this study, we make two further simplifications to the original model of OI15. Firstly,
we calculate the electron temperature Te by solving the equations of momentum andenergy conservation rather than by using the solution to the full Boltzmann equation.
The rate coeﬃcients for gas-phase recombination and plasma adsorption onto grains arethen evaluated by approximating the velocity distribution function with a Maxwellian
with temperature Te . The approach greatly simplifies the analytic expressions of therate coeﬃcients that otherwise involve confluent hypergeometric functions (see Section
3 of OI15). Such an approach was originally proposed by Hershey (1939) for calculating
the mobility of heavy ions at a high electric field, and OI15 followed this approach tocompute the ion temperature Ti . In this study, we apply this approach to both Ti
and Te . Secondly, we neglect the impact ionization of neutral molecules by electricallyheated electrons by assuming that the electron energy in MRI turbulence is well below
the ionization potential of the neutrals (∼ 10 eV). The results of our calculations showthat this assumption holds in most parts of protoplanetary disks.
We denote the mean drift velocity and mean kinetic energy of a charged species α(= i for ions, e for electrons) by ⟨vα ⟩ and ⟨ϵα ⟩, respectively. In a weakly ionized gas
with an applied electric field E, the momentum and energy of the charged species

36

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

are determined by the balance between the neutral gas drag and acceleration by theelectric field (Hershey, 1939). Explicitly, the solution of the momentum and energybalance equations can be written as (Equations (A9) and (A10) of OI15)⟨vα ⟩ =

mα + mnqα E∆tα ,mα mn(3.13)

3(mα + mn )3(qα E∆tα )2 ,⟨ϵα ⟩ = kT +22(mα mn )2(3.14)

where qα , mα , and ∆tα are the charge, mass and mean free time of the plasma particles(e.g., qe = −e and qi = e, where e is the elementary charge). Since the magnetic field isneglected in this study, the mean drift velocity is parallel to the electric field. In a weakly
ionized gas, the plasma mean free time is determined by neutrals gas particles,∆tα = (nn ⟨σαn vαn ⟩)−1 ,(3.15)

where vαn is the relative velocity between a plasma particle and a neutral particle, and
σαn is the momentum-transfer cross section for the plasma–neutral collision. For electrons, σen is approximately constant at low energies (Yoon et al., 2008), and therefore
we may approximate ⟨σen ven ⟩ as σen ⟨ven ⟩. For ions, ⟨σin vin ⟩ is approximately constantowing to the polarization force between ions and neutrals (Wannier, 1953). Equations (3.13) and (3.14) are exact only when ∆tα is constant, but still hold in a good
accuracy even when ∆tα is velocity-dependent (Wannier, 1953).

The plasma temperature Tα is defined so that 3kTα /2 is equal to the kinetic energyof random motion, ⟨ϵα ⟩ − mα ⟨vα ⟩2 /2. Using Equations (3.13) and (3.14), Tα can be
written as(mα + mn )2(qα E∆tα )2 .
(3.16)
3km2α mn√
1/2
For electrons, we approximate ⟨ven ⟩ in ∆te with ⟨ve2 ⟩ = 3kTe /me . This allows us
Tα = T +to solve Equation (3.16) with respect to Te , and we obtain1
Te = T  +
2

√
1 2
+
4 3

(

E
Ecrit

)2


,

(3.17)

3.1. DISK AND IONIZATION MODELS
where√
Ecrit ≡

37

6me kT nn σenmne(3.18)

is the critical field strength above which electron heating becomes significant. We haveassumed me ≪ mn in deriving Equation (3.17). For ions, Equation (3.16) directly
gives
(
)2 )
2
kT
2(mi + mn )2 me σenE,= T 1+
m2i m2n⟨σin vin ⟩2 Ecrit
(
(
)(
)2 )
T
E
= T 1 + 7.6 × 10−7
,
100K
Ecrit
(

Ti

(3.19)

where we have set ⟨σin vin ⟩ = 1.6×10−9 cm3 s−1 (Nakano & Umebayashi, 1986) and σen =10−15 cm2 (Yoon et al., 2008) in the second expression, and used mi = 29 amu.

3.1.3

Ionization Balance and Accuracy of Simplified Approach

We calculate the plasma densities in a protoplanetary disk taking into account graincharging. The equations that describe the ionization balance in a dusty disk are (Equations (32), (33) and (35) of OI15)
ζnn − Krec (Te )ni ne − Kde (ϕ, Te )nd ne = 0,

(3.20)

ζnn − Krec (Te )ni ne − Kdi (ϕ, Ti )nd ni = 0,

(3.21)

ni − ne + Znd = 0,(3.22)

where ne and ni are, respectively, the number density of electrons and positive ions; Krecis the gas-phase recombination rate; Kde and Kdi are the adsorption rates of electrons
and ions onto grains; Z is the grain charge number; and ϕ is the coulomb potential on
grain surface. ϕ is related to Z as
ϕ=

eZ
.
a

(3.23)As the collisional frequency, Krec and Kde depend on the electron temperature Te , while
Kdi depends on the ion temperature Ti . Kde and Kdi also depend on the coulomb
potential of a grain surface ϕ. For HCO+ , the recombination rate Krec is given by

38

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)(Ganguli et al., 1988)(
Krec (Te ) = 2.4 × 10

−7

Te
300 K

)−0.69cm3 s−1 .

(3.24)

Approximating the ion velocity distribution by a Maxwellian with mean velocity ⟨vi ⟩and temperature Ti , Kdi is given by (Shukla & Mamun, 2002, OI15)
(
)
2
m
⟨v
⟩
2kT
i
i
i
exp −
Kdi (ϕ, Ti ) = πa2
πmi
2kTi
(
)]
(
)
kTi + 2e|ϕ|| ⟨vi ⟩ |+| ⟨vi ⟩ | 1 +erf √.
mi ⟨vi ⟩2
2kTi /mi[√(3.25)In this study, we also approximate the electron velocity distribution by a Maxwellian
with temperature Te . The drift velocity ⟨ve ⟩ can be neglected here since the drift speed√| ⟨ve ⟩ | is generally much smaller than the random speed ∼ kTe /me owing to thesmallness of me /mn (see Golant et al., 1980; Lifshitz & Pitaevskii, 1981). The electron
adsorption rate coeﬃcient Kde is given by the simple expression (Shukla & Mamun,
2002)√
Kde (ϕ, Te ) = πa2

)
 (
eϕ

 1+
, ϕ > 0,8kTe 
kTe
(
)
×
eϕ
πme 

 exp
, ϕ < 0.
kTe

(3.26)It should be noted that Equations (3.26) and (3.25) assume perfect sticking of ionsand electrons onto grain surfaces. This is a good approximation as long as the plasma
temperatures are well below 100 eV (see Section 3.2.2 of OI15 for more discussion).

Equations (3.20)–(3.22) determine ne , ni and Z at each location in a disk as a functionof E. We solve these equations using the procedure presented by Okuzumi (2009, their
Section 2.2; see also Section 3.2.4 of OI15).

To test the accuracy of our simplified approach, we reproduce the current–field relation including plasma heating (the nonlinear Ohm’s law of OI15) with adopting the
calculation steps in OI15. Current density is generally given by
J(E) = qe ne ⟨ve ⟩ + qi ni ⟨vi ⟩ .

(3.27)

3.1. DISK AND IONIZATION MODELS

39

Current Density [esu cm-2 s-1]

10-2
-3

10

-4

10

10-5
10-6
10-7
10-8 -11
10

10-10
10-9
10-8
10-7
Electric Field Strength [esu cm-2]

10-6

Figure 3.1: Test of the simplified plasma heating model presented in Section 3.1.3.
The solid curve shows the J–E relation for ‘model C’ of OI15 derived using the exact
electron velocity distribution (see Figure 10 of OI15), while the dashed curve shows our
reproduction based on the simplified approach.

Including plasma heating, the number densities depend on the electric fields strength
E. To obtain the current density, we first calculate plasma temperatures Te and Ti
from Equations (3.17) and (3.19) in an applied electric field E. We then calculate thenumber densities of plasma ne and ni from the ionization balance (Equation (3.22)).
We finally obtain the current density using Equations (3.13) and (3.27). In Figure 3.1,
we compare our result with the result of OI15 for the parameter set ‘model C’ of OI15.
We find that our calculation reasonably reproduces the previous result even at high fieldstrengths (E ≳ 10−9 esu cm−2 ) where electron heating is significant. The maximum
relative diﬀerence between the two results is 37%.

40

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

3.2

Active, Dead, and E-heating Zone

3.2.1

Conditions for MRI Growth

In the limit of ideal MHD, the criterion for the MRI is given by (Balbus & Hawley,1991)
λideal < H,

(3.28)vAzΩ(3.29)

where
λideal ≡ 2π

is the characteristic wavelength of the most unstable axisymmetric MRI modes, and√vAz = Bz / 4πρ and Bz are the vertical components of the Alfven velocity and magnetic field, respectively. Equation (3.28) expresses that the MRI operates when thelengthscale of the MRI modes is smaller than the vertical extent of the disk. Whenviewed as a function of z, λideal increases with z because ρ decreases toward the disk
surface. If we use Equation (3.4), the above MRI criterion can be rewritten in terms ofheight asz<√
2 ln (βc /8π 2 ) H ≡ Hideal ,(3.30)

where the height Hideal defines the upper boundary of the MRI active zone.

Inclusion of a finite Ohmic resistivity η introduces another criterion for MRI growth.
The criterion can be expressed in terms of the Elsasser number (Turner et al., 2007)
Λ≡

2
vAz
.
ηΩ

(3.31)The instability grows when
Λ>1and decays when Λ < 1 (e.g., Sano & Miyama, 1999).

(3.32)

3.2. ACTIVE, DEAD, AND E-HEATING ZONE

3.2.2

41

Zoning Criteria

Here we describe how to determine turbulent state at a position in protoplanetarydisks. Electron heating aﬀects on the MRI turbulence when the ionization fraction is
suﬃciently decreased. We express the condition that the heating takes place and aﬀect
MRI turbulence, and then summarize three turbulent states of MRI and steps of zoning
a disk into the state.
For electron heating to take place, the field must be suﬃciently amplified before MRI
turbulence reaches a fully developed state that means the stop of MRI growth. Muranushi et al. (2012) performed a local unstratified resistive MHD simulation and foundthat the fully developed current density is√
Jmax = fsat

ρ
cΩ,
2π

(3.33)

where fsat ≈ 10 according to the results by Muranushi et al. (2012). Here, we assume
fsat to be fsat = 10 and the maximum current density is Jmax . Thus, when the current
density reaches Jmax before electric field reaches the criterion for electron heating Ecrit ,MRI turbulence does not cause the electron heating.
As we will describe later in this section, we use current density to decide whether electron heating take place or not. Therefore, we transform the condition for suppressing
MRI into a form using current density. We adopt Λ = 1 (Equation (3.32)) as the criterion for suppressing MRI which is triggered by electron heating. Using the electric
conductivity σc and the relation η = c2 /4πσc , the condition for sustaining MRI turbu2
lence Λ ≳ 1 leads to a condition σc ≳ c2 Ω/(4πvAz). Under the Ohm’s law J(E) = σc E,

the condition can be rewritten as a lower limit to the current density
J(E) ≳ JΛ=1 (E),(3.34)

whereJΛ=1 (E) ≡ σc (Λ = 1)E =

c2 Ω
E.
2
4πvAz

(3.35)

Using the above criteria, we can classify a region in protoplanetary disks into three

42

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)
Select position (r, z)!
Set E = 0!
MRI Growth ?!
(Λ > 1 ?)!

Dead zone!
NO!

YES!

(No MRI turbulence)!

Increase E!
MRI Suppressed ?!
(J(E) = JΛ=1 ?)!
NO!

E-heating zone!
YES! (Self-regulated MRI turbulence)!

MRI Fully Developed ?
(J(E) = Jmax ?)NO!

Active zone!
YES! (Fully developed MRI turbulence)!

Figure 3.2: Flow chart showing key steps of zoning a protoplanetary disk into thedead, active, and e-heating zones.
diﬀerent zones corresponding to three turbulent state of MRI.
1. Dead zone. Because of the low ionization fraction, Ohmic dissipation suppress all
the unstable MRI mode. Suppressed MRI does not generate turbulence and also
current density. We will refer to the region where MRI is completely suppressedas the “dead zone”. In this case, the condition of Ohmic dissipation (Equation(3.35)) is satisfied with no MRI turbulence.
2. E-heating zone. Electric fields of MRI turbulence become suﬃciently high forelectron heating to be caused. The Ohmic dissipation is amplified by the electron
heating after the MRI grows. We will refer to the region where electron heatingaﬀects MRI turbulence as the “e-heating zone”, where the “e” refers to both
“electric field” and “electron.” In this case, current density falls down the criticalcurrent density of Ohmic dissipation (Equation (3.35)).
3. Active zone. MRI sustains fully developed turbulent state because the gas is
suﬃciently ionized so that Ohmic dissipation is not eﬃcient. We will refer to theregion where vigorous MRI turbulence is sustained as the “active zone” in thisstudy. In this case, the current density J reaches and sustains its maximum value
Jmax before electron heating reduces the MRI turbulence.
We summarize the calculation steps for zoning the disk region under some assump-3.3. LOCATION OF THE E-HEATING ZONE

43

tions. We assume that the electric field strength correspond to the activity of MRI
turbulence since developed MRI generates strong electric fields. The growth of MRI
implies increasing electric fields, and the decay of MRI implies decreasing electric fields.
Furthermore, we also assume that magnetic fields are not varied by the MRI growth forsimplicity. Under these assumptions, we determine the turbulent state at the position
with following steps (see Figure 3.2): First, we select a calculated position in the region
satisfying Equations (3.30) of a disk. We then calculate values at the position withsetting E = 0. When MRI is initially suppressed by Ohmic dissipation, i.e., Λ < 1 at
E = 0, the positions belong to the dead zone. During satisfying unstable condition,i.e., Λ > 1, the electric field strength E is increased from E = 0 with iterating until the
turbulent state at the position is determined. We calculate current density J(E) andassess some conditions in E. When MRI turbulence causes electron heating and Ohmic
dissipation become eﬃcient, i.e., J(E) = JΛ=1 , the position belongs to the e-heatingzone. When MRI turbulence is fully developed, i.e., J(E) = Jmax , the position belongs
to the active zone. We conduct the above steps in the whole region in a disk, and zone
a protoplanetary disk into the dead, active, and e-heating zones.

3.3

Location of the E-heating Zone

We here predict the location of the e-heating zone in protoplanetary disks using themethodology described in Section 3.2.2. We conduct a parameter study varying the
midplane plasma beta βc , grain size a, dust-to-gas mass ratio fdg , and surface density
scaling factor fΣ . Following Sano et al. (2000), we select the MMSN (fΣ = 1 andq = 3/2) with a = 0.1 µ, fdg = 0.01, and βc = 1000 as the fiducial model. We start out
with this fiducial model in Section 3.3.1, and discuss the dependence on the parametersin the subsequent subsections. A summary of the parameter study is given in Table 3.1.
We also describe ion heating in Section 3.3.5.

44

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

Table 3.1: Sizes of the Dead and E-heating Zones for Various Parameter Setsβc

a (µm)

fdg

102
103
104
105
103
103
103
103
103
103
103
103
103
103
103
103

0.1
0.1
0.1
0.1
0.1
1
10
100
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1

10−2
10−2
10−2
10−2
10−2
10−2
10−2
10−2
10−1
10−2
10−3
10−4
10−2
10−2
10−2
10−2

0

10

fΣ

Outer radius (AU)
Dead zone E-heating zone
1
18
74
1
24
82
1
34
82
1
56
82
1
24
82
1
11
39
1
8
19
1
8
11
1
52
151
1
24
82
1
12
41
1
8
20
10
55
149
3
36
114
1
24
82
0.3
14
44

Height [AU]

9
6
3
0

20
30
40
50
60
70
Distance From Central Star [AU]80

90

Figure 3.3: Cross-section view of the fiducial protoplanetary disk indicating the locationof the dead, e-heating and active zones (red, green-shaded, and blue regions, respectively). The dashed line shows the gas scale height H, while the dotted line shows thecritical height Hideal below which the MRI criterion in the ideal MHD limit is satisfied(see Equation (3.30)).

3.3. LOCATION OF THE E-HEATING ZONE

3.3.1

45

Fiducial Disk Model

Figure 3.3 shows the two-dimensional (radial and vertical) map of the dead, active, ande-heating zone in the fiducial disk model. The MRI criterion in the ideal MHD limit(Equation (3.28)) is satisfied at altitudes below z = Hideal ≈ 2.3H (see Equation (3.30)).
The region above this height is MRI-stable with the MRI modes suppressed by too
strong magnetic tension. The dead zone is located inside 24 AU from the central starand near the midplane where the gas is shielded from ionizing irradiation. The sizeof the dead zone for this disk model is consistent with the prediction by Sano et al.
(2000) (see their Figure 7(b)), although their dead zone is slightly thicker than ours
because of the neglect of X-ray ionization. We find that the e-heating zone extends fromthe outer edge of the dead zone out to 82 AU from the central star. This means that
MRI turbulence can develop without aﬀected by electron heating only in the outermost
region of r ≳ 80 AU.
To illustrate how our zoning criteria work in this particular example, we plot in Figure 3.4 the relation between the current density J and electric field E in the midplaneat 15 AU, 45 AU and 90 AU, which represent the dead, e-heating, and active zones,
respectively. Recall that for fixed E, MRI turbulence grows if J(E) > JΛ=1 and decays otherwise (Equation (3.34)). At 15 AU, J(E) falls below JΛ=1 for all values of E,
implying that the MRI is unable to grow at this location. At 45 AU, the MRI growth
condition is satisfied during the initial growth stage of E ≪ 10−11 esu cm−2 , but breaks
down before J reaches Jmax because of the decrease in J(E) due to electron heating.
This implies that MRI turbulence is allowed to grow in the initial stage but saturates
at a level lower than that for fully developed turbulence. At 90 AU, J(E) reaches
Jmax before electron heating sets in, implying that fully developed MRI turbulence issustained here.
In contrast to electron heating, ion heating is found to be negligible at all locationsin the fiducial disk model. In the e-heating zone, the electric field strength at thesaturation point is typically ≲ 102 Ecrit (see the center and right panels of Figure 3.4),
which is an order of magnitude lower than the field strength required for ion heating,
103 Ecrit .

46

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

10-4

10

Current Density J [esu cm-2 s-1]

Current Density J [esu cm-2 s-1]

-3

-4

10

-5

10

-6

10

-7

10

-8

10-5
10-6
10-7
10-8
10-9

10

10-13 10-12 10-11 10-10 10-9

10-8Electric Field Strength E [esu cm-2]

10-14 10-13 10-12 10-11 10-10 10-9
Electric Field Strength E [esu cm-2]-2 -1

Current Density J [esu cm s ]

10-5
10-6
10-7
10-8
10-9
10-10
10-14 10-13 10-12 10-11 10-10 10-9
Electric Field Strength E [esu cm-2]

Figure 3.4: Relations between the current density J and electric field strength E in
the midplane at 15 AU (upper left panel), 45 AU (upper right panel), and 90 AU (lower
panel), which represent the J–E relations in the dead, e-heating, and active zones,
respectively. The thick solid line shows the current–field relation J(E), the dotted line
the maximum current density of MRI, Jmax (Equation (3.33)), the vertical gray solidline the criterion for electron heating, Ecrit (Equation (3.18)), and the dashed line thecritical current density JΛ=1 below which the MRI decays owing to Ohmic dissipation(Equation (3.35)). The black dots on the J–E relations indicate the saturation points
at which either fully developed (J(E) = Jmax ) or self-regulated (J(E) = JΛ=1 ) MRI
turbulence is sustained.

8

Height [AU]420Height [AU]βc = 102

6

0

8

10

20
30
40
50
60
70
Distance From Central Star [AU]80βc = 104

6
4
2
0

0

10

20
30
40
50
60
70
Distance From Central Star [AU]

80

8

90

47
βc = 103

6
4
2
0

90

Height [AU]Height [AU]3.3. LOCATION OF THE E-HEATING ZONE

0

8

10

20
30
40
50
60
70
Distance From Central Star [AU]

80

90

80

90

βc = 105

6
4
2
0

0

10

20
30
40
50
60
70
Distance From Central Star [AU]Figure 3.5: Same as Figure 3.3, but for diﬀerent values of the midplane beta βc .

3.3.2

Dependence on the Magnetic Field Strength

Figure 3.5 shows how the size of the dead and e-heating zones depend on the midplane
plasma beta βc . Recall that a higher βc corresponds to a weaker magnetic field B
threading the disk. As we increase βc , the dead zone expands because the Elsasser
number Λ ∝ B 2 decreases. On the other hand, we find that the boundary between thee-heating and active zones is less sensitive to the choice of βc . As can be inferred fromthe middle and right panels of Figure 3.4, this boundary is approximately determinedby the condition that the current density J(E) reaches Jmax at a local maximum lying
at E ≈ Ecrit . Since both Ecrit and Jmax are independent of B and hence of βc , so is theboundary between the e-heating and active zones.

3.3.3

Dependence on the Grain Size and Dust-to-Gas Mass
Ratio

The size and amount of dust grains in disks are important parameters in the ionization
model as they eﬃciently remove plasma particles from the gas phase. Obviously, thesequantities change as the grains coagulate, settle, or are incorporated by even larger
solid bodies like planetesimals. We here explore how the change of these parametersaﬀect the size of the dead and e-heating zones.
To begin with, we show in Figure 3.6 the location of the dead, active, and e-heating zones
with the dust-to-gas ratio fdg fixed to 0.01 but with the grain size a varying between

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)8

Height [AU]420Height [AU]a = 0.1 µm

6

0

8

10

20
30
40
50
60
70
Distance From Central Star [AU]80

a = 10 µm

6
4
2
0

0

10

20
30
40
50
60
70
Distance From Central Star [AU]

80

90

8

a = 1.0 µm

6
4
2
0

90

Height [AU]Height [AU]

48

0

8

10

20
30
40
50
60
70
Distance From Central Star [AU]

80

90

80

90

a = 100 µm

6
4
2
0

0

10

20
30
40
50
60
70
Distance From Central Star [AU]Figure 3.6: Same as Figure 3.3, but for diﬀerent values of the grain size a.

0.1 µm and 100 µm. We can see that the e-heating zone shrinks with increasing grain
size. On increasing a by a factor of 10, the outer radius of the e-heating zone decreases
by a factor of ≈ 2. Qualitatively, this is simply because the ionization fraction of thegas increases with decreasing total surface area of the grains. Equation (3.20) showsthat the electron abundance xe = ne /nn in equilibrium is inversely proportional to thetotal surface area of grains per unit volume 4πa2 nd as long as adsorption of plasma
particles onto the grains dominate over gas-phase recombination. When dust grains
aggregate, their total surface area decreases inversely proportional to a, and hence theelectron abundance increases linearly with a. The resulting increase in the electric
conductivity causes a shift of the J–E curve toward higher J, enabling the curve tocross the J = Jmax line at smaller orbital radii. We also find that the outer radius of
the dead zone decreases at a similar rate to that of the e-heating zone when we go froma = 0.1 µm to 1 µm. However, the decrease in the dead zone size stops beyond this grain
size because gas-phase recombination takes over plasma adsorption onto dust grains.
As a consequence, the e-heating zone becomes narrower and narrower as a increases
beyond 10 µm.
Decreasing the dust-to-gas mass ratio fdg has a similar eﬀect to increasing the grainradius because the total surface area of the grains is linearly proportional to fdg . Thiscan be seen in Figure 3.7, where we show the location of the dead and e-heating zones
for a = 0.1 µ with fdg varying between 10−1 and 10−4 . We see that the outer radiiof the active and e-heating zone decrease by a factor of ≈ 2 when fdg is decreased bya factor of 10. This trend is similar to what we have seen when increasing the grain

3.3. LOCATION OF THE E-HEATING ZONE
16

fdg = 1012-1Height [AU]Height [AU]

16

8
4
0

0

20

40
60
80
100
120
Distance From Central Star [AU]

140

-2

8
4
0

20

40
60
80
100
120
Distance From Central Star [AU]

140

160

140

160

16

fdg = 1012-3Height [AU]Height [AU]fdg = 10

12

0

160

16

8
4
0

49

0

20

40
60
80
100
120
Distance From Central Star [AU]

140

-4

8
4
0

160

fdg = 10

12

0

20

40
60
80
100
120
Distance From Central Star [AU]Figure 3.7: Same as Figure 3.3, but for diﬀerent values of the dust-to-gas mass ratio
fdg .
16
Height [AU]Height [AU]16fΣ = 10

12
8
4
0

0

20

40
60
80
100
120
Distance From Central Star [AU]

140

4
0

20

40
60
80
100
120
Distance From Central Star [AU]

140

160

40
60
80
100
120
Distance From Central Star [AU]

140

160

16
Height [AU]Height [AU]16fΣ = 1

12
8
4
0

8

0

160

fΣ = 3

12

0

20

40
60
80
100
120
Distance From Central Star [AU]

140

160

fΣ = 0.3

12
8
4
0

0

20

Figure 3.8: Same as Figure 3.3, but for diﬀerent values of the surface density scaling
factor fΣ .
radius by the same factor.

3.3.4

Dependence on the Disk Mass

Finally, we examine how the size of the e-heating zone depends on the disk mass. Figure
3.8 shows the location of the e-heating zone for diﬀerent values of fΣ . Here, we fix thedust-to-gas mass ratio fΣ so that both the gas and dust densities scale with fΣ . We
find that the e-heating zone expands toward larger orbital radii and higher altitudes asfΣ increases. In the horizontal direction, the expansion is mainly due to the increased
amount of dust grains with increasing fΣ . As we have explained in 3.3.3, the ionization
fraction of the gas scales inversely with 4πa2 nd , and hence with fΣ . Therefore, increasingfΣ by a factor has the same eﬀect as increasing fdg by the same factor as long as theCHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)!&!!!"!#$%!

!(

!%

!%

!$

!"

!#

!"

</)=>+!789:!"!B!!

?3-@A/!;<=>

!+

!#"

!$" !%" !&" !'" !(" !)",-./0123!4567!831/509!:/05!;<=>

!*"

!+"
!$

!';
"#$!"!#$&#!

!'"

!?
!#

!;
!"

!"!@!!

50

!'

!"

!#"

!$"
!%"
!&" !'"" !'#"()*+,-./!0123!4/-+1,5!6+,1!789:

!'$"

!'%"

Figure 3.9: Two-dimensional distribution of the ion temperature Ti for the βc = 100
model (upper panel) and fdg = 0.1 model (lower panel). The solid lines show theboundary of the e-heating zone, while the dotted lines show Hideal .
ionization rate ζ is unchanged (which is approximately true at ∼ 100 AU where cosmic
rays penetrate down to the midplane). This is exactly what we see in Figures 3.7 and
3.8, where the e-heating zone expands to 150 AU when either fdg or fΣ is increasedby the factor of 10 from the fiducial value. By contrast, the vertical expansion of thee-heating zone is caused by the attenuation of X-rays that occurs at higher altitudes
with increasing gas column density.

3.3.5

Ion Heating

We observe ion heating in two cases where βc = 100 and where fdg = 0.1. Figure 3.9
plots the distribution of the ion temperature Ti in the saturated state for these cases.
In the case of βc = 100 (the upper panel of Figure 3.9), Ti is 3–4 times higher than thetemperature in a region slightly outside the e-heating zone. In this case, the Elsasser
number Λ exceeds unity even after electron heating reduces Λ. This allows the electricfield strength to reach the critical value for ion heating (≈ 103 Ecrit ) in the vicinity ofthe e-heating zone. In the case of fdg = 0.1 (the lower panel of Figure 3.9), ion heatingtakes place near the upper boundary of the e-heating zone. However, the region is very
narrow, and the temperature rise is less than 2T . Therefore, in this case, ion heating

100

100

10-1

10-1

10-2

10-2

αMRI

αMRI

3.4. SATURATION OF TURBULENCE IN THE E-HEATING ZONE

10-3
10-4
10-5

51

10-3
10-4

0

10 20 30 40 50 60 70 80 90
Distance From Central Star [AU]

10-5

0

10 20 30 40 50 60 70 80 90
Distance From Central Star [AU]Figure 3.10: Radial distribution of αMRI (Equation (3.40)) for the fiducial model(left panel) and βc = 104 (right panel). The solid black line shows αMRI,mid including
electron heating on the mid-plane, and the solid blue line shows ᾱMRI including electron
heating integrated in the z-direction. The dashed black line shows αMRI,mid without
including electron heating on the mid-plane, and the dashed blue line shows ᾱMRI
without including electron heating integrated in the z-direction.
might be practically negligible.

3.4

Saturation of Turbulence in the E-heating Zone

We have shown in Section 3.3 that self-regulation of the MRI due to electron heating
can occur over a large region of protoplanetary disks. Then the question arises how
strongly the e-heating will suppress the MRI turbulence in the e-heating zones. This
question can only be fully addressed with MHD simulations including magnetic diﬀusionand electron heating in a self-consistent manner, which is far beyond the scope of thisstudy. In this section, we attempt to estimate the saturation level of MRI turbulencefrom simple scaling arguments.
As usual, we quantify the strength of turbulence with the Shakura–Sunyaev α parameter
α = Trϕ /P , where P = ρc2s is the gas pressure and Trϕ is the rϕ component of turbulent
stress. In MRI-driven turbulence, Trϕ is generally dominated by the turbulent Maxwell
stress −δBr δBϕ /4π (Hawley et al., 1995; Miller & Stone, 2000), where δBr and δBϕ
are the radial and azimuthal components of the turbulent (fluctuating) magnetic fields.

52

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

Therefore, we evaluate the α parameter for MRI turbulence as
αMRI ≈ −δBr δBϕ.
4πρc2s(3.36)In reality, the Reynolds stress (Fleming & Stone, 2003; Okuzumi & Hirose, 2011) or thecoherent component of the Maxwell stress (e.g., Turner & Sano, 2008; Gressel et al.,2011) can dominate over the turbulent Maxwell stress at locations where the MRI is
significantly suppressed. However, we do not include these components in our αMRI
because they do not reflect the local MRI activity at such locations (see the references
above).
Next we relate the amplitude of turbulent magnetic fields to the amplitude of the electric
current density J = |J | using the Ampére ’s law J = (c/4π)∇ × δB. We neglect largescale, coherent components in B since the electric current is inversely proportional to
the length scale of fields. We assume that the magnetic field in MRI-driven turbulence
is dominated by the azimuthal component δBϕ and varies over a length scale ∼ λideal ,where λideal is the wavelength of the most unstable MRI modes already introduced in
Equation (3.29). Then, from the Ampére ’s law, one can estimate the magnitude of thecurrent density asc|∇ × B|
4π
√
c vAz
ρ
δBϕ
≈
δBϕ =
cΩ
,
4π Ω
4π
Bz

J =

(3.37)

where we have replaced the derivative ∇ with wavenumber 2π/λideal = Ω/vAz . If we
use the maximum current Jmax for fully developed MRI turbulence (Equation (3.33)),
Equation (3.37) results in a simple scaling relation√ J
δBϕ≈ 10 2.
Bz
Jmax(3.38)For fully developed MRI turbulence where J ≈ Jmax , the above equation predicts
δBϕ /Bz ∼ 10, in agreement with the results of MHD simulations (e.g., Hawley et al.,
1995; Sano et al., 2004).
Now let us consider situations where e-heating is so eﬀective that the growth of the MRI

3.4. SATURATION OF TURBULENCE IN THE E-HEATING ZONE

53

is saturated at J ≈ JΛ=1 ≪ Jmax . Assuming δBz ≲ Bz0 for this case, we have√JδBϕ ≈ 10 2Bz0.
Jmax

(3.39)

This equation predicts the amplitude of δBϕ as a function of Bz0 and J/Jmax . MHD
simulations show that δBr ≈ −(0.4 . . . 0.6)δBϕ in MRI turbulence (Hawley et al., 1995;
Sano et al., 2004). Assuming that this scaling also holds in our case, we have δBr δBϕ ≈2−100Bz0(J/Jmax )2 . Finally, substituting this into Equation (3.36), we obtain the scaling

relation between αMRI and J/Jmax ,
αMRI

(
)2
100B02J≈4πρc2s Jmax
(
)−1 (
)2
β0
J
≈ 0.2
,
1000
Jmax

(3.40)

2
where β0 ≡ 8πρc2s /Bz0
= βc exp (−z 2 /2H 2 ) is the plasma beta (not necessarily at the

midplane) associated with the net vertical field Bz0 . Formally, the derivation leading
to Equation (3.40) breaks down when MRI is so active that δBz ≫ Bz0 and J ≈Jmax . Nevertheless, we find that Equation (3.40) reproduces the results of ideal MHD
simulations with a reasonably good accuracy. Equation (3.40) predicts that αMRI ≈ 2
for β0 = 102 and αMRI ≈ 0.02 for β0 = 102 when J = Jmax . These are consistent withthe results of isothermal simulations by Sano et al. (2004) showing that the Maxwell
component of α is ∼ 1 for β0 = 102 and ∼ 0.01 for β0 = 104 (see their Table 2, column(10)). Therefore, we will apply Equation (3.40) to both the e-heating zone and activezone.
The left panel of Figure 3.10 show the radial distribution of αMRI for the fiducial disk
model predicted from Equation (3.40). Here we plot the midplane value αMRI,mid ≡αMRI (z = 0) and the density-weighted average in the vertical direction,∫ HidealᾱMRI ≡

−HidealαMRI (z ′ )ρ(z ′ )dz ′Σ,(3.41)

where we have assumed αMRI = 0 in the magnetically dominated atmosphere at |z| >Hideal . The former quantity measures the MRI activity at the disk midplane, while thelatter quantity is more closely related to the vertically integrated mass accretion rate

54

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)(Suzuki et al., 2010). For the fiducial disk model, we find that αMRI,mid ∼ 10−5 and 10−3
at the inner and outer edge of the e-heating zone (20 AU and 80 AU), respectively. These
values are more than two orders of magnitude lower than the value αMRI,mid = 0.2 inthe active zone (r ≳ 80 AU). This implies that the MRI is “virtually dead” deep inside
the e-heating zone. We also find that αMRI,mid changes discontinuously at the boundary
between the e-heating and active zones. The reason is that when the saturated statechanges at the point, J/Jmax also changes from unity to one order of magnitude becauseof the N-shaped current–field relation (see middle and right panels of Figure 3.4). The
vertical average ᾱMRI decreases more slowly with decreasing r, because the upper layerof the disk remains MRI-active (see Figure 3.3). This picture is qualitatively similar tothe classical layered accretion model of Gammie (1996). In right panel of Figure 3.10,
we also plot the radial distribution of αMRI,mid and ᾱMRI for a disk with βc = 104 . We
find that αMRI,mid in e-heating zone is almost unchanged from the fiducial disk. The
reason is that increase of (J/Jmax )2 ≈ 10 cancels out the depletion of βc−1 ≈ 10−1 inEquation (3.40). Therefore, αMRI,mid remains low saturation level.
In summary, our simple estimate predicts that MRI turbulence can be significantly
suppressed in the e-heating zone. In this sense, the e-heating zone acts as a extended
dead zone. However, our estimate relies on the hypothetical scaling between the andturbulent Maxwell stress and J/Jmax , which is as yet justified by MHD simulations. 1
In order to test our prediction, we will perform resistive MHD simulations including
electron heating in future work.

3.5

Charge Barrier against Dust Growth in the Eheating Zone

So far we have focused on the role of electron heating on the saturation of MRI turbulence. As pointed out by OI15, electron heating also has an important eﬀect on the1

However, there are some support for Equation (3.40) from MHD simulations including ambipolar
diﬀusion, not Ohmic dissipation. Bai & Stone (2011) reported the Maxwell component of α (their Table2) and the cumulative probability distribution of J (Figure 6) for three simulation runs with β0 = 400
and with diﬀerent values of ambipolar diﬀusivity. Their results show that αMaxwell ≈ 0.17, 0.029, and
0.0041 for models with J/Jmax ≈ 1, 0.3, and 0.1 (median values), respectively. These are consistent
with Equation (3.40) predicting that αMRI ≈ 0.5, 0.045, and 0.005 for these values of J/Jmax .

3.5. CHARGE BARRIER AGAINST DUST GROWTH IN THE E-HEATING ZONE55
growth of dust grains. In an ionized gas, dust grains tend to be negatively charged
because electrons collide and stick to dust grains more frequently than ions. The resulting Coulomb repulsion slows down the coagulation of the grains through Brownian(thermal) motion. This “charge barrier” is also present in weakly ionized protoplanetary disks, in which dust grains tend to be charged as in a fully ionized gas when theirsize is larger than 1 µ (Okuzumi, 2009; Matthews et al., 2012). The important role ofelectron heating in this context is that heating electrons further promote the negative
charging of the grains, because the grain charge in a plasma is linearly proportional tothe electron temperature (e.g., Shukla & Mamun, 2002). In this section, we explore
how this aﬀects dust coagulation in the e-heating zone.

For simplicity, let us assume that dust grains have the single radius a and charge Z.
The grains can collide with each other if the condition
Ecol > Eelc

(3.42)

is satisfied (Okuzumi, 2009). Here, Ecol is the kinetic energy of the relative motion of
two colliding grains, andEelc ≈(eZ)22a(3.43)

is the Coulomb repulsion energy of the grains just before contact. We focus on small
dust grains near the midplane and assume that the relative motion is dominated byBrownian motion and turbulence-induced motion. Then, the kinetic energy of relative
motion can be expressed as
Ecol = EBrown + Eturb ,(3.44)

where EBrown and Eturb are the kinetic energy of Brownian motion and turbulenceinduced motion, respectively. Brownian motion is the thermal motion of grains, andEBrown is approximately expressed as
1
EBrown ≈ µu2th ,2
where the thermal velocity of grains uth is expressed as uth =(3.45)√
8kT /πm and the

reduced mass of grains µ is expressed as µ = m2 /(m + m) = m/2. The relative energy

56

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)of turbulence-induced motion is expressed as
1
Eturb ≈ µ(∆uturb )2 ,2

(3.46)

where ∆uturb is the relative velocity of the grains excited by turbulence. For small
grains, ∆uturb is approximately given by (Weidenschilling, 1984; Ormel & Cuzzi, 2007)∆uturb ≈√αdisp Re1/4 cs Ωτs ,(3.47)

where αdisp ≡ ⟨δv 2 ⟩ /c2s is the velocity dispersion of the gas ⟨δv 2 ⟩ normalized by c2s , Re
is the Reynolds number of turbulence, and√τs = ρ• a/( 8/πcs ρ)(3.48)

is the stopping time of the grains (we have adopted Epstein’s drag law for τs ). The
Reynolds number is expressed as Re = αdisp c2s Ω−1 /νmol , where νmol is the molecularviscosity. We estimate α with and without electron heating, using Equation (3.40)
presented in Section 3.4. Turbulence dominates the collisional energy when αdisp ishigh and/or a is large. For the moment, we simply assume αdisp = αMRI , where αMRI
is the normalized local Maxwell stress introduced in Equation (3.36). This assumption
holds when the Reynolds stress in the e-heating zone is comparable to the Maxwell
stress. In reality, the Reynolds stress in the e-heating zone might be higher than theMaxwell stress for a reason discussed later. Therefore, the estimate of Eturb presented
here should be taken as a lower limit.
To obtain Z and αMRI , we calculate the ionization fraction (Section 3.1.3), determine
the turbulent state (Section 3.2.2), and estimate the MRI-turbulent viscosity (Section
3.4) with changing grain radius a at a location. We then obtain Ecol and Eelc by abovementioned method. It should be noted that grains have single size and changing grainradius means changing the size of all grains at the location. Thus, the turbulent stateat the location also depends on a.
In Figure 3.11, we plot the ratio Eelc /Ecol as a function of a at 35 AU in the midplane forthe fiducial disk model. The ratio quantifies the eﬀectiveness of the charge barrier: the3.5. CHARGE BARRIER AGAINST DUST GROWTH IN THE E-HEATING ZONE572

10E
elc //E
col
Eelc
Ekin

101
0

10

10

-1

10-2

10

-3

0.1

1
Grain Size [µm]Figure 3.11: Eﬀectiveness of the charge barrier against grain growth as a function of
the grain size at the midplane 35 AU in the fiducial model. The solid line (red) shows
Eelc /Ecol including electron heating, and the dashed line (blue) shows Eelc /Ecol without
including electron heating. The horizontal dotted line shows Eelc /Ecoll = 1, above whicha strong Coulomb repulsion between the grains suppresses their mutual collision crosssection. Here it is assumed that αdisp (= ⟨δv 2 ⟩ /c2s ) is equal to αMRI , the normalized
local Maxwell stress given by Equation (3.36) (but see also Figure 3.12).

collisional cross section of two equally charged grains is significantly suppressed when
Eelc /Ecol ≫ 1. We find that electron heating significantly enhances the charge barrierfor submicron-sized grains. If electron heating is not included, this location belong to
the dead zone and the active zone with grain size being ≲ 0.05 µm and ≳ 0.05 µm,
respectively. In this case, Eelc /Ecol is much lower than unity in all a. Thus we canconclude that dust grains at this location can grow without the charge barrier. On theother hand, if electron heating is included, this location belongs to the e-heating zone
when 0.05 µm ≲ a ≲ 1.4 µm (see also Figure 3.6). In the e-heating zone, grains arecharged by heated electrons, leading to increase of Eelc , and MRI turbulence as collisionalsource is well suppressed, leading to decrease of Eturb . Consequently, Eelc /Ecol is larger
than unity when 0.08 µm ≲ a ≲ 0.5 µm. In particular, Eelc /Ecol takes its maximum
value of 40 at a = 0.2 µm corresponding to EBrown = Eturb . Both the suppression of
turbulence and grain charge would enhance the charge barrier.

58

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

10

2

E
elc /EcolEelc/ Ekin

101
100

10

-1

10-2
10-3

0.1

1
Grain Size [µm]Figure 3.12: Same as Figure 3.11, but we here evaluate αdisp = ⟨δv 2 ⟩ /c2s using Equation (3.49) assuming that sound waves propagate from upper MRI active layers to the
midplane as it is observed for the conventional dead zone.
There are at least two mechanisms that could drive further growth of dust in the eheating zone. One is vertical turbulent mixing of dust particles as already pointed
out by Okuzumi et al. (2011). In general, the charge barrier is less significant at higher
altitudes where dust particles have a higher collision energy due to vertical settling (anddue to if MRI is active there). Electron heating, which was not considered by Okuzumi
et al. (2011), does not change this picture because it is also ineﬀective at high altitudes.
Micron-sized grains in the e-heating zone can easily be lifted up to such high altitudes if
only weak turbulence is present there (Turner et al., 2010, see also dust scale height Hd
in Section 3.6.1). The lifted grains are allowed to collide and grow there until they fallback to the e-heating zone. In this way, small grains in the e-heating zone are able tocontinue growing on a timescale much longer than vertical diﬀusion timescale. Okuzumi
et al. (2011) showed that the charge barrier is overcome on a timescale of 105 –106 yr,
but they did not consider the amplification of grain charging due to electron heating.
How much the growth is delayed in the presence of electron heating should be studiedin future work.
Another potentially important mechanism is dust stirring by random sound waves. It

3.5. CHARGE BARRIER AGAINST DUST GROWTH IN THE E-HEATING ZONE59
is known that the Reynolds stress in a dead zone exceeds the Maxwell stress because ofsound waves propagating from upper MRI-active layers (e.g., Fromang & Papaloizou,
2006; Turner et al., 2010; Okuzumi & Hirose, 2011). If this is also the case for our
e-heating zone, the assumption αdisp = αMRI would significantly underestimate theparticle collision energy in the e-heating zone. To estimate this eﬀect, we now calculate αdisp using an empirical formula for the gas velocity dispersion in the dead zone
(Okuzumi & Hirose, 2011),
⟨

δv

2

⟩

(
≈ 0.78ᾱMRI c2s exp

)
z2
,
2H 2

(3.49)

where ᾱMRI is the density-weighted vertical average of αMRI defined by Equation (3.41).
Equation (3.49) expresses the amplitude of random sound waves inside a dead zone.
Figure 3.12 shows Eelc /Ecol in this case and is obtained in the same way as in Figure
3.11 but we here use Equation (3.49) for αdisp in Equation (3.47). The use of Equation (3.47) for the sound wave-driven collision velocity assumes that the time correlationof the waves’ velocity fluctuations exponentially decays on the timescale of Ω−1 as in
the Kolmogorov turbulence. We find that Eelc /Ecol now falls below unity at all grain
sizes. Thus, sound waves traveling from MRI-active layers, if they exist, could help
dust overcome the charge barrier in the e-heating zone. However, the argument madehere is not conclusive because the induced collision velocity depends on the assumed
time correlation function, or equivalently power spectrum, of the random sound waves.
If the power spectrum of the waves has only a small amplitude at high frequencies(to which small dust particles are sensitive) compared to the turbulent spectrum, thewave-induced collision velocity would be lower than given by Equation (3.47). The spectrum of velocity fluctuations in the e-heating zone should be studied in future MHD
simulations.

60

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

3.6

Discussion

3.6.1

Dust Diﬀusion

We have assumed so far that the dust-to-gas mass ratio is vertically constant. This
assumption breaks down when dust particles settle toward the midplane. If this is
the case, the dust-to-gas ratio would decrease at high altitudes, and consequently thee-heating zone would shrink in the vertical direction as expected from Figure 3.7.

However, as we will show below, dust settling is negligible even in the e-heating zone
because even weak turbulence is able to diﬀuse small grains to high altitudes. Youdin& Lithwick (2007) analytically derived dust scale height Hd in the sedimentationdiﬀusion equilibrium. If the particle stopping time τs is much smaller than the Keplerian timescale Ω−1 , which is true for small particles, the dust scale height can beapproximately written as(
Hd ≈ H 1 +Stαdisp,z)−1/2,(3.50)

where St = τs Ω is the so-called Stokes number and αdisp,z = ⟨δvz2 ⟩ /c2s is the verticalcomponent of the velocity dispersion normalized by c2s . Equation (3.50) implies that
dust settling takes place (Hd < H) when St > αdisp,z . Under the disk model employedin this study, St can be expressed as(−8

St = 3 × 10a0.1 µm

)
fΣ−1

( 2 )
( r )3/2
z
exp
.
1 AU
2H 2

(3.51)

Therefore, for a = 0.1 µ, dust settling in the e-heating zone (r ∼ 10–100 AU) occurs
only if αdisp,z ≲ 10−5 –10−6 . In the e-heating zone, αMRI ∼ 10−5 –10−3 at the midplane(see Figure 3.10), and therefore we may safely neglect dust settling even if the Reynolds
stress is as small as the Maxwell stress (αdisp,z ∼ αMRI ). A larger a does not change
this conclusion, because we then would have a higher αMRI or the e-heating zone would
vanish.

3.6. DISCUSSION

3.6.2

61

Eﬀects of Grain Size Distribution and Porosity

We have characterized dust grains with a single particle size a assuming that the size
distribution of dust grains is narrow. Under this assumption, the e-heating zone covers
only a small part of protoplanetary disks when the particles grow to millimeter sizes(see Figure 3.6). However, caution is required in applying our results to more generalcases where particles have a size distribution. In such cases, the smallest grains tend to
dominate the total surface area of dust (which controls the ionization balance), whereas
the largest grains tend to dominate the total mass of dust, simply because smaller grains
have a larger area-to-mass ratio. Therefore, it is not obvious what the typical particle
size is in these cases.
Here we discuss more quantitatively how we can apply the results of single-size calculations to cases with a size distribution. Let us assume that the particle size distribution
is given by the power-law form
dnd
3ρfdg
=
a−3.5
√
da
8πρ• amax(3.52)

with amin < a < amax (amin ≪ amax ), where dnd /da is the number density of dust
particles per unit particle radius, and amin and amax are the minimum and maximum
particle sizes, respectively. The distribution is normalized so that that the total particle∫
mass density md (dnd /da)da becomes equal to ρfdg . Equation (3.52) applies when theparticle size distribution is determined by fragmentation cascade (Dohnanyi, 1969) andis also known to reproduce the size distribution of interstellar dust grains (Mathis et al.,1977). The quantity we are interested in is the total surface area of the particles as it
mainly determines the ionization balance in a gas–dust mixture (e.g., Sano et al., 2000).
This can be calculated as∫ amax
4πa2
amin

dnd
3ρfdg
1
da ≈
.
√
da
ρ•
amin amax(3.53)√
Note that the factor 1/ amin comes from the fact that the integration in Equation (3.53)
is dominated by the smallest particles (because a2 (dnd /da)da ∝ d(a−0.5 )), whereas the√factor 1/ amax from the fact that the total mass is dominated by the largest particles.

62

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

By contrast, if all dust particles have a single size asingle , their total surface area is
4πa2single nd,single = 3ρfdg /(ρ• asingle ). Comparing this with Equation (3.53), we find thatthe total surface area of particles whose size distribution is given by Equation (3.52) isequal to that of single-size particles if
asingle =√amin amax .

(3.54)

Since the total surface area approximately determines the ionization state, Equation (3.54)
may be used to generalize the results presented in this study to the cases where theparticle size distribution obeys Equation (3.52).
Observations of millimeter dust emission from protoplanetary disks suggest that thelargest dust particles in the disk have a size of centimeters (e.g., Testi et al., 2003;
Natta et al., 2004; Rodmann et al., 2006; Ricci et al., 2010). Assuming amax = 1 cm
and amin = 0.1 µm, we obtain asingle = 30 µm. In this case, we expect from Table 3.1
that the e-heating zone extends to ∼ 15 AU. Thus, even if cm-sized grains exist in
protoplanetary disks and the total mass of grains is dominated by such large grains,the e-heating zone can be present in the disks.
For the same reason, large dust particles can alone provide a large e-heating zone if thedust particles are highly fluﬀy aggregates of tiny grains. Okuzumi (2009) showed thatthe ionization balance is insensitive to the particle radius when the fractal dimensionis ≈ 2, for which the total surface area of the aggregates is approximately conserved
during the aggregation process.

3.6.3

Hall Eﬀect and Ambipolar Diﬀusion

The plasma heating model employed in this study neglects the eﬀects of magnetic
fields on the motion of plasma particles. In terms of non-ideal magnetohydrodynamics,
this is equivalent to neglecting ambipolar diﬀusion and Hall eﬀect (see, e.g., Wardle,
1999). A full treatment of these non-Ohmic eﬀects introduces to the model additional
complexities arising from the relative angle between the magnetic and electric fields(Okuzumi, Mori, & Inutuska, in prep.), which is beyond the scope of this work. In

3.6. DISCUSSION

63

Table 3.2: Am and Ion Abundance xi in E-heating Zone for Various Parameter Setsβc

a (µm)

fdg

102
103
104
105
103
103
103
103
103
103
103
103
103
103
103
103

0.1
0.1
0.1
0.1
0.1
1
10
100
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1

10−2
10−2
10−2
10−2
10−2
10−2
10−2
10−2
10−1
10−2
10−3
10−4
10−2
10−2
10−2
10−2

fΣ

Am in e-heating zone
Inner edge Outer edge
1
0.14
0.56
1
0.17
0.62
1
0.23
0.62
1
0.41
0.62
1
0.17
0.62
1
0.21
0.72
1
0.43
0.82
1
0.54
0.72
1
0.16
0.57
1
0.17
0.62
1
0.24
0.71
1
0.41
0.80
10
0.32
1.26
3
0.22
0.82
1
0.17
0.62
0.3
0.16
0.65

xi in e-heating zone
Inner edge Outer edge
1.9×10−12 4.7×10−11
3.4×10−12 5.9×10−11
7.4×10−12 5.9×10−11
2.5×10−11 5.9×10−11
3.4×10−12 5.9×10−11
1.7×10−12 2.8×10−11
2.3×10−12 1.3×10−11
2.6×10−12 5.6×10−12
8.5×10−12 1.2×10−10
3.4×10−12 5.9×10−11
2.1×10−12 2.8×10−11
2.3×10−12 1.3×10−11
1.8×10−12 2.5×10−11
2.5×10−12 3.9×10−11
3.4×10−12 5.9×10−11
5.4×10−12 9.5×10−11

this subsection, we only briefly discuss how plasma heating and these non-ideal MHD
eﬀects could aﬀect each other.
Ambipolar diﬀusion can suppress MRI in low density regions of protoplanetary disks(e.g., Blaes & Balbus, 1994; Hawley & Stone, 1998; Kunz & Balbus, 2004; Desch, 2004;
Bai & Stone, 2011; Simon et al., 2013a,b). If MRI is eﬀectively suppressed in the eheating zone, electric fields may not suﬃciently grow to cause electron heating. The
eﬀectiveness of ambipolar diﬀusion is characterized by the ambipolar Elsasser number
Am = γi ρi /Ω (e.g., Blaes & Balbus, 1994; Lesur et al., 2014), where γi = ⟨σin vin ⟩ /(mn +mi ) and ρi = mi ni . According to MHD simulations including ambipolar diﬀusion,
MRI-driven turbulence behaves as in the ideal MHD limit if Am ≫ 1, while ambipolar
diﬀusion suppresses turbulence if Am ≪ 1 (e.g., Bai & Stone, 2011). Table 3.2 lists
the values of Am as well as the ion abundance xi = ni /nn at the inner and outer edgesof the e-heating zone before electron heating sets in (E = 0). We find that Am ≈0.2–0.7, implying that ambipolar diﬀusion would moderately aﬀect MRI turbulencein the e-heating zone. Therefore, MHD simulations including both electron heating
and ambipolar diﬀusion are needed to assess which eﬀect determines the saturation
amplitude of MRI turbulence in these outer regions of the disks.

64

CHAPTER 3. ELECTRON HEATING ZONE : MORI & OKUZUMI (2016)

The Hall eﬀect is also important at r ∼ 10–50 AU (see Figure 1 of Turner et al., 2014).
The Hall eﬀect can either damp or amplify magnetic fields, which depends on therelative orientation between the disk’s magnetic field and rotation axis and on the signof the Hall conductivity (e.g., Bai, 2014; Wardle & Salmeron, 2012). At relatively highgas densities (nn ≳ 1010 cm−3 ), the Hall conductivity is usually positive (Wardle & Ng,
1999; Nakano et al., 2002; Salmeron & Wardle, 2003), but can become negative when
the number density of electrons is significantly lower than that of ions. Interestingly, ourpreliminary investigation shows that the Hall conductivity can indeed become negative
as the electron number density is decreased by electron heating (Okuzumi et al., in
prep.). This suggests that electron heating might reverse the role of the Hall term.
Whether this occurs under conditions relevant to protoplanetary disks will be studiedin future work.

Chapter 4
Suppression of MRI by the Electron
Heating
In this chapter, we perform MHD simulations to investigate an eﬀect of the electron
heating on the nonlinear evolution of MRI. Our previous study (Mori & Okuzumi, 2016,or see previous chapter) have investigated the extent where the electron heating takes
place. They found that the electron heating is eﬀective in a large region of protoplanetary disks. Electron heating evoked by a strong electric field of MRI turbulence reduces
the ionization fraction. When ionization fraction is suﬃciently low, MRI is stabilized.
Thus, the electron heating has potential to suppress the MRI. We have also estimatedhow electron heating can weaken the MRI turbulence. Then, we found the turbulence
can be largely suppressed. However, the estimation is based on only a scaling relation
between magnetic turbulence stress and current density. Thus, the detailed eﬀects ofelectron heating on the turbulence strength is still not obvious. In addition, the nonlinear evolution of MRI is also unknown. This paper first performs MHD simulations
including electron heating. Our goal in this work is to quantify the eﬀect of electron
heating on MRI. We investigate whether the scaling relation exists or not and if so,
how much the electron heating suppresses the turbulence strength depending on current density. In this chapter, for simple analysis of the increasing resistivity by electron
heating, we neglect the nonideal MHD eﬀects other than Ohmic dissipation.
First, in Section 4.1, we present the numerical setup and procedure in our simulations.
65

66

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATINGIn Section 4.2, we then show some results and present the interpretations. Finally,in Section 4.3, we analytically derive a relation between current density and Maxwell
stress.

4.1

Numerical Setup and Procedure

We adopt a local shearing box whose center is located at a distance from the central starr0 . The coordinates x, y, and z are radial, azimuthal, and vertical axis, respectively.
We numerically solve the following MHD equations (the equation of continuity, theequation of motion, and Faraday’s equation)∂ρ+ ∇ · (ρv) = 0,
(4.1)
∂t
(
)
∂v
1
B2
1
+ (v · ∇) v = −2Ω × v + 3Ω2 xex − ∇ P ++(B · ∇B) , (4.2)
∂t
ρ
8π
4πρ
∂B
= −c∇ × E,
(4.3)
∂twhere v is the gas velocity, ρ is the mass density, P is the gas pressure, Ω is the angularvelocity at r0 , B is the magnetic field, E is the electric field, and c is the speed oflight. In this chapter, we assume the isothermal change of gas and use the isothermal
equation of state,
P = c2s ρ.

(4.4)

Electric field E is given byE=

4πη(E ′ )1
J − v × B.
2
c
c

(4.5)

Here, we have used the Ohm’s law in comoving frame,
c2
E ′,
J =
′
4πη(E )′

(4.6)

4.1. NUMERICAL SETUP AND PROCEDURE

67

and the Lorentz transformation of current density J,
J′ = J,(4.7)

where the ′ means values in the comoving frame. In this chapter, since we use ourresistivity model instead of usual constant resistivity, the dependence of the resistivity
on electric field strength is explicitly expressed. Adopting the same formula of MO16,
we use the electric field in the comoving frame. In our simulations, for electric field, we
express it only in a comoving frame in order to simply derive electric fields, E ′ = ηJ.
In Section 4.2, note that we omit the superscript ′ of electric field strength.
The relation between current density and electric field including electron heating (nonlinear Ohm’s law) has N-shape (see Figure 10 of OI15; see also Figure 3.4 ). In MHD
formulation, electric field is obtained from the current density via the Ohm’s law. However, under the nonlinear Ohm’s law, we obtain multiple values of electric fields whena current density is given. Thus, in the MHD simulation, we cannot exactly treat thenonlinear Ohm’s law. In this work, we focus on only the resistivity increased by electron
heating but does not address an eﬀect of the complex N-shaped Ohm’s law. An exact
treatment of Ohm’s law is beyond the scope of this work.
To include the increase of resistivity by electron heating, we introduce a analyticalresistivity model depending on E. From above the reason, mimicking the nonlinear
Ohm’s law, we give the simple analytical resistivity model where the resistivity increases
but does not drops current density with increase of E during electron heating. Figure
4.1 shows the schematic diagram. Specific formula of the resistivity is expressed as






or

c2 ′
E,
E ′ < EEH ,
4πη
0
J′ =
( ′ )γ−1
2

c
E


E ′ , E ′ > EEH ,
4πη0 EEH

(4.8)

 4πη
0 ′

J ′ < JEH , c2 J ,( ′ )1/γ −1
E′ =

4πη0
J


J ′ , J ′ > JEH .
2
c
JEH

(4.9)

where γ is a constant value suﬃciently less than unity, EEH is an electric fields strengthCHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

Current Density!

68

Nonlinear Ohm’s Law!
Model!

Slope γ!

JEH!

Electric ﬁeld Strength!

EEH!

Figure 10 in OI15!

Our Model!

Figure 4.1: Schematic diagram of modeling our J–E relation after the nonlinear
Ohm’s law. We neglect the N-shape of the nonlinear Ohm’s law. Instead, we givethe shallow inclination of J–E relation during electron heating. In this model, theresistivity increases with increase of E.

when electron heating start to take place, and EEH =as

4πη0
JEH . Especially, η is written
c2

 4πη
0

J ′ < JEH , c2 ,( ′ )1/γ −1
η=

4πη0
J


, J ′ > JEH .
2c
JEH(4.10)In this chapter, we take γ to be 0.1 and JEH to be arbitrary parameter. Using thismodel, we solve the MHD equation.
We use a shearing box and thus set a uniform shear flow, (vx0 , vy0 , vz0 ) = (0, − 23 Ωx, 0).
The simulation box size is (Lx , Ly , Lz ) = (H, 2πH, H), where Lx , Ly , and Lz are
lengths of simulation box in each directions and H is the gas scale height. We use theperiodic boundary conditions at each boundary x = −0.5H, 0.5H; y = −πH, πH; z =−0.5H, 0.5H because we assume the simulated phenomena occurs in local region.
Following Hawley et al. (1995), we define the computational units of length Lu , time Tu ,and mass Mu as H, 10−3 Ω−1 , and ρ0 H 3 , respectively, where the subscript “u” means
units of our simulation results. Therefore, the unit of velocity is given by Vu = 103 cs ,the unit of density ρu = ρ0 , and the unit of pressure Pu = 106 ρ0 cs . The sound speed,
the initial density, and the initial pressure of numerical values (cs /Vu , ρ0 /ρu , and P0 /Pu )are, respectively, given by cs /Vu = 10−3 , ρ0 /ρu = 1, and P0 /Pu = 10−6 . Let the unit of4.1. NUMERICAL SETUP AND PROCEDURE
magnetic field strength be
Bu =√8πP0 .

69(4.11)

We obtain the unit of current density from Ampére’s equation,
c Bu
.
4π H

(4.12)

E
η J
=
,
Eu
ηu Ju

(4.13)

Ju =To write the Ohm’s law E = (4πη/c2 )J as

we define the unit of electric fields strength as
Eu =

Bu ηu
cs
= 103 Bu ,cHc(4.14)

where we use ηu = L2u /Tu = 103 Hcs . Using these units, in this work, we define Jmax
as
Jmax = 10Ju .

(4.15)

Here, Jmax indicates just a typical value of the resulting maximum value of current density in fully developed MRI (Muranushi et al., 2012). Note that the resulting maximum
current does not necessarily equal to 10Ju .

We set initial setup as Figure 4.2. We give the uniform vertical magnetic fields whosestrength is−1/2

Bz0 = β0Bu ,(4.16)

where the initial plasma beta vales β0 is
β0 =

8πP0
.
2
Bz0

(4.17)In order for MRI to grow in the beginning, we give the initial Elsasser number Λ0 as
Λ0 = 10. Here, in this work, the Elsasser number is defined by
Λ=

vA2
.
ηΩ

(4.18)

70

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

J!
10Ju!

Jmax!

JEH!

EEH!

E!

Figure 4.2: Schematic diagram of initial setup of our calculation. Initially, for MRI
to grow, we give the Elsasser number Λ0 = 10. Then, for the current growth to be
suppressed, we give JEH , at which electron heating begin to occur, to be less than Jmax
Because an eﬃciency of ohmic dissipation is changed at Λ = 1, we can approximately
expect MRI growth in Λ > 1 and MRI decay in Λ < 1. Thus we can also expectthat final saturated current density converges on Λ = 1. In order to investigate thedependence of the results of a saturated current density, we changes the parameter
JEH , which expresses where in a J–E plane, electron heating starts up. We randomly
gives the initial amplitude of pressure perturbations as δP/P0 = 2.5 × 10−4 and theone of velocity perturbations as |δv|/cs = 5 × 10−5 . The initial perturbation values are
two orders of magnitude less than the ones of (Hawley et al., 1995). Because we treat
very small fluctuations in our simulation, to properly treat them, we note that initial
perturbations should be less than amplified fluctuations.
In this study, we use Athena, open source code which uses Godunov’s scheme. We perform MHD simulations with a unstratified local sharing box. The fluid is compressible,
isothermal, and magnetized.
In this work, we plot all values as the volume-averaged quantities as ⟨...⟩ and express also
the time-averaged quantities as ⟨⟨...⟩⟩. For example, the volume average of a variable ais written as∫
a dV⟨a⟩ = ∫,
dV

(4.19)

4.1. NUMERICAL SETUP AND PROCEDURE

71

where the integration region is the simulation box −0.5H < x < 0.5H, −πH < y < πH,and −0.5H < z < 0.5H. The time average of a is written as∫a dt⟨⟨a⟩⟩ = ∫dt

(4.20)

where we average a value from 20 orbits to 30 orbits in this work. Here, we express
time with “orbits” which indicates one orbital time, 2π/Ω. The integration interval
corresponds to final state in our simulation, and thereby the time averaged value is
expected to means the saturated value.
The Courant-Friedrichs-Lewy number, which is a ratio of fluid velocity to computablevelocity, is fixed at 0.4. In order for the number to be fixed, time intervals per one time
step is variable since the size of a mesh is constant. We take the number of meshes tobe (Nx , Ny , Nz ) = (64, 128, 64) in each direction. For MRI in a uniform vertical field,
we have to note that at least vertical size of one cell (i.e. resolution) should be less thanthe most unstable wavelength to properly treat MRI growth. The most unstable modeis written asvA,zλm.u., ideal = 2π
= 2π
Ω

√

Therefore, because Nz requires Nz > 11.2

2
H = 0.089
β0( β0 )1/2
104

(

β0
104

)−1/2
H.

(4.21), Nz = 64 is reasonable resolution in

our calculation.
The most interesting values in our calculation is the Maxwells stress, which greatly
aﬀects the disk evolution. Moreover, if we simply assume the magnetic energy to be thekinetic energy, the Maxwell stress also controls the velocity fluctuation which aﬀects
dynamics of solid particles. The Maxwell stress ⟨−Bx By ⟩ /(4π) is often expressed as a
Shakura-Sunyaev α parameter which is defined as the accretion stress normalized thepressure. When the Maxwell stress dominates in the accretion stress, α is written asx By ⟩⟩
α = ⟨⟨−B. Although pressure P varies, in this work, we use a dimensionless value
4πP

which is defined as the Maxwell stress normalized by the initial pressure,αM =⟨⟨−Bx By ⟩⟩.
4πP0(4.22)

72

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING
Table 4.1: Summary of results.

Run ID
34
35
36
37
38
39
40
41
42

β0
104
104
104
104
104
104
104
104
104

JEH
0.001
0.003
0.01
0.03
0.1
0.3
1.0
3.0
∞

Λ0
10
10
10
10
10
10
10
10
10

⟨⟨

⟩⟩
B 2 /(8πP0 )(1.0 ± 0.0) × 10−4
(1.0 ± 0.0) × 10−4
(1.1 ± 0.0) × 10−4
(1.6 ± 0.0) × 10−4
(7.6 ± 0.0) × 10−4
(1.8 ± 1.2) × 10−3
(7.0 ± 2.4) × 10−3
(1.3 ± 0.1) × 10−2
(6.2 ± 2.1) × 10−2⟨⟨αM ⟩⟩(2.4 ± 0.3) × 10−8
(2.6 ± 0.0) × 10−7
(2.8 ± 0.0) × 10−6
(2.6 ± 0.0) × 10−5
(2.8 ± 0.0) × 10−4
(7.6 ± 2.1) × 10−4
(4.2 ± 0.7) × 10−3
(1.1 ± 0.0) × 10−2
(5.6 ± 0.7) × 10−2⟨⟨ρvx δvy ⟩⟩ /(8πP0 )(2.0 ± 0.4) × 10−10(1.5 ± 0.2) × 10−9(1.5 ± 0.2) × 10−8
(1.5 ± 0.1) × 10−7
(1.7 ± 0.2) × 10−6
(9.9 ± 9.7) × 10−4
(2.7 ± 2.7) × 10−3
(2.9 ± 1.5) × 10−3
(8.3 ± 2.9) × 10−3⟨⟨J⟩⟩ /Ju
(1.7 ± 0.1) × 10−3
(5.2 ± 0.0) × 10−3
(1.7 ± 0.0) × 10−2
(5.2 ± 0.0) × 10−2
(1.7 ± 0.0) × 10−1
(4.1 ± 0.3) × 10−1
(1.2 ± 0.0)
(2.8 ± 0.0)
(8.5 ± 0.9)

We have confirmed that the pressure P can vary within a factor of two at most in our
calculations. In addition, we define a dimensionless value with respect to resistivity
η,

( 2 )
vA,z0
Λη =
ηΩ

(4.23)

where vA,z0 is the Alfvén velocity consisting of initial magnetic field strength, vA,z0 =√Bz0 / 4πρ0 . Λη expresses the Elsasser number assuming that magnetic field strengthand density do not vary. In Figure 3.4 of Section 3.2, we plot a line of Λη = 1 in a J–E
plane to show the threshold of the resistive MHD. The magnetic field is not amplified
during a linear growth, and thus Λη exactly indicates the Elsasser number Λ.
In our previous work (MO16), we have found out the scaling relation between theMaxwell stress (or αM ) and the current density. In the next section, in order to verify
the scaling relation, we plot the Maxwell stress and the current density.

4.2

Results of Shearing-Box MHD Simulations

In this section, we present the representative simulation results. We set the initial
parameters as β0 = 104 and Λ0 = 10. We give JEH of 8 parameters as JEH /Ju =3 × 10−4 , 1 × 10−2 , 3 × 10−2 , 1 × 10−1 , 3 × 10−1 , 1, 3, and ∞, where JEH /Ju = ∞ means
the ideal MHD case. The initial setup and result are summarized in Table 4.1.
Figure 4.3 shows the finale state (t = 30 orbits) of the representative our simulation
results of JEH /Ju = 0.03 and the ideal MHD case, JEH /Ju = ∞, for comparison. We
also show the x–z and y–z plane at the final state of the each parameters in Figure 4.4.

4.2. RESULTS OF SHEARING-BOX MHD SIMULATIONSOur J-E Relation

73

Ideal MHD

Figure 4.3: Outline of representative our MHD shearing box simulations. The left
figure shows the distribution of magnetic field strength |B| with our J–E relation of
JEH /Ju = 0.03. The right figure shows the result in ideal MHD case for comparison.
These figures show the final state of the calculations. Note that the color range of thetwo is largely deferent.
We find the ordered magnetic field and the non-turbulent steady state, the laminar flowor the stationary state, with our J–E relation model, whereas the magnetic turbulenceis generated in the ideal MHD case because of the MRI nonlinear growth. In the nonturbulent state, the magnetic field is stretched in the y-direction, and thereby sinusoidal
structure of the magnetic field forms in y–z plane (the upper left panel in the upper
block in Figure 4.4). Because azimuthal magnetic field By varies in the z-direction,
the layered contour of magnetic field strength appears. Such an anti-parallel magnetic
field structure forms a current sheet which is maximum in calculating box. In ourmodel, the growth of the current sheet is suppressed by the increasing resistivity. Thus,
the initial uniform magnetic field is stretched but its growth stops at the sinusoidal
structure.
In Figure 4.5, we shows the tracks of the volume-averaged current density ⟨J⟩ andelectric field strength ⟨E⟩ in order to check our model. The current densities initiallygrow on the line of Λ0 = 10 and then branch oﬀ the line after the current densities
reaches each JEH . We confirm that our model works according to our expectation of J–
E relation. We also confirm that the Elsasser number of initial magnetic field strength

74

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

Figure 4.4: Magnetic field distribution in the x–z plane (−0.5H ≤ x ≤ 0.5H, −0.5H ≤z ≤ 0.5H; upper block including four panels) and y–z plane (−0.5H ≤ y ≤0.5H, −0.5H ≤ z ≤ 0.5H; lower block including four panels) at the final state (t = 30
orbits) of the each calculation. The each panel in blocks shows the case of JEH /Ju = 0.03
(upper left), 0.3 (upper right), 3 (lower left), and the ideal MHD (lower right). The
arrows show directions of the magnetic field in x-z plane, but does not express themagnitude. The colors show the magnetic field strength |B| , but note that the color
scale varies at each initial parameters.

4.2. RESULTS OF SHEARING-BOX MHD SIMULATIONS
102
1

10

0

〈J〉 [Ju]

10

10

75
Jmax
Λη = 10Λη = 1Λη = 0.1
34
35
36
37
38
39
40
41
42

-1

10-2
10-3
10

-4

10-5

10-11 10-10 10-9 10-8 10-7 10-6 10-5
〈E〉 [Eu]

Figure 4.5: Tracks in the current-field plane. The vertical axis shows volume-averaged
current density ⟨J⟩, and the horizontal axis volume-averaged electric field strength ⟨E⟩.
The colors show run ID which means variation of JEH . The dashed lines of white grey,
grey, and dark grey show the initial Elsasser number Λη = 10, 1, and 0.1. The circles
are plotted at equal intervals (∆t = 1 orbits) to express the speed in the J–E plane.
The white filled circles show that saturated vale of ⟨J⟩ and ⟨E⟩ of each JEH .

Λη is approximately saturated at Λη = 0.1.

In Figure 4.6, we shows the time variation of the Maxwell stress varying the parameter
JEH which control saturated current density. We found that the saturated Maxwell
stress depends on JEH . Thus saturated Maxwell stress depends on the saturated current density. Lower current density provides lower Maxwell stress. Therefore, we canconclude that MRI growth stops if electron heating begin to take place at much low
electric field strength. We also find the saturated MRI behavior can be divided into thetwo case, the fluctuated case (JEH ≲ 0.1) and the stationary case (JEH ≳ 0.3). According to Figure 4.3, the two case would correspond to the turbulent state and laminarstate. Thus, we can find the threshold of current density at which the MRI completely
dies away. Hereafter, we call the threshold by Jturb /Ju ≈ 0.3.

In MO16, we used a scaling relation between the Maxwell stress and current density toinvestigate the possibility of suppressing MRI by electron heating. The scaling relation

76

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

Maxwell Stress [P0]

10-1
10

-2

10

-3

10

-4

34
35
36
37
38
39
40
41
42

10-5
10-6
10-7
10-8
10-9

10-10

0

5

10

15
20
Time [Orbits]

25

30

Figure 4.6: Time variation of αM with variating JEH . The colors show the diﬀerent
value of JEH , and the black solid line shows the ideal MHD case.

100

34
35
36
37
38
39
40
41
42

Maxwell Stress [P0]

10-2
10-4
10-6
10-8

10-10
10-12

0

5

10

15
20
Time [Orbits]

25

30

Figure 4.7: Figure 4.6 with plot of the analytical scaling relation presented in MO16.
The dashed lines show the scaling relation.

4.2. RESULTS OF SHEARING-BOX MHD SIMULATIONS
100

10

-2

〈〈-Bx By〉〉/4π [P0]

10

-1

10

2

0.005 (〈〈J〉〉/Ju)

-3

10-4

77
34
35
36
37
38
39
40
41
42

10-5
10-6
10

-7

10

-8

10-9 -3
10

10-2

10-1
〈〈J〉〉 [Ju]

100

101

Figure 4.8: Dependence of volume-averaged ⟨αM ⟩ to volume-averaged ⟨J⟩ (coloreddots). The colors show the same as Figure 4.3. The dashed line shows the fitting
formulae in low current region, ⟨αM ⟩ = 0.005(⟨J⟩ /Ju )2 or 0.5(⟨J⟩ /Jmax )2 .
is written as(
αM ≈ 0.02

β0
104

)−1 (

J
Jmax

)2
.

(4.24)

However, the equation is supported only by the scaling law. In addition, even theexistence of this relation is not clear because many previous studies have not focused
the relation of current density to the Maxwell stress. Thus, in order to verify thisrelation, we plot the formula over the resulting Maxwell stress in Figure 4.7. The
relation is expressed as a function only of the initial plasma beta β0 and current density
J. Thus, if the relation is correct, the relation should be similar to the result at eachtime. However, Figure 4.7 shows that the relation does not follow the resulting Maxwell
stress, although the relation follows the plot during the MRI linear growth. At the final
state, the saturated Maxwell stress of the two is diﬀerent by a factor of 3–10. In Section
4.3, we will investigate the cause of this discrepancy and will propose a new predictive
formula for α as a function of J. We improve the relation between the Maxwell stressand current density (Equation (4.24)) using a more sophisticated scaling law (Equation(4.34) in Section 4.3).
To obtain an empirical formula of the stress–current relation in saturated state, taking

78

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

〈〈Λ〉〉

104
10

3

10

2

10

1

34
35
36
37
38
39
40
41
42

100
10 (〈〈J〉〉/Ju)2

10-1
10-2 -3
10

10-2

10-1
〈〈J〉〉 [Ju]

100

101

Figure 4.9: Dependence of time- and volume-averaged Elsasser number ⟨⟨Λ⟩⟩ on timeand volume-averaged ⟨⟨J⟩⟩.
the time average of the Maxwell stress and current density, we show dependence of thesaturated Maxwell stress ⟨⟨αM ⟩⟩ on saturated current density ⟨⟨J⟩⟩ in Figure 4.8. In
this chapter, since the interval of the time integration is from 20 orbits to 30 orbits,
the time average means saturated values. We confirm an obvious correlation between
the two quantities. Fitting the correlation with a quadratic function, we obtain theempirical formula of the relation, 0.005(⟨⟨J⟩⟩ /Ju )2 or 0.5(⟨⟨J⟩⟩ /Jmax )2 . We can seethat the Maxwell stress largely depends on current density. This formula means that alower current leads to lower turbulence strength.. The dependence on current density,
αM ∝ J 2 , is consistent with Equation (4.24) presented in MO16. Thus, in Equation 4.24,
the dependence of current density is consistent with the results, whereas the magnitude
is diﬀerent from the results. Note that since in this calculation the initial plasma beta
β0 is constant, the dependence on β0 does not appear in this fitting formula.
Here we address what determines the boundary of the laminar and turbulent state,
namely what determines Jturb . In Figure 4.8, we can see that the fitting formula isslightly incorrect at high current density (⟨J⟩ ≳ Jturb ),. This may concern the nonlinear
evolution of the Elsasser number. Figure 4.9 shows dependence of the volume- andtime-averaged Elsasser number ⟨⟨Λ⟩⟩ on current density. Here, the Elsasser number isdefined as Λ = vA2 /(ηΩ) = |B|2 /(4πρ0 ηΩ). As we see, at J > Jturb , we can see that

4.2. RESULTS OF SHEARING-BOX MHD SIMULATIONS
5

34
35
36
37
38
39
40
41
42

-1

200 (〈〈J〉〉/Ju)

〈〈β〉〉

104

103

10

-5

6e-7 (〈〈J〉〉/Ju)-110-6〈〈η〉〉 [ηu]

10

79

10

34
35
36
37
38
39
40
41
42

-7

102

101 -3
10

10-2

10-1
〈〈J〉〉 [Ju]

100

101

10-8 -3
10

10-2

10-1
〈〈J〉〉 [Ju]

100

101

Figure 4.10: Dependence of time- and volume-averaged plasma beta ⟨⟨β⟩⟩ and resistivity⟨⟨η⟩⟩ on time- and volume-averaged ⟨⟨J⟩⟩. The colors shows the same as Figure 4.9.
the threshold current density Jturb approximately corresponds to the current density
at Λ = 1. Thus, we can naturally understand that the MRI does not generate themagnetic turbulence because MHD is resistive at Λ ≲ 1. On the other hands, the MRI
generates the turbulence at Λ ≳ 1.
The increase of the Elsasser number is understood by the following. The Elsasser
number can be written asvA2
B 2 ρc2s1 ηuΛ=== 2 × 10−3.
ηΩ
4πP ηΩβ ηThus, Λ is proportional to β −1 η −1 . Figures 4.10 shows dependence of the volumeand time-averaged plasma beta ⟨⟨β⟩⟩ (left panel) and resistivity ⟨⟨η⟩⟩ on ⟨⟨J⟩⟩(rightpanel). We can see that when ⟨⟨Λ⟩⟩ begins to increase, ⟨⟨Λ⟩⟩ also begins to increase.
⟨⟨η⟩⟩ begins to decrease when ⟨⟨Λ⟩⟩ reaches unity. Thus, we can expect that ⟨⟨Λ⟩⟩ isincreased by ⟨⟨β⟩⟩ in 0.03Ju < J < Jturb . ⟨⟨η⟩⟩ decreases due to changing state of flow
and thereby ⟨⟨Λ⟩⟩ increases with ⟨⟨J 2 ⟩⟩.
To see why the MRI is saturated in steady non-turbulent state when the resistivity
by the electron heating, in Figure 4.11, we shows the time variation of the critical
wavelength. The critical wavelength is the shortest wavelength in unstable MRI modeand defined as)
1 vA,z (λcrit,z ≈ 2π √1 + Λ−1,z
3 Ω(4.25)

which is given by the dispersion relation in a steady state (Sano & Miyama, 1999). The

80

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

Critical Wavelength/H

101

10

1
34
35
36
37
38
39
40
41
42

0

10-1

10-2

0

5

10

15
20
Time [Orbits]

25

30

Figure 4.11: Critical wavelength λcrit,z (lower panel) in vertical direction and normalizedby H. The colors show the same as Figure 4.3.

critical wavelength is approximately equal to simulation box size H. This indicatesthat unstable wavelength become larger and larger during MRI growth, and then all
MRI unstable mode dies when the unstable wavelength cannot exist. Thus, in thissimulation, the final state is determined by the box size.

4.3

Predictive Formulae

In this section, we derive a formula which predicts the Maxwell stress from current
density. Simulation of MRI can treat significantly shorter time than the disk lifetime.
Thus, such the predictive formula enable us to introduce the eﬀect of turbulence as αparameter. In Subsection 4.3.1, we derive a fitting formula which predicts the Maxwell
stress in both the laminar and turbulent state. Then, in Subsection 4.3.2, focusing thatthe fluctuation can be assumed to be first order, we derive an analytical expressionof the laminar state that exactly gives the physical quantities. The expression also
provides an interpretation of the scaling relation in the laminar state.

4.3. PREDICTIVE FORMULAE

4.3.1

81

Scaling Relation

Here we improve the scaling relation presented in MO16 and find the practical predictive
formula. First, we determine an fitting formula of the laminar state from a scaling law,and then we determine that in the ideal MHD case. After that, we link the two formulae
together by taking the average to reproduce simulation results.

We first derive a scaling relation between the Maxwell stress and the current density,
noting that MHD is resistive when MRI is saturated by the increased resistivity. To
relate current density to Maxwell stress, we take ∇ to be the typical wavenumber k inc∇×B. Here, we assume that k is the critical wavenumber
the Ampére’s equation J = 4π

in the vertical direction,k = kz,crit =√ vAz03
η(4.26)

The current densities of x and y components are described asckz,crit By ,
4π
c
Jy ≈
kz,crit Bx .
4π

Jx ≈ −(4.27)(4.28)

Therefore, using the current density, Maxwell stress ⟨Bx By ⟩ /(4π) normalized P0 canbe expressed as
αM

(

)2
4π
c
)
(
)2 (
1
η
Jx Jy
= 200
,
2
3 vAz0 H
Jmax)( 2 )()2 (ηΩ
1 vAz0
Jx Jy
,
= 200
2
2
3
c2s
vAz0
Jmax
)
(
)
(
1
Jx Jy
1
−2
= 400
Λη
,
2
3 βz0
Jmax
1
=
4πP0

Jx Jy
2
kz,crit)((4.29)

where the current densities is normalized by the maximum current density in fully√ρdeveloped MRI, Jmax = 10 2πcΩ,(Muranushi et al., 2012). J would be dominated by
Jx because By is generally larger than Bx (e.g., Hawley et al., 1995; Sano et al., 2004),

82

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

and therefore we assume that Jx ≈ J and Jy ≈ J(Bx /By ).
1
αM = 400
3

(

1
βz0

)

2
−2 Bx JΛη,
2
By Jmax(4.30)

As we see below, we here assume the ratio of Bx to By in the saturated laminar state(Equation (4.41)),Bx= −2Λη,sat = −0.1
By

(

βz0
104

)−1/2
.

(4.31)

Then, we obtain(
αM,lami = 1.3 × 10

−3

βz0
104

)−1.5

(
Λ−2
η

)2

J
Jmax

.

(4.32)

This expression provides the Maxwell stress once current density and Λη is given.
On the other hand, in the case of the ideal MHD, the scaling relation has been obtainedin MO16 (or see Equation (4.24))(−2

αM,turb ≈ 3 × 10

βz0
104

)−1 (

J

)2

Jmax

,

(4.33)

where we have modified a factor of 1.5 to fit with our results.
Finally, we take link the two formulae together by taking the average. For best fitting,
we square a sum of root values,()2(αM, turb )1/2 + (αM, lami )1/2
(
)1/2 2 (()−1 )1/2 (
(
)−1.5
)2
βz0
βz0
J
−2
−3
−2


=
3 × 10+ 1.3 × 10
Λη
(4.34)
104
104
Jmax

αM =

In Figure 4.12, we plot (4.34) over Figure 4.6. We can see that the equation surprisingly
reproduce the resulting Maxwell stress at each time. The formula requires the initialplasma beta β0 , current density J, and Λη . β0 is usually given by initial magnetic field.
Λη is analytically given by (4.38) in the laminar state. Even in the turbulent state, Λη
can be estimated to be approximately 0.1. In the J–E plane, Λη corresponds to thecontour whose slope is proportional to J. Thus, once Λη is given, J also is determined.

4.3. PREDICTIVE FORMULAE
10

83

0

34
35
36
37
38
39
40
41
42

Maxwell Stress [P0]

10-2
10-4
10-6
10-8

10-10
10-12

0

5

10

15
20
Time [Orbits]

25

30

Figure 4.12: Figure 4.6 with plot of Equation (4.34). The dashed lines show the equation.
Using this formula, we can predict the Maxwell stress from the J–E relation at alocation in protoplanetary disks.
We have to note that these results is based just on the simple analytic J–E relation. In
reality, the electron heating rises the N-shaped nonlinear Ohm’s law. Of course, finally,
the real eﬀect of electron heating on MRI is investigated by considering the realistic
resistivity model with a MHD simulation. However, this work first unveils the relationof Maxwell stress and current density and the state when MRI completely dies away bysuppressing current density. The importance of current density in MRI is stressed. We
eventually conclude that the reliability of the suppression of MRI by electron heating
has been shown by this work.

4.3.2

Analytical Solution of Laminar State

Next, we present a solution of the laminar state that provides a better understandingof the empirical formula of Equation (4.32). As we see below, the laminar state can beobtained by an analytical manner. A point of the manner is to assume that linearized
equations are available because the fluctuations of the final state is suﬃciently small
to be first order displacement even in the “nonlinear” stage. Actually, the nonlinearityof the simulation is imposed on the evolution of resistivity. In addition, some features

84

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATINGof the laminar state is also available. At the laminar state, one wavelength of themagnetic field vertically lies in simulation box, as we see in Figure 4.3. Thus, we
can use the vertical wavelength kz = 2π/H. Since the laminar state is steady, we
can evidently neglect time-evolution terms in basic equations. In the concrete, we put
ω = 0. Since ω and k are given in the dispersion relation ω = ω(k, η), we obtainthe saturated resistivity in the steady laminar state, ηsat . ηsat also gives the saturated
Elsasser number Λ. Moreover, since η is a function of current density J in our model(see Equation (4.10)), we can obtain also the saturated current density, Jsat . Thus,
finally, we will obtain an analytical solution of the Maxwell stress in the saturated
laminar state because the current is closely related to magnetic field.

First, let us derive resistivity in the saturated state, ηsat . We use the dispersion relation
derived from the linearized equations in Sano & Miyama (1999),
4
2
2
vAz0
kz4 + Ω2 ηsatkz2 − 3Ω2 vAz0
kz2 = 0,(4.35)

where we have assumed incompressible fluid, a sheet approximation, and a uniform vertical magnetic field. The uniform magnetic field is justified by the fact that fluctuationof Bz is much smaller than Bz0 . Solving Equation (4.35) with respect to ηsat , we obtainη as function of kz ,√ηsat =

2
3vAz0
v4
− Az0
.
kz
Ω2

(4.36)

Using kz = 2π/H in Equation (4.36), η is analytically given by

η=

2
vAz0

√

Ω

2
vAz0

√

βz0
−1≈
βcrit
Ω

βz0
,
βcrit

(4.37)

where βcrit is defied as βcrit = 8π 2 /3 and we have used βz0 > βcrit ≈ 26.3. Thus,
2
dividing Equation (4.37) by vAz0/Ω, we obtain the Elsasser number in the saturated

state asv2
Λη, sat = Az0 =
ηΩ

√
βcrit
= 0.051
βz0

(

βz0
104

)−1/2
.

(4.38)We show dependence of the time- and volume-averaged Λη, sat on ⟨⟨J⟩⟩ in Figure 4.13
in order to verify this formula. We confirm that at J ≲ Jturb , the formula is consis-

4.3. PREDICTIVE FORMULAE

85

101

0

Λη sat

10

34
35
36
37
38
39
40
41
42

10-1

10-2 -3
10

10-2

10-1
〈〈J〉〉 [Ju]

100

101

Figure 4.13: Dependence of Λη, sat on ⟨⟨J⟩⟩. The dashed line shows Equation (4.38).

tent with the result. When the formula well reproduces the result, the fluctuation issmall because this formula is derived only by assumption of the linearization of basic
equations. Because the nonlinearity of the fluctuations does not appear, we can expect
this formula to strongly predict Λη, sat and other physical quantities. Note that in usual
the dispersion relation with ω = 0 leads to the shortest unstable wavenumber (critical√wavenumber), kcrit = 3vAz0 /η. Thus, the assumption of kcrit = 2π/H means that MRI
growth stops at the laminar state because unstable wavenumbers can not exist in thesimulation box. We can also derive Equation (4.37) from a equation 2πH = kcrit .
Next, we give the Maxwell stress as a function of current density J. We obtain arelationship between Bx and By from Sano & Miyama (1999),
Bx =v2
By = Az0ηsat Ω(

2
vAz0(−2δBy ) = −2Λη,sat By ,ηsat Ω

1 1 Ω2−2
4 4 vAz0
kz2)(
Bx = Λη,sat

1 1 Ω2−2
4 4 vAz0
kz2(4.39))
Bx .

(4.40)

These relations give the ratio between Bx and By as, dividing Equation (4.39) byBy ,Bx= −2Λη,sat = −0.1
By

(

βz0
104

)−1/2
.

(4.41)

86

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

Here, we define the Maxwell stress normalized by P0 in the laminar state as αM, lami .
From Equation (4.39), αM is given by
αM, lami = Λη,sat

1
2B 2 .
4πP0 y(4.42)

Ampére’s equation J = c/(4π) (∇ × B) can current density to magnetic field. Substituting kz for ∇, we obtain the relation between Jx and By ,Jx ≈ −

ckz By .
4π

(4.43)

Because a current sheet is formed by anti-parallel By , Jx dominates J. Thus, we canreasonably assume thatJ ≈ Jx .

(4.44)

Thus, using Equations (4.38), (4.42), (4.43), and (4.44), we obtain αM, lami as a functionof Jsat ,αM, lami

)2
(
4
Jsat
= 2 2 Λη,satk H
Ju()−1/2 (
)2
βz0
Jsat
≈ 0.52,104
Jmax(4.45)

where we have used kz = 2π/H and normalized Jsat by the maximum current density
Jmax = 10Ju which is expected to be the saturated current in fully developed MRI
turbulence. This formula can be shown from the scaling relation of Equation (4.32),
substituting Λ−2η of Equation (4.38) into Equation (4.32).

Then, we present maxwell stress in our model mimicking the real nonlinear Ohm’s law.
The J–E relation of our model at J > JEH is given by the power function with theindex γ (Equation (4.8)), J ∝ E γ . Since we have already obtained the saturated ηsatin Equation (4.37), the specified J–E relation can gives the saturated current density.
Once saturated current density Jsat is obtained, Equation (4.49) gives the Maxwell
stress αM which is the most interesting value. Thus, we here derive a specific formulaof αM when the J–E relation is given by the power function at saturated state. In the4.3. PREDICTIVE FORMULAE

87

power function of the J–E relation, ηsat is given by, (Equation (4.10)),(
ηsat (J) = η0

Jsat
JEH

) γ1 −1.

(4.46)

2
Multiplying this equation by Ω/vAz0, the inverse is given by(Λη,sat = Λ0

Jsat
JEH

)1− γ1(4.47)

Using Equation (4.44), we obtain(Λη, satΛ0

Jx = JEH

γ
) γ−1

.

(4.48)

Substituting this equation for Jsat in Equation (4.49), we obtain(

αM, lami

2γ
) 1−γ
Λ0
≈ 0.52Λη, sat2γ()−1/2 (
)2 (( )()1/2 ) 1−γ
βz0
JEH
Λ0
βz0
≈ 0.52
195
4
10
Jmax
10
104

βz0
104

)−1/2 (

JEH
Jmax

)2 (

(4.49)(4.50)We have used γ = 0.1. Here, using γ = 0.1, the Maxwell stress in our model is writtenas(αM, lami ≈ 1.68

βz0
104

)−1/2 (

JEH
Jmax

)2 (

βz0
104

)−7/18 (

Λz, 0
10

)2/9
.

(4.51)

This predictive formula of the Maxwell stress in the laminar state is available when
J–E relation during electron heating is given by a power function with the index γ =0.1.

In usual, a nonlinear state of MRI can not be expressed by such the analytic expression
because of the nonlinearity of the fluctuations. Thus, when electron heating eﬀectively
completely dies away MRI, the accretion stress will be exactly given by these quantities.

88

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING
100
Vertical Velocity Dispersion [cs]

34
35
36
37
38
39
40
41
42

10-1
10-2
10-3
10-4
10-5
10-6

10
10

0

5

10

15
20
Time [Orbits]

25

30

0

34
35
36
37
38
39
40
41
42

-1

〈〈vz〉〉 [cs]

10-2
10-3
10-4
10-5
10-6 -3
10

10-2

10-1
〈〈J〉〉 [Ju]

100

101

Figure 4.14: Time variation of volume-averaged vertical velocity ⟨vz ⟩normalized bysound speed cs (upper panel) and dependence of time- and volume-averaged resistivity⟨⟨vz ⟩⟩ on time- and volume-averaged ⟨⟨J⟩⟩ (lower panel). The colors are the same as in
Figure 4.5.

4.4

Discussion

4.4.1

Vertical Diﬀusion of Dust

Here we describe the implication for dust dynamics in the disks. Under the classical planetesimal formation theories, the sedimentation of dust forms a dusty layer on
midplane which is gravitationally unstable. Then, the dust layer causes gravitational
instability and forms planetesimals. This model has been focused in terms of avoiding
meter size barrier. However, disk turbulence easily stir up the dust layer and diﬀuse it. Therefore, weak disk turbulence may help the planetesimal formation.

4.4. DISCUSSION89

To investigate the potential of the laminar state to help dust settling, we investigate
vertical velocity at each parameter. The upper panel of Figure 4.14 shows the time
variation of vertical velocity, and the lower panel of Figure 4.14 shows dependence ofthe time-averaged one on ⟨J⟩. We can see that vertical velocity drops at ⟨J⟩ /Ju ∼ 10−1
where the state of gas motion changes from turbulence to laminar flow. To settle down,
dust grains requires
St >

δvz2
,
c2s

(4.52)

where St is the Stokes number, and δvz2 is velocity dispersion. Under the minimum-masssolar nebula in a hydrodynamic equilibrium, St can be expressed as(
St = 1 × 10

−5

a
1 µm

)(

r )3/2
exp
10 AU

(

)
z2
.
2H 2

(4.53)

When the disk flow is laminar or J < Jturb , even 0.1 µm sized dust grains settles on
the midplane dust layer. The laminar state in e-heating zones may help formation of adust layer and planetesimal formation there.

4.4.2

Weak Disk Turbulence Suggested by Recent Observations

Electron heating would suppress the disk turbulence in a large region, but the direct
observability of our theory is still unknown. Therefore the theory might be supportedwith a indirect manner. If electron heating is widely eﬀective in protoplanetary disks,
the disk turbulence is expected to be weak. Unfortunately, an observation of disk
turbulence requires a high spectral resolution to distinguish Keplerian velocity (∼ 100cs )and non-Keplerian velocity (∼ 0.1cs or ≲ 0.01cs , weak turbulence). Thus, the disk
turbulence in protoplanetary disk is actually vague, but the turbulence strength has
been expected to be α ∼ 0.1–0.01 from the accretion stress allowing lifetime of thedisks.
However, recent disk observations might possibly imply the weak disk turbulence. The
disk around HL Tau, which is thought to be the most typical protoplanetary disk
surrounding a T Tauri star (a young solar-type star), is recently observed by ALMA

90

CHAPTER 4. SUPPRESSION OF MRI BY THE ELECTRON HEATING

observatory having extremely high resolution, and then the amazing detailed figure
is unveiled with the high resolution (ALMA Partnership et al., 2015). The disk hasaxisymmetric many rings and gaps approximately within 100 AU from the star. Becausesuch the young system (≤ 1-2 Myr; Briceño et al., 2002) is thought to be dynamic, many
researchers discuss the formation mechanism. Because gas giant planets are well known
to make a gap, the gaps might be formed by gas ginats. However, according to a coreaccretion model of classical planet formation, gas giants take time to get their envelope
from the surrounding disk gas (several Myr). In addition, if each of the gaps is cavedby a gas giant, planets are orderly arranged with the separation of each planets of ∼ 10
AU. Thus, the orbital stability has been discussed (ALMA Partnership et al., 2015).
The other possibilities to explain the rings and gaps is some (gravitational) instabilities(Takahashi & Inutsuka, 2014) and dust growth with/without a planet (Zhang et al.,2015; Okuzumi et al., 2015).
It goes without saying that the gap formation mechanism is debated, but the disk
properties of HL Tau are of great interest too. Pinte et al. (2016) reproduced thesimilar observational image with the radiative transfer simulation and obtained thedust and gas properties. According to the paper, such a clear gap requires for the
disk to be geometrically thin, which means the weak turbulence (α ≲ a few 10−4 ).
That is because if the disk is thick, the back of the disk should be observed. Kataoka
et al. (2015) constrained the maximum dust particle size with a size distribution andobtained relatively small maximum particle size (∼ 150 µm), using the dust continuum
data of ALMA and the polarization data of the other observatories. Because even the
maximum size of 3 mm assumed by Pinte et al. (2016) leads to α ≲ a few 10−4 , themaximum size of Kataoka et al. (2015) should lead to much smaller α value.
Moreover, Flaherty et al. (2015) observed a disk around A-type star, HD163296, andobtained the spectral map which limits on non-thermal gas velocity dispersion, i.e.
turbulence. They constraints the velocity dispersion of turbulence to be less than
∼ 0.03cs which corresponds to α ≲ 10−3 outside ∼ 30 AU, from CO(3-2) transition.
The value is one order of magnitude less than typical α value of fully developed MRI
turbulence α ∼ 10−2 . However, we have to note that the observation with the transitionline observe a sparse region of disks (≲ 5×10−3 g cm−2 or 3–5 H in their model) and the4.4. DISCUSSION

91

precise information of disk turbulence in a dense region has still not been obtained. In
future, a observatory having higher spectral resolution will clearly reveal the turbulence
distribution of the disks.

Chapter 5
Summary and Conclusion

5.1

Summary

In this work, we investigated the quantitative eﬀect of electron heating on the MRI with
MHD simulation including the increasing resistivity by the electron heating. First, we
presented the zones where electron heating can suppress MRI turbulence and showedthat the zones dominate large regions that have been is in a vigorous turbulent state
so far. However, how much electron heating suppresses the MRI was still unknown.
We perform MHD simulation to address the nonlinear evolution of MRI including theresistivity by electron heating. We introduce a simple analytic resistivity model thatenables us to treat the increasing resistivity. We summarize the results as follows:• We find an clear relation between magnetic turbulence strength and its current
density in Figure 4.8 in Section 4.2. The relation means that a lower current leads
to lower turbulence strength.. Thus, we confirmed the electron heating suppresses
MRI turbulence when current density is well suppressed by electron heating.
• We find that when turbulence completely dies away, laminar accretion flow iscaused by ordered magnetic field in Figure 4.3 in Section 4.2. The growth ofmagnetic field stops at the sinusoidal structure because the resistivity increases
at the current sheet,• Based on the simulation results and the scaling relation between the accretion
93

94

CHAPTER 5. SUMMARY AND CONCLUSION
stress and current density, we obtain a formula that successfully predicts theaccretion stress in the presence of electron heating once current density is given. In
a protoplanetary disk, the current density can be estimated by using the currentdependent resistivity and the saturated resistivity. Thus, we can predict theaccretion stress in extensive regions where the electron heating occur. This helps
us to construct realistic planet formation theory.

This paper have revealed the importance of electron heating on the nonlinear evolutionof MRI. Disk turbulence in protoplanetary disks critically impacts on the disk evolution
and planet formation. A role of electron heating should be elucidated in order toconstruct realistic planet formation models.
In this work, we have neglected the stratified density structure. However, our results
might depends on the simulation box size in a vertical direction. We have to notethat these results is based just on the simple analytic J–E relation. In reality, theelectron heating rises the N-shaped nonlinear Ohm’s law. Of course, finally, the realeﬀect of electron heating on MRI is investigated by considering the realistic resistivity
model with a MHD simulation. However, this work first unveils the relation of Maxwell
stress and current density and the state when MRI completely dies away by suppressing
current density. The importance of current density in MRI is stressed. We conclude
that the reliability of the suppression of MRI by electron heating has been shown by
this work.

5.2

Future WorkWe will address the following in future work.

Development of General Nonideal MHD Eﬀects Amplified Electron Heating
In this study, we have neglected the all nonideal MHD eﬀects amplified the electron
heating except for the Ohmic dissipation. We will develop an applicable manner to Halleﬀect and ambipolar diﬀusion. This work is done with Satoshi Okuzumi.

5.2. FUTURE WORK

95

Inspection of MRI Activity in Disk around HL Tau
Recently, Pinte et al. (2016) suggested weak turbulence in the disk around HL Tau. So
far young disk is thought to be dynamic because of high accretion rate onto the starand free fall to the disk from the envelope. The accretion rate of the disk surrounding
HL Tau was estimated 8.7 × 10−7 M⊙ yr−1 as an upper limit (Beck et al., 2010). This
corresponds to α ∼ 0.01 at 100 AU, whereas Pinte et al. (2016) suggested α ≲ 10−4 toreproduce the observation. This seems to be inconsistent, but we have to note that thetwo α parameters shows diﬀerent physics. α in mass accretion rate comes from accretionstress, whereas α in Pinte et al. (2016) comes from vertical velocity dispersion. Thus,
the gap of the two suggest high Schmidt number, which is a ratio of accretion viscosity
to diﬀusion coeﬃcient. Because our work suggested that laminar state leads to highSchmidt number (see Section 4.4.1), the laminar state of our work is good for thispoint.
We will first investigate the possibility that MRI in the disk is stabilized by nonideal
MHD including suppression of MRI by electron hating. First, using some suggested disk
profile (density, temperature and so on), we evaluate dimensionless number indicatingeﬃciency of nonideal MHD eﬀects. In that case, the eﬀects includes Ohmic dissipation,
Hall eﬀect, and Ambipolar diﬀusion, but not include electron heating. Then, we strict
the parameter (magnetic field strength, X-ray luminosity, and so on) to reproduce theweak turbulence in a whole region of HL Tau’s disk. We will also evaluate the αparameter suggested from several nonideal MHD simulations. After that, we introduce
the electron heating eﬀect. Comparing the two restricted parameters with respect toaccretion rate and diﬀusion coeﬃcient, we will investigate the plausibilities.

Formation of Ring Structure in Disk around HL Tau
In addition to the weak turbulence has a potential to explain a formation mechanism of
the ring structure. The observed ring width is approximately 10 AU. Although several
formation mechanisms are proposed, there is no consensus. One possible formation
mechanism is secular gravitational instability (Youdin, 2005a,b, 2011; Takahashi &Inutsuka, 2014), which is gravitational instability caused by gas friction between gas

96

CHAPTER 5. SUMMARY AND CONCLUSION

and dust. Takahashi & Inutsuka (2014) found the most unstable wavelength in theradial direction is typically 10 AU in weak turbulence. If the instability grows globally,
the rings whose width is 10 AU forms.
We will investigate possibility of secular gravitational instability in disk around HLTau. The weak turbulence has already been suggested (Pinte et al., 2016), but the
distribution of turbulence strength is unknown. We calculate turbulence strength simply
assuming the existence of dead zone and e-heating zone. This work will be performed
with Takahashi Sanemichi.

Gravitational Instability Far from the Star
In protoplanetary disks which is one kind of accretion disks, the disk gas accretes onto
the star. If the disk turbulence presents, the accretion is driven by turbulent viscosity closely related to turbulence strength. Thus, in the weak turbulence region, e.g.
dead zone, accreted gas is accumulated and eventually causes gravitational instability(Sano et al., 2000; Martin & Lubow, 2011). The dead zones, where MRI is completely
suppressed, lies within ∼ 20 AU in typical disk at early stage. Since electron heating
enlarges the zone where MRI is completely or moderately suppressed, gravitational instability would occurs in a larger region. This might suggests gas giant far from thestar which is formed by the gravitational instability. The formation mechanism of suchgas giants is still mysterious.
We will investigate how electron heating enlarges the gravitationally unstable region.
The stress–current relation presented in this work gives accretion stress in each position.
Assuming the steady and gravitationally stable accretion disk driven by the Maxwell
stress, we calculate dependence of accretion rate on surface density within the gravitationally unstable surface density, i.e. Toomre’s gravitational unstable condition. Next,
assuming a gas accretion from the outer disk by eﬀective turbulent viscosity, we givethe constant accretion rate. If the stable disk is realizable, there is a surface density
realizing the constant accretion rate. In other words, if there is no surface density realizing the constant accretion rate, the disk goes unstable. We will investigate whether
the disk can be stable or not at each position.

5.2. FUTURE WORK

97

Simulation with the Real Nonlinear Ohm’s Law
We will try to consider the real dependence of resistivity on E. In MHD equation
system, the growth of electric field strength can not be explicitly treated. In addition,
when dσ/dE < 0, electric field strength is unstable. The time scale is remarkably fast,i.e. ∼ 10−4 –10−2 s−1 . Thus, the real nonlinear Ohm’s law is not available without any
idea.
We will develop the physically exact solution by focusing the existence of physics ofmore slower time. Actually, the nonlinear Ohm’s law is obtained by assuming thecharge equilibrium. However, in such the short time interval, the equilibrium is broken.
Thus, we use typical time scale of charge reaction.

Simulation in Stratified Box
In this work, we perform MHD simulation with unstratified disk. Although density
varies in the vertical direction in a hydrostatic equilibrium in reality, the unstratified
disk is thought to be justified by assuming that the perturbation mode is limited in gas
scale height H. Thus, if unstable mode is limited within typically scale hight even in
the stratified simulation box, our results is reliable in the stratified case. However, ourresults would be controlled by the size of simulation box (see Figure 4.11). Thus, theour result will be modified. In order to investigate the stratified case, we will performthe similar MHD simulation again in the stratified simulation box.

98

CHAPTER 5. SUMMARY AND CONCLUSION

Acknowledgement
First of all, I would like to express my gratitude to my advisor, Dr. Satoshi Okuzumi,
for his constant encouragement and patient support since I have been thought fromhim. His guidance over the sea have always properly helped me in the current year.
Dr. Shigeru Ida tolerantly accept me to go away from Titech as a visiting student,
nevertheless he generously gives me financial support. Besides my advisors, I would liketo also thank Dr. Shu-ichiro Inutsuka, Dr. Takeru Suzuki, and Dr. Hiroshi Kobayashi
for their kindly accepting me as a visiting student at Nagoya University. I am gratefulto Dr. Takayuki Muranushi for his kindly giving me the Athena code modified by myrequest and his friendly technical supports for it. I am also grateful to Dr. Neal Turner,
Dr. Mario Flock and Dr. Yasuhiro Hasegawa for the very fruitful discussions andthe kind support during my visiting JPL. I thank also my colleague in our laboratory,
Takahiro Ueda and Yuhito Shibaike, for stimulating discussions and also my roommatesin Ta laboratory of Nagoya University. Computations were performed on “maddux” atTa laboratory of Nagoya University.

Bibliography
Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756
ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
Bai, X.-N. 2014, ApJ, 791, 137
Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737
Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
Beck, T. L., Bary, J. S., & McGregor, P. J. 2010, ApJ, 722, 1360
Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163
Blum, J., & Wurm, G. 2008, ARA&A, 46, 21
Briceño, C., Luhman, K. L., Hartmann, L., Stauﬀer, J. R., & Kirkpatrick, J. D. 2002,
ApJ, 580, 317
Carballido, A., Cuzzi, J. N., & Hogan, R. C. 2010, MNRAS, 405, 2339
Carballido, A., Stone, J. M., & Pringle, J. E. 2005, MNRAS, 358, 1055
Chen, F. F. 1977, INTRODUCTION TO PLASMA PHYSICS 『プラズマ物理入門』,
ed. 内田岱二郎 (丸善株式会社)
Davydov, B. 1935, Sowjet.
Desch, S. J. 2004, ApJ, 608, 509
Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531
99

100

BIBLIOGRAPHY

Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, The AstrophysicalJournal, 765, 114
Ferguson, E. E. 1973, Atomic Data and Nuclear Data Tables, 12, 159
Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99
Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908
Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735,
122
Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343
Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751
Gammie, C. F. 1996, ApJ, 457, 355
Ganguli, B., Biondi, M. A., Johnsen, R., & Dulaney, J. L. 1988, Phys. Rev. A, 37, 2543
Glassgold, A. E., Najita, J., & Igea, J. 1997, ApJ, 480, 344
Golant, V. E., Zhilinsky, A. P., & Sakharov, I. E. 1980, Fundamentals of plasma physics(Wiley New York)
Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051
Gressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291
Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
Hawley, J. F., & Stone, J. M. 1998, ApJ, 501, 758
Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
Hershey, A. V. 1939, Physical Review, 56, 916
Inutsuka, S., & Sano, T. 2005, ApJ, 628, L155
Johansen, A. 2009, Lecture notes on the magnetorotational instability
Kataoka, A., Muto, T., Momose, M., Tsukagoshi, T., & Dullemond, C. P. 2015, ArXiv
e-prints, arXiv:1507.08902
Kataoka, A., Muto, T., Momose, M., et al. 2015, 1504.04812

BIBLIOGRAPHY

101

Kunz, M. W., & Balbus, S. A. 2004, MNRAS, 348, 355
Landau, L. D., & Lifshitz, E. M. 1953, Physical Kinetics (New York: Pergamon)
Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56
Lifshitz, E. M., & Pitaevskii, L. P. 1981, Physical Kinetics (Pergamon Press)
Martin, R. G., & Lubow, S. H. 2011, ApJ, 740, L6
Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
Matthews, L. S., Land, V., & Hyde, T. W. 2012, ApJ, 744, 8
Mayor, M., & Queloz, D. 1995, Nature, 378, 355Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398
Mori, S., & Okuzumi, S. 2016, ApJ, 817, 52
Muranushi, T., Okuzumi, S., & Inutsuka, S. 2012, ApJ, 760, 56
Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199
Nakano, T., & Umebayashi, T. 1986, MNRAS, 218, 663
Natta, A., Testi, L., Neri, R., Shepherd, D. S., & Wilner, D. J. 2004, A&A, 416, 179
Okuzumi, S. 2009, ApJ, 698, 1122
Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65
Okuzumi, S., & Inutsuka, S. 2015, ApJ, 800, 47
Okuzumi, S., Momose, M., Sirono, S., Kobayashi, H., & Tanaka, H. 2015, ArXiv eprints, arXiv:1510.03556
Okuzumi, S., Tanaka, H., & Sakagami, M. 2009, ApJ, 707, 1247
Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M. 2011, ApJ, 731, 96
Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413
Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25
Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15

102

BIBLIOGRAPHY

Rodmann, J., Henning, T., Chandler, C. J., Mundy, L. G., & Wilner, D. J. 2006, A&A,
446, 211
Safronov, V. S. 1969, Evolution of the protoplanetary cloud and formation of the earthand planets.
Salmeron, R., & Wardle, M. 2003, MNRAS, 345, 992
Sano, T., Inutsuka, S., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321
Sano, T., & Miyama, S. M. 1999, ApJ, 515, 776
Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486
Shukla, P. K., & Mamun, A. A. 2002, Introduction to Dusty Plasma Physics (Instituteof Physics Publishing Ltd., Bristol)
Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ,
775, 73
Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ,
764, 66
Simon, J. B., Hawley, J. F., & Beckwith, K. 2009, ApJ, 690, 974
Spitzer, Jr., L. 1941, ApJ, 93, 369
Suzuki, T. K., Muto, T., & Inutsuka, S. 2010, ApJ, 718, 1289
Takahashi, S. Z., & Inutsuka, S. 2014, ApJ, 794, 55
Testi, L., Natta, A., Shepherd, D. S., & Wilner, D. J. 2003, A&A, 403, 323
Turner, N. J., Carballido, A., & Sano, T. 2010, ApJ, 708, 188
Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
Turner, N. J., & Sano, T. 2008, ApJ, 679, L131
Turner, N. J., Sano, T., & Dziourkevitch, N. 2007, ApJ, 659, 729
Umebayashi, T. 1983, Progress of Theoretical Physics, 69, 480
Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617

BIBLIOGRAPHY

103—. 1990, MNRAS, 243, 103—. 2009, ApJ, 690, 69
Velikhov, E. P. 1959, JETP
Wadehra, J. M. 1984, Phys. Rev. A, 29, 106
Wannier, G. H. 1953, Bell System Technical Journal, 32, 170
Wardle, M. 1999, MNRAS, 307, 849
Wardle, M., & Ng, C. 1999, MNRAS, 303, 239
Wardle, M., & Salmeron, R. 2012, MNRAS, 422, 2737
Weidenschilling, S. J. 1977, MNRAS, 180, 57—. 1984, Icarus, 60, 553
Weidenschilling, S. J., & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H.
Levy & J. I. Lunine, 1031–1060
Wolk, S. J., Harnden, Jr., F. R., Flaccomio, E., et al. 2005, ApJS, 160, 423
Yoon, J.-S., Song, M.-Y., Han, J.-M., et al. 2008, Journal of Physical and Chemical
Reference Data, 37, 913
Youdin, A. N. 2005a, ArXiv Astrophysics e-prints, astro-ph/0508659—. 2005b, ArXiv Astrophysics e-prints, astro-ph/0508662—. 2011, ApJ, 731, 99
Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7
Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513,
A57

Appendix A
Appendix

A.1

Time Evolution of Simulation Results

We enumerate results of section 4.2 with respect to time-dependent value.

105

106

APPENDIX A. APPENDIX

Figure A.1: Time evolution of magnetic field strength in the case of JEH /Ju = 0.03.
Each panel shows time, t = 0 (upper left), 10 (upper right), 20 (lower left), and 30(lower right) orbits. The colors shows the magnetic fields strength |B| normalized by
the maximum value.

A.1. TIME EVOLUTION OF SIMULATION RESULTS

Figure A.2: Same as Figure A.1, but in the case of JEH /Ju = 0.3.

107

108

APPENDIX A. APPENDIX

Figure A.3: Same as Figure A.1, but in the case of JEH /Ju = 3.

A.1. TIME EVOLUTION OF SIMULATION RESULTS

Figure A.4: Same as Figure A.1, but in the ideal MHD case.

109

APPENDIX A. APPENDIX
-1

10

-2

10

-3

10

-4

34
35
36
37
38
39
40
41
42

10-5
10-6
10-8
10-9

10-10
0

5

10

15
20
Time [Orbits]

25

10-2
Vertical Kinetic Energy [P0]

10

0

34
35
36
37
38
39
40
41
42

10-1

10-7

10-11

10
Azimuthal Kinetic Energy [P0]

10

34
35
36
37
38
39
40
41
42

-3

10-4
10-5
10-6

10-2

30

10-7

0

5

10

15
20
Time [Orbits]

25

100

30

34
35
36
37
38
39
40
41
42

Kinetic Energy [P0]

Radial Kinetic Energy [P0]

110

10-1

10-8
10-9

10-10
10-11
10-12

0

5

10

15
20
Time [Orbits]

25

10-2

30

0

5

10

15
20
Time [Orbits]

25

30

Figure A.5: Time-evolution of kinetic energy in x (upper left), y (upper right), z (lowerleft) and the total (lower right) direction. The colors show computational ID whichmeans diﬀerence of JEH .
0

10-2
10-4
10

-6

10-8

10-10
10

34
35
36
37
38
39
40
41
42

10-2
10-4
10

-6

10-8

10

0

5

10

15
20
Time [Orbits]

25

10-2

30

Vertical Magnetic Energy [P0]

34
35
36
37
38
39
40
41
42

10-3

10-4

0

10-10

-12

10-14

10
Azimuthal Magnetic Energy [P0]

34
35
36
37
38
39
40
41
42

-12

10-14

0

5

10

15
20
Time [Orbits]

25

100

Magnetic Energy [P0]

Radial Magnetic Energy [P0]

10

10

30

34
35
36
37
38
39
40
41
42

-1

10-2

10-3

0

5

10

15
20
Time [Orbits]

25

30

10-4

0

5

10

15
20
Time [Orbits]

25

30

Figure A.6: Time-evolution of kinetic energy in x (upper left), y (upper right), z (lowerleft) and the total (lower right) direction. The colors show computational ID whichmeans diﬀerence of JEH .

1

10

0

10

34
35
36
37
38
39
40
41
42

-1

10-2
10-3

0

5

10

15
20
Time [Orbits]

25

30

Vertical Current Density [Ju]

34
35
36
37
38
39
40
41
42

0

-1

10-2

34
35
36
37
38
39
40
41
42

10-1
10-2

10-6

0

5

10

15
20
Time [Orbits]

25

102
10

1

10

0

30

34
35
36
37
38
39
40
41
42

10-1

10-3

10

10-4

-2

10-3

10-5
10-6

10

0

10-5

101

10

1

10-4

10-5

10

10

111

10-3

10-4

10-6

Azimuthal Current Density [Ju]10

Current Density [Ju]Radial Current Density [Ju]A.1. TIME EVOLUTION OF SIMULATION RESULTS

0

5

10

15
20
Time [Orbits]

25

30

10-4

0

5

10

15
20
Time [Orbits]

25

30

Figure A.7: Time-evolution of current density in x (upper left), y (upper right), z(lower left) and the total (lower right) direction. The colors show computational ID
which means diﬀerence of JEH .

10-5

34
35
36
37
38
39
40
41
42

Resistivity [ηu]

10-6

10-7

10-8

0

5

10

15
20
Time [Orbits]

25

30

Figure A.8: Time variation of η. The colors shows the same as Figure A.5.

112

APPENDIX A. APPENDIX

10-4

Emax(β=104)
34
35
36
37
38
39
40
41
42

Electric Field Strength [Eu]

10-5
10-6
10-7
10-8
10-9

10-10
10-11
10-12

0

5

10

15
20
Time [Orbits]

25

30

Figure A.9: Time variation of E. The colors shows the same as Figure A.5.

105

Elsasser Number Λ

10

34
35
36
37
38
39
40
41
42

4

103
102
101
100

10-1
10-2

0

5

10

15
20
Time [Orbits]

25

30

Figure A.10: Time variation of Λ. The colors shows the same as Figure A.5.

A.1. TIME EVOLUTION OF SIMULATION RESULTS

10-5

113

34
35
36
37
38
39
40
41
42

Resistivity [ηu]

10-6

10-7

10-8

0

5

10

15
20
Time [Orbits]

25

30

Figure A.11: Time variation of η. The colors shows the same as Figure A.5