Magnetic Field in Star Formation
Kohji Tomisaka
National Astronomical Observatory Japan
August 20, 20221

Copyright ©2022 by Tomisaka, K.
version 1.0: August 20, 2022
This is prepared for a lecture presented in Bootcamp associated with XVIIIth Rencontres du Vietnam conference “Star Formation in Different Environments” held at The International Center ofInterdisciplinary Science Education (ICISE), Quy Nhon, Vietnam, August 21-27, 2022. If you
find any problems in this note, I would appreciate very much if you let me know it by e-mail:kohji_tomisaka@yahoo.co.jp. This was compiled at 00:59 of 18th Aug, 2022.
Download from https://astrophys.sakura.ne.jp/Lecture Notes/Magnetic Field/1.0.pdf

Contents
1 Introduction
1.1 Lorentz Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Equation of Motion and Maxwell Stress Tensor . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Gravity Flux Tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Induction Equation and Flux Freezing . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Magnetic Reynolds Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Flux Freezing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Hydrostatic Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Magnetohydrostatic Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
2
2
3
4
4
4
6

2 Cloud Supported by Magnetic Field
2.1 Virial Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Virial Analysis of Isothermal Cloud (B=0) . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Magnetohydrostatic Cloud (B̸=0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Observation of Magnetic Field Strength . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Zeeman Splitting Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Davis–Chandrasekhar-Fermi Method . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Observation of Mass-to-Flux Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
7
8
9
10
11
12

3 Angular Momentum Transfer by Magnetic Field
13
3.1 Magneto-driven Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4 Magnetic Diﬀusion
15
4.1 Ionization Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.2 Ambipolar Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.3 Dissipation of Magnetic Field by Ambipolar Diffusion, Hall Effect, and Ohmic Dissipation ∗ 19
5 Dynamical Collapse of Spherical Isothermal Clouds
22
5.1 Preprotostellar Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2 Protostellar Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6 Example of Numerical Simulation
24
6.1 Evolution of Rotating Cloud with Magnetic Field . . . . . . . . . . . . . . . . . . . . . 24
6.1.1 Preprotostellar Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.1.2 Protostellar Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.1.3 U-type and I-type Outflows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.2 Angular Momentum Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
7 Aligned vs. Misaligned Rotation

32

A Maxwell Equation

40

B Fluid Dynamics
41
B.1 What is fluid? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
B.2 Equation of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
B.3 Lagrangian and Eulerian Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
B.4 Continuity Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42i

B.5 Expression for Momentum Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.6 Energy Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.6.1 Polytropic Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.6.2 Energy Equation from the First Law of Theromodynamics . . . . . . . . . . . .
C Polytrope

43
44
44
44
46

D Bonnor Ebert Mass
48
D.1 Hydrostatic Equilibrium in Tenuous Gas . . . . . . . . . . . . . . . . . . . . . . . . . . 48
D.2 Bonnor-Ebert Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
E Free-fall Time52

F Alfvén Wave

54

G Magnetic Braking
55
G.1 Aligned Rotator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
G.2 Perpendicular Rotator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
H Self-Similar Solution of Dynamical Contraction
57
H.1 Derivation of Equation for Self-Similar Variables∗ . . . . . . . . . . . . . . . . . . . . . 57
H.2 Self-Similar Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
H.3 Inside-out Collapse Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59ii

1

Introduction

1.1

Lorentz Force

Plasma containing magnetic field B feels Lorentz force as
FB =

j×B
,
c

(1.1)

where j and c represent respectively the electric current density and the speed of light. In this lecture
note, I use cgs Gauss unit to describe electro-magnetic equations (in MKSA unit the Lorentz force isexpressed F = j × B). Ampère’s law with displacement current is expressed as
∇×B =

1 ∂E
4π
j+
,
c
c ∂t

(1.2)

(see Appendix A). Neglecting the displacement current (the second term of the righthand side), theAmpère’s law is written ascj=∇ × B.
(1.3)
4πThis yields the Lorentz force expressed only with the magnetic field as
FB =(∇ × B) × B
.
4π

(1.4)

Equation (1.4) shows an important character of the magnetized material that the Lorentz force does
not work in the direction of B but works in the perpendicular direction both to B and to ∇ × B.
Using a vector analysis formula∇A2 = 2(A · ∇)A + 2A × (∇ × A),(1.5)

the Lorentz force of equation (1.4) is rewritten as
FB = −

1
1∇|B|2 +
(B · ∇)B,
8π
4π

(1.6)

where the first term is called magnetic pressure force and the second is magnetic tension force.
Figure 1 displays a typical configuration in which the magnetic pressure and the magnetic tension
are working. When there exists a gradient in the strength of magnetic field, the magnetic pressure
force appears [Fig.1 (left)]. On the other hand, when the magnetic field changes in the direction of
the magnetic field, the magnetic tension force works [Fig.1 (right)].

Figure 1: (Left): magnetic pressure works toward the direction of −∇|B|2 . (Right): the magnetic
tension force works similar to the force when a magnetic field line is regarded as an elastic
string.

11.2

Equation of Motion and Maxwell Stress Tensor

Considering the Lorentz force the equation of motion for the (fully ionized) plasma is given as


∂v
1
1
dv
≡ρ
+ (v · ∇) v = −∇p + ρg −∇|B|2 +
(B · ∇)B,
(1.7)
ρ
dt
∂t
8π
4π
where g represents the gravity and −∇p represents the pressure force. Using the continuity equation∂ρ+ ∇ · (ρv) = 0,
∂t

(1.8)

equation (1.7) is rewritten∂ρvi∂
1 ∂|B|21+(ρvi vj + pδi,j ) = ρgi −+∂t∂xj
8π ∂xi
4π



∂
Bj
∂xj


Bi ,

(1.9)

where i and j take 1, 2, or 3 corresponding to x, y, and z and δi,j represents the Kronecker’s δ. We
used the Einstein’s summation convention as
X ∂Aj
∂Aj
=
.
(1.10)∂xj∂xjj

The quantity appeared in the second term of lefthand side equation (1.9), ρvi vj + pδi,j is called
momentum flux density tensor. Lorentz force is reduced to
1 ∂|B|2
1
∂
∂
1
1
2
−
+
Bj
Bi =
− |B| δi,j +Bi Bj ,
(1.11)
8π ∂xi
4π
∂xj
∂xj
8π
4π
1
Bi Bj −
where we used ∇ · B = 0 or Bi (∂Bj /∂xj ) = 0. The quantity appeared in this equation 4π
12δ
iscalledMaxwell
stresstensor.
Equation(1.9)
iswritten
as
|B|
i,j
8π


∂ρvi
∂
1
∂
1
2
+
Bi Bj −
|B| δi,j ,(1.12)(ρvi vj + pδi,j ) = ρgi +∂t∂xj∂xj 4π8π

1.2.1

Gravity Flux Tensor1

For the gravity, gravity force of ρg is also expressed using the gravity flux tensor as


1 ∂
1 2
ρgi = −
gi gj − |g| δi,j .
4πG ∂xj
2(1.13)

This kind of expression is used in MHD simulation codes such as Athena++. That is, gravity flux
tensor is expressed as
1
1
gi gj − |g|2 δi,j .
(1.14)−
4πG
2This is understood as follows: writing down the tensor, equation (1.13) is equivalent to

 
 1 2

gx − gy2 − gz2
gx gy
gx gz
gx


2

1
∂
∂
∂
1
2
2
2


gy gz
gx gy
ρ  gy  = −2 −gx + gy − gz
 .
4πG ∂x ∂y ∂z
1
2
2
2
gz
gx gz
gy gz
2 −gx − gy + gz
(1.15)
1

Skip this section unless you are interested in the formulation of numerical schemes.

2Figure 2: The change of flux anchored to a surface dS. One part of the change comes from thetemporal change in magnetic flux density. The other comes from the change of the boundaryof the surface due to gas motion. Time variation of the area is given by vdt × ds.

For example, x-component of this equation is written as





1
∂ 1 2
∂
∂
2
2
ρgx = −
(gx gy ) +(gx gz )g − gy − gz +4πG ∂x 2 x
∂y
∂z




 
∂gy
1
∂gx
∂gx ∂gz=−gx ∇ · g −
−
gy +
−
gz .
4πG
∂x
∂y
∂z
∂x

(1.16)(1.17)

When ∇ × g = 0, since ∇ · g = ∇ · (−∇ψ) = −4πGρ, the righthand side of equation (1.17) is identicalto ρgx .
Thus, if we use both the Maxwell stress tensor and the gravity flux tensor, the momentum density
equation is shown as




∂ρvi
∂
1 ∂
1 2
∂
1
1
2
+
(ρvi vj + pδi,j ) = −gi gj − |g| δi,j +Bi Bj −
|B| δi,j . (1.18)∂t∂xj
4πG ∂xj2∂xj 4π8πThat is, the time derivative of momentum density is determined with divergence of the momentum
flux tensor, the gravity flux tensor, and the Maxwell stress tensor.

1.3

Induction Equation and Flux Freezing

Time evolution of the magnetic field is governed by the Faraday’s induction equation [eq.(A.2)],∇×E =−

1 ∂B
.
c ∂t

Using the electrical conductivity σ, the electric field E is related to the electric current j as


v×B
′
j = σE = σ E +,c(1.19)(1.20)

where E ′ = E + (v × B)/c represents the electric field observed with the plasma moving with avelocity v (Lorentz transformation). This is the simplest form of generalized Ohm’s law for ionizedmaterial. Using the Ohm’s law [eq. (1.20)], the induction equation (1.19) becomes∂B
c2 2= ∇ × (v × B) +
∇ B,
∂t
4πσ

(1.21)

where we used a vector analysis formula of ∇ × (∇ × B) = ∇(∇ · B) − ∇2 B = −∇2 B and assumed
σ =const.
3

1.3.1

Magnetic Reynolds Number

The ratio of the first term to the second one of the rhs of equation (1.21) gives the ratio of induction
to diffusion, which is called the magnetic Reynolds number,Rm =|∇ × (v × B)|
c224πσ ∇ B=

U B/L
c2 B/4πσL2=UL
c2 /4πσ≡UL,η(1.22)

where η represents the magnetic diffusivity and U and L are typical size and speed of a system. In
the fully ionized plasma, since the electrical conductivity σ is high (resistivity is low) the system has
high magnetic Reynolds number, Rm ≫ 1. In this case, the magnetic diffusion is negligible comparedwith the induction, and equation (1.21) is simply reduced to∂B= ∇ × (v × B) .
∂t

(1.23)In this case, global structure of the plasma is governed by equations (1.8), (1.7), and (1.23). Thissystem is called ideal magnetohydrodynamics (ideal MHD). When the magnetic diffusion term is notignorable, the system is called as the non-ideal or resistive MHD.
1.3.2

Flux Freezing

An important nature of the ideal MHD is derived from equation (1.23), that is, the matter and themagnetic field is coupled with each other (Alfvén’s theorem [1942]). Consider Ra gas parcel threadedwith the magnetic field B. The magnetic flux anchored to the gas parcel ΦB ≡ B · dS changes withaccordance to
Z
Z
Z
d
∂B
dΦB
=
B · dS =· dS + B · v × ds,(1.24)dt
dt∂t
where s = ∂S. The first term of the rhs of equation comes from the time derivative of the magnetic
flux density, while the last term represents the change of the
R integral regionR due to fluid motion (seeFig. 2). ∵ dS
=
(vdt)
×
ds.
Using
the
Stokes’theorem∇ × A · dS = A · ds, the last term isRR
reduced to (B × v) · ds = ∇ × (B × v) · dS. Finally, we obtain

Z 
∂B
dΦB
=
− ∇ × (v × B) · dS = 0,(1.25)
dt∂t
where we used equation (1.23). This means that the magnetic flux anchored to a gas parcel does notchange in the ideal MHD regime.

1.4

Hydrostatic Balance

Consider a plane-parallel disk extending infinitely in the x-y plane (Fig.3). Using isothermal soundspeed cis = (kT /µmp )1/2 , static balance of a gas obeying the isothermal equation of state, p(z) =(kT /µmp )ρ(z) = c2is ρ(z) is given as
dp(z)
dρ(z)
dψ(z)
= c2is
= −ρ
dz
dz
dz

(1.26)

where ψ(z) represents the gravitational potential. When the gas is self-gravitating, Gauss’s theorem
for the gravitational potential is written as−dψ(z)= gz = −2πGσdz
4(for z > 0),(1.27)

Figure 3: Schematic view of plane-parallel sheet in which pressure gradient and the self-gravity isbalanced each other. The gravity at z is proportional to the column density σ contained
between −z and z.

where gz and σ represent respectively the gravitational acceleration and the column density within±|z| as
Z
Z
z

z

σ=

ρdz = 2ρdz.

(1.28)dρ= −2πGρσ,dzd 1 dρ
4πG∴= − 2 ρ.
dz ρ dzcis(1.29)−z0Thus, equation (1.26) is reduced toc2isThis equation (1.29) gives the solution
z,
cis /(2πGρc )1/2

 ρ 1/2zc
σ(z) = 2,
cis tanh
2πG
cis /(2πGρc )1/2


ρ(z) = ρc sech2(1.30)

where ρc represents the density at the midplane of the disk, ρ(z = 0). Since only when |z| → ∞ gas
density goes to zero, the gas disk extends to infinity, although its total column density is finite.
When the disk is confined by the external pressure pext , the disk extending to a point where thepressure p(z) is balanced with pext . The density at the surface is given as
ρs =

pext
.
c2is

Integrating equation (1.29) from the center z = 0 to the surface,
Z σs
πG 2
2
cis (ρs − ρc ) = −πGσdσ = −σ .
2 s
0(1.31)(1.32)

When density contrast between the center and the surface is large ρc ≫ ρs ,
ρc =

πG 2
σ .
2c2is s

The central density ρc is proportional to the total column density σs squared.
5(1.33)

1.5

Magnetohydrostatic Balance

Next, we consider a B-field parallel to the isothermal self-gravitating disk (see Fig. 1) asBx (z)
B =  0 .
0(1.34)

Balance is now expressed asdp(z)1 dBx2 (z)
dψ(z)
+
= −ρ
,
(1.35)dz
8π dz
dz
which contains the magnetic pressure term of Lorentz force equation (1.4). Similar to equation (1.29),
magnetostatic balance is obtained byc2isdρ
1 dBx2 (z)+= −2πGρσ.
dz 8π dz

(1.36)

Integrating this,1
πG 2(Bs2 − Bc2 ) = −
σ .
8π
2 s
Neglecting the surface density and its magnetic flux density (ρc ≫ ρs and Bc2 ≫ Bs2 ),c2is (ρs − ρc ) +c2is ρc +1 2 πG 2
B =
σ .
8π c
2 s

(1.37)(1.38)

This shows that thermal pressure c2is ρc and the magnetic pressure Bc2 /8π are balanced against thegravitational pressure (πG/2)σs2 .
Comparing equations (1.33) and (1.38), to support the same column density σs , the magnetic
model requires smaller central density compared with the hydrostatic model 2 . Magnetic pressure
gradient −(1/8π)(dBx2 /dz) has an effect to support the gas against the self-gravity.

2
Reducing the magnetic ﬁeld strength Bc ↘ does not lead magnetohydrostatic balance impossible. Reducing themagnetic ﬁeld Bc ↘ results in increasing the central density ρc ↗. However, be careful that this is a story of theplane-parallel isothermal disk.

6

2

Cloud Supported by Magnetic Field

2.1

Virial Analysis

Hydrodynamic equation of motion using the Lagrangean time derivative [eq. (B.3)] is dvρ= −∇p − ρ∇Φ.
dt

(2.1)

Multiplying the position vector r and integrating over a volume of a cloud, we obtain the virial
relation as
1 d2 I= 2(T − T0 ) + W,
(2.2)
2 dt2
where
Z
I = ρr2 dV,
(2.3)

1 2
3
3
T =
pth + ρv dV = p̄Vcl ,22
2
is a term corresponding to the thermal pressure plus turbulent pressure,Z
3
T0 =pth r · ndS = p0 Vcl
2
S
is the inertia of the cloud,Z 

comes from a surface pressure, p0 , andW =−Z
ρr · ∇ΦdV = −3 GM 2
5 R

(2.4)

(2.5)

(2.6)

is the gravitational energy. To derive the last expression in each equation, we have assumed the cloudis spherical and uniform3 . Here we use a standard notation as the radius R, the volume Vcl = 4πR3 /3,
the average pressure p̄, and the mass M .

2.2

Virial Analysis of Isothermal Cloud (B=0)

Figure 4: (Left): external pressure p0 is plotted as a function of cloud radius R [eq. (2.9)]. (Right):
the same as the left panel but for the magnetic cloud [eq. (2.13)].

To obtain a condition of mechanical equilibrium, we assume d2 I/dt2 = 0. Equation (2.2) becomes
4π p̄R3 − 4πp0 R3 −3

3 GM 2
= 0.
5 R

(2.7)

The numerical factor depends on the density distribution ρ(r). The numerical factor of equation (2.6), 3/5, is forthe uniform density model.

7

Assuming gas is isothermal p = c2is ρ, the average pressure is written asp̄ = c2is ρ̄ = c2is

3M
.
4πR3

(2.8)

Using equation (2.8) to eliminate p̄ from equation (2.7), the external pressure is related to the massand the radius as
3c2 M
3GM 2
p0 = is 3 −
,
(2.9)
4πR
20πR4
which is displayed in Figure 4 (left). Keeping M constant and increasing R from zero, p0 increases
first, but it takes a maximum, p0,max = 3.15c8is /(G3 M 2 ), and finally declines. This indicates that thesurface pressure must be equal to or smaller than p0,max for a cloud to be in the equilibrium. In otherwords, keeping p0 and changing R, it is shown that M has a maximum value to have a solution. The
maximum mass is equal to
c4is
Mmax = 1.77.
(2.10)
1/2
G3/2 p0
Clouds massive than Mmax cannot be supported against the self-gravity. This corresponds to the
Bonnor-Ebert mass, which is the maximum mass of isothermal spherical gas cloud confined by theexternal pressure [eq. (D.9)]4 , although the numerical factors are slightly different.

2.3

Magnetohydrostatic Cloud (B̸=0)

Consider here the effect of the magnetic field. In the magnetized medium, the Lorentz force
F =

1
1
1
(∇ × B) × B = − ∇(B · B) +(B · ∇)B
4π
8π
4π(2.11)

works in the ionized medium, where we used vector formula of ∇(A · B) = (B · ∇)A + (A · ∇)B + A ×(∇ × B) + B × (∇ × A). The first term of equation (2.11), which is called the magnetic pressure,has an effect to support the cloud against the self-gravity.
The virial analysis is also applicable to the magnetohydrostatic clouds. The terms related to themagnetic fields are
Z
Z
Z
B2
B2
M =
dV + (r · B)B · ndS −r · ndS
8π
S 8π
S


Z
Φ2B
Φ2
B 2 − B02
1≃dV ≃ 2− B ,(2.12)
8π
6π
RR0
where ΦB represents
a magnetic flux and it is assumed to be conserved if we change the radius, R,
R
that is ΦB ≡ B · ndS = πB0 R02 = πBR2 . Equation (2.7) becomes
4π p̄R3 − 4πp0 R3 −3 GM 21 Φ2+ 2 B = 0,
5 R
6π R

Φ2

(2.13)

where we ignored the term 6π1 2 RB0 . The last two terms have the same dependence on R and thus they
are rewritten as3G
M 2 − MΦ2 ,(2.14)
5R
where MΦ is defined as 3GMΦ2 /5 = Φ2B /6π 2 .
This shows the effects of the magnetic fields:
4

In the context of stellar structure, the same mass is called Schönberg-Chandrasekhar limit.

8

1. B-fields effectively reduce the gravitational mass as M 2 −MΦ2 = M 2 −5Φ2B /(18π 2 G). This plays
a part to support a cloud.
2. However, even when a cloud contracts (decreasing its radius from R0 to R), the ratio of thegravitational to the magnetic terms keeps constant since these two terms are proportional to∝ R−1 . Thus, if the magnetic term does not work initially, the gravitational term continues topredominate over the magnetic term.
If M < MΦ , a sum of last two terms in equation (2.13) is positive. Since the second term of rhs ofequation (2.9) is positive, there is one R which satisfies equation (2.13) irrespective of the externalpressure p0 . While, if M > MΦ , there is a maximum allowable external pressure p0 . Therefore,
M = MΦ gives a criterion whether the magnetic fields work to support the cloud or not. Morerealistic calculations [Mouschovias (1976a, 1976b), Tomisaka et al. (1988)] give us a criterion
G1/2dm
G1/2 σ1== 0.17 ≃
,
dΦB
B
2π

(2.15)

where σ and B means the column density and the magnetic flux density. A cloud with a massΦB
2πG1/2
is sometimes called magnetically supercritical, while that with
M>

M<

ΦB
2πG1/2(2.16)(2.17)

is magnetically subcritical.
More precisely speaking, the criterion showed in equations (2.16) and (2.17) should be applied fora cloud which has a much larger mass than the Bonnor-Ebert mass. That is, even without magnetic
fields, the cloud less-massive than the Bonnor-Ebert mass has a hydrostatic configuration shownin Figure 28 (left). The cloud with central density of ρc = 10 has a stable density distribution.
Magnetohydrostatic clouds with different magnetic fluxes are calculated by Mouschovias & Spitzer(1976), Mouschovias (1976a, 1976b) and Tomisaka et al. (1988)]. Mass of the cloud is obtained
against the central density, which is an increasing function of magnetic flux (Fig. 5). The maximum
allowable mass (critical mass) supported by some magnetic flux increases with the magnetic flux. To
fit the numerical results, Tomisaka et al. (1988) obtained an expression for the critical mass when
the cloud has a mass-to-flux ratio dm/dΦB , the isothermal sound speed cis , and the external pressure
p0 as
(

2 )−3/2
c4is
1/2π
.
(2.18)
Mcr = 1.3 1 −1/2
G1/2 dm/dΦB |r=0
p0 G3/2
This shows that the critical mass is a decreasing function of the mass-to-flux ratio or increasing
function of the magnetic flux. And the critical mass becomes much larger than the Bonnor-Ebert
1/2
mass ≃ c4is /(p0 G3/2 ) only when the mass-to-flux ratio at the center of the cloud is reaching 1/2π
at which the term in the curry bracket goes to zero. Hereafter, we call here the cloud/cloud corewith mass larger than the critical mass Mcr a supercritical cloud/cloud core. The cloud/cloud coreless-massive than the critical mass is subcritical.

2.4

Observation of Magnetic Field Strength

To determine whether a cloud is magnetically subcritical or supercritical, we have to measure thestrength of the magnetic field.
9Figure 5: Magnetohydrostatic cloud mass is plotted against the central density. Respective curves
correspond to clouds with different magnetic flux. Plasma beta of the external medium(β0 ) specifies the magnetic flux. Model with β0 = ∞ corresponds to that shown in Fig. 28(right).

2.4.1

Zeeman Splitting Method

If the atoms and molecules are permeated by magnetic field B, in the normal Zeeman splitting, thetransition splits into three lines
ν1 = ν0 − β|B|,(2.19)

ν2 = ν0 ,(2.20)

ν3 = ν0 + β|B|,(2.21)

where β = eℏ/2me = 1.3996 Hz µG−1 is the Bohr magneton constant. This comes from the fact thatthe upper and lower levels are splitted into Um = ml βB using the magnetic quantum number ml . Hi
λ = 21cm, OH (maser emissions and thermally excited emissions), H2 O maser emissions are used tomeasure the magnetic field strength. Figure 6 (right) shows the B-field strength of the line-of-sightcomponent is plotted against the hydrogen atomic column density (NH ≡ NHi + 2NH2 ). This shows
the correlation between the magnetic field strength and gas column density.

10Figure 6: (Left): mass-to-flux ratio (λc ≡ µc ) is plotted against the column density (NH2 ). The y-axis
represents the mass-to-flux ratio normalized with the critical one. Dots are obtained fromthe Zeeman splitting measurements (Blos ) and stars are based on the Davis–ChandrasekharFermi method (Bpos ). Taken from Crutcher (2004). (Right): magnetic strength plottedagainst the hydrogen column density. The B-field strength is obtained by Zeeman observations of Hi atom, OH, and CN molecules. Taken from Crutcher (2012).

2.4.2

Davis–Chandrasekhar-Fermi Method

Davis (1951) and Chandrasekhar & Fermi (1953) proposed an idea to obtain the strength of interstellar magnetic field using the character of the Alfvén wave. (The Alfvén wave is a transverse wave
propagating in plasma described by MHD, in which magnetic tension force works as the restoring
force. See Appendix F.) Applying the same method to the molecular cloud, Crutcher et al. (2004)
obtained the magnetic field strength in the plane of the sky in three prestellar clouds as ≃ 80 µG(L183), ≃ 140 µG (L1544), and ≃ 160 µG (L43). The method is as follows: Using the phase velocityof transverse Alfvén wave (Alfvén speed), VA = B/(4πρ)1/2 , the wave is described as
y = A cos [k(x − VA t)] ,(2.22)

where y is a displacement in the plane of the sky and perpendicular to the mean magnetic field. Thisgives
 2  2
∂y
2 ∂y
VA
=
,
(2.23)∂x∂t
where X represents the time average. The rhs of the equation, which is the one-dimensional velocity
dispersion in y direction, seems to be equal to the velocity dispersion in the line-of-sight ∆V as


∂y
∂t

2
= ∆V 2 .

(2.24)The lhs of the equation is estimated not by time average but by spatial average. Dispersion ofobserved polarization direction ∆ϕ (rad) gives the dispersion in ∂y/∂x as


∂y
∂x

2
= ∆ϕ2 ,

11

(2.25)

in which the lhs is measured by the wave-pattern of the Alfvén wave. Therefore, the magnetic fieldstrength on the plane of the sky is given as
Bpos = (4πρ)1/2

2.5

∆V
.
∆ϕ

(2.26)

Observation of Mass-to-Flux Ratio

Figure 6 shows the mass-to-flux ratio normalized with the critical value (Crutcher 2004; 2012)
µc ≡
where



M
Φ

(M/Φ)obs,(M/Φ)crit


=
crit

1
2πG1/2(2.27)(2.28)

is the critical mass-to-flux ratio which is obtained numerically from equation (2.15). Condition ofµc > 1 represents a magnetically supercritical cloud while µc < 1 represents a magnetically subcritical
cloud.
Figure 6 indicates cloud cores are more or less found near the critical mass-to-flux ratio. They
are not distributed either in the regions µc ≫ 1 (very supercritical) or µc ≪ 1 (very subcritical).

123Angular Momentum Transfer by Magnetic Field

Angular momentum changes by a torque exerted on the object, N = F × R. Considering theaxisymmetric system, the torque is given as Nz = Fϕ r, where we use the cylindrical coordinate(r, z, ϕ) in this section. We divide the magnetic field and the electric current into poloidal ones
Bp = Br e r + Bz e z , j p = j r e r + j z e z .

(3.1)

Bt = Bϕ e ϕ , j t = j ϕ e ϕ ,(3.2)

and toroidal ones
where er , ez , and eϕ represent the unit vectors of the cylindrical coordinate.
The ϕ-component of the Lorentz force comes from the poloidal electric current and the poloidal
magnetic field as
1(3.3)
Fϕ e ϕ = j p × B p .
cUsing the Ampère’s law [eq.(1.3)], the poloidal electric current is given as∂Bϕc 1 ∂rBϕ
jp =
er −
ez .
(3.4)
4π r ∂r∂z
Thus, for the torque to work in the axisymmetric system, both poloidal and toroidal components of
B-field are necessary.

3.1

Magneto-driven Model

Figure 7: Magneto-centrifugal wind model. (Left): ordinary effective potential for a gas parcel withconserving the specific angular momentum, j = r2 ω = const. The gas parcel is assumed tolocated originally at (r0 , 0). The x- and y-axes mean normalized distance from the central
star as r/r0 and z/r0 . The darker color represents the deeper potential. (Right): theeffective potential for a gas parcel with conserving the angular rotation speed ω = ω0 = constof the disk [eq. (3.6)].

To accelerate gas till supersonic speed, the magnetic force is considered to be an important player.
Blandford & Payne (1982) have proposed the “magneto-centrifugal wind model” for the accelerationin two-dimensional axisymmetric configuration.
13

Figure 7 (right) explains the magneto-centrifugal wind mechanism, which corresponds to Fig. 1 ofBlandford & Payne (1982). Consider a disk under the gravity of a central star with mass M∗ . Assume
that gas in the disk at the radius r0 rotates with a Kepler speed5 ofvϕ (r0 ) =

GM∗
r0

1/2
.

(3.5)

Consider a gas parcel whose original position was (r = r0 , z = 0). The effective potential for theparticle rotating with a constant angular rotation speed of ω0 = vϕ (r0 )/r0 is written as
"
  #
GM∗
1 2 2
1 r 2
GM∗
r0
ϕ(r, z) = −− r ω0 = −
+
,
(3.6)
r0
(r2 + z 2 )1/2 2(r2 + z 2 )1/2 2 r0
3/2

where ω0 = (GM∗ )1/2 /r0 . Be careful that the ordinary effective potential is calculated for constantangular momentum. The second term is ordinarily written as +(1/2)(j 2 /r2 ) using the specific angularmomentum j ≡ rv = r2 ω. However, equation (3.6) is for constant angular rotation speed. The gas
will co-rotate with the same angular rotation speed as the disk, if the magnetic field is sufficiently
strong. Equation (3.6) represents the effective potential for the gas co-rotating with the disk at r0 .
Figure 7 (right) plots the isopotential contour lines for the effective potential [eq. (3.6)]6 . The contour
line passing the point (r0 , 0) represents ϕ = −3GM∗ /2r0 . The angle between the isopotential line of
ϕ = −3GM∗ /2r0 and the r-axis is equal to ±60◦ (see the problem below). Further, it can be obtainedthat the effective potential increases upwardly while it decreases with reaching the central star orradially outwardly.
This leads to an important conclusion as follows:
1. If the angle between the magnetic field line and the r-axis, θmag , is larger than 60◦ (θmag > 60◦ ),gas parcel moving along the field line must climb the effective potential. It needs extra energy
to depart from the disk.
2. If θmag < 60◦ , gas parcel moving along the field line slides down the effective potential. It can
depart from the disk without extra energy.
Be careful that this is valid if the magnetic field is sufficiently strong. If the Kepler disk is threadedwith the magnetic field lines with a shallow angle, and if the angular momentum is sufficiently
transferred along the magnetic field line, in this case, the gas is accelerated by the extra centrifugal
force and finally escapes from the potential well of the central star. This mechanism is called the“magneto-centrifugal wind” mechanism.
Problem 1
Show that the angle between the isopotential line which passes the point (r0 , 0) and the r-axis
is equal to ±60◦ . Calculate the Taylor expansion of the effective potential near the point (r0 , 0)and show that to satisfy ϕ(r0 + ∆r, ∆z) = ϕ(r0 , 0), (∆z/∆r)2 must be equal to 3.


We have shown that the magnetic field has an effect to extract angular momentum of rotating gas.

5
6

Rotation speed centrifugally balanced with the gravity of the point mass source is called Kepler speed.
In Fig.7 (left), eﬀective potential for constant angular momentum is plotted, for comparison.

14

4

Magnetic Diﬀusion

In §2, we have seen that the mass-to-flux ratio controls whether the self-gravity is balanced withmagnetic force. In this section, we see the magnetic diffusion processes which change the mass-to-flux
ratio of the gas.

4.1

Ionization Rate

Figure 8: Ionization fractions X/nH are plotted against the number density nH . The ionization rateof an H2 molecule by cosmic rays outside the cloud is taken ζ0 = 1 × 10−17 s−1 . 20% of C
and O and 2% of metallic elements are assumed to remain in the gas phase and the rest in
grains. M+ represents the metallic ions such as Mg + , Si+ , and Fe+ , while m+ represents+
molecular ions other than H+
3 , e.g., HCO . Taken from Nakano, Nishi, & Umebayashi(2002).

In the dense clouds, the ionization fraction is low. Since the uv/optical radiations from stars can
not reach the cloud center, potential ionization comes from the cosmic ray particles. In this case therate of ionization per volume is given as ζnn , where ζ represents the ionization rate due to the cosmic
ray and typically we take ζ ∼ 10−17 s−1 . In Figure 8, the ionization fraction for various density isshown (Nakano, Nishi, & Umebayashi 2002). This clearly shows that the fraction of ions decreases−1/2
approximately in proportion to nH
for nH ≲ 108 cm−3 . This is understood as follows: Equilibrium+balance between one kind of ion m and its neutral m0 ,m0 ↔ m+ + e− ,(4.1)

is considered. The recombination (reaction from the right to left) rate per unit volume is expressed
αne nm+ , while the ionization rate (left to right) per unit volume is ζnm0 , where α means the recombination rate coefficient. Nakano (1984) obtained in the range of 102 cm−3 ≲ n ≲ 108 cm−3
ρi = Cρ1/2
n ,
15

(4.2)

where the numerical factor C = 4.46 × 10−16 g1/2 cm−3/2 .

4.2

Ambipolar Diﬀusion

The ionization fraction of the gas in the cloud is low. Neutral gas and the ions are coupled via mutual
collisions and ions are coupled with the magnetic field. Thus, the neutral molecules, major componentof the gas in mass, are coupled with the magnetic field indirectly via ionized ions. Assuming ρn ≫ ρiand ρ = ρn , equation of motion for the molecular gas isρdvn= −c2is ∇ρ − ρ∇ψ + ρrω 2 er + F ,
dt

(4.3)

where the forces appeared in the rhs represent, respectively, the pressure force, the self-gravity, thecentrifugal force, and the force exerted on the neutral component per unit volume through two-bodycollision with ions. The friction force has a form
F = ρi ρn γ(vi − vn ).

(4.4)

We have assumed ρ ≃ ρn since the mass density of the charged component is much smaller than thatof the neutral one. On the other hand, equation of motion for ions is
ρi

dvi
1
=
(∇ × B) × B − F ,
dt
4π

(4.5)

where we ignored the pressure and self-gravity forces compared with the Lorentz force. If the momentof inertia of the ions is ignored [lhs of eq. (4.5)=0], equations (4.3) and (4.5) give the equation of motion
similar to that of the one-fluid as follows:ρdv1= −c2is ∇ρ − ρ∇ψ + ρrω 2 er +(∇ × B) × B,
dt
4π

(4.6)

where we used the fluid velocity v ≡ (ρn vn + ρi vi )/(ρn + ρi ) ≃ vn .
Ignoring the moment of inertia of ions, equation (4.5) indicates that the friction force should bebalanced with the Lorentz force asF =1(∇ × B) × B.
4π

(4.7)

Using equations (4.2) and (4.4), this gives the drift velocity of ions against the neutrals as
vD ≡ vi − vn ,F=,
ρi ρn γ1=(∇ × B) × B,
4πρi ρn γ
1=(∇ × B) × B.
3/2
4πCγρn(4.8)

The term (∇ × B) × B/4π represents the Lorentz force. Equation (4.8) indicates the drift velocity
vD = vi − vn is toward the direction of the Lorentz force. In magnetized cloud, the Lorentz force
works outwardly, the drift velocity vD also points outwardly. Viewing from the ions and magnetic
field lines, vn − vi = −vD points inwardly. Thus, the mass in a flux tube varies with time; that of the16

Figure 9: Relation central density ρc against mass-to-flux ratio M/Φ. Plane-parallel and axisymmetricmodels are illustrated in the left and right panels, respectively. Crosses represent the relation
obtained from numerical calculation by Mouschovias, Paleologou, & Fiedler (1985). Taken
from Tomisaka, Ikeuchi, & Nakamura (1990).

cloud center increases while that of the cloud surface decreases. The ambipolar diffusion is a process
to change the mass distribution against the magnetic flux tube.
It is useful to derive the timescale of ambipolar diffusion. The characteristic timescale of magnetic
flux loss is given as τA = L/|vD |, where L represents the size of the cloud. This gives
4πCγρ3/2 L|(∇ × B) × B|
3/2 
−2 
2
n
B
L
6
≃ 3 × 10 yr.
104 H2 cm−330 µG
0.1 pcτA ≃(4.9)(4.10)

This ambipolar diffusion timescale is several times longer than the dynamical contraction timescaleor the free-fall timescale as
3π 1/2tﬀ =(4.11)
32Gρ̄
−1/2

n
6
.
≃ 0.5 × 10 yr(4.12)
104 H2 cm−3
Thus, the magnetic flux of the cloud is not be lost during the dynamical contraction (≃ free-fall).
It is shown that the ambipolar diffusion process changes the mass distribution against magnetic
flux. How is the structure affected by the increase of the mass-to-flux ratio at the cloud center?
Here, we will derive the relation between the mass-to-flux ratio and central density and answer thequestion. Numerical calculation is necessary for the axisymmetric (realistic) cloud. However, this canbe obtained for the case of a plane-parallel disk in which the magnetic field is running parallel to the
disk. For the plane-parallel disk the hydrostatic balance is achieved as


B2
∂ψ
∂
p+
= −ρ
= −2πGρσ,(4.13)
∂z
8π
∂z
17

Rz
where σ = −z ρdz. Equation (4.13) is integrated to obtain the relation between the variables at thecenter z = 0 as


Bs2
Bc2
πG 2
− ps +
σ ,(4.14)
pc +
=
8π
8π
2 s
where we used the relation dσ/dz = 2ρ [this is already mentioned at eq. (1.37)]. Quantities with
subscript c and s represent those at the center and at the surface, respectively. Dividing both side ofequation by ps = c2is ρs , this is rewritten as
 
 
 
1 Bc 2 ρs ρc 2
1 Bs 2 ρs
ρc
πGσs2
+
−
=
+ 1.
(4.15)
ρs 8π ρcc2s ρs8π ρsc2s2ρs c2s
1/2

In the second term of lhs, yc ≡ (ρc /Bc )(cs /ρs ) is a nondimensional mass flux ratio at the center andthis is related to the center-to-surface density ratio ρ′c ≡ ρc /ρs . Using these nondimensional variables,
equation (4.15) is expressed as(ρ′ )21
ρ′c + c 2 −= ρ̄′c ,(4.16)
8πyc
8πys2
where ρ̄′c ≡ πGσs2 /(2ρs c2s ) + 1 represents the central density necessary for the disk to be supported
without magnetic field (y → ∞). Figure 9 (left) illustrates this relation. This shows that if the central
mass-to-flux ratio increases (moving upward), the central density increases. In low density, ρ′c ≪ ρ̄′c ,equation (4.16) indicates that ρ′c ∝ yc ∝ ρc /Bc . This means that on the low density regime, thedensity increases but the strength of magnetic field does not increase. Equation (4.16) also indicatesthat yc increase much as yc → ρ̄′c ys when ρ′c → ρ̄′c . This relation between mass-to-flux ratio andthe central density is well fitted to the numerical calculation of plane-parallel disk cloud driven by
the ambipolar diffusion by Mouschovias, Paleologou, & Fiedler (1985). Does this indicate that the
mass-to-flux ratio increases much when the central density increases, for example, from 104 cm−3 to106 cm−3 ?
From this point, quasistatic evolution of the magnetized cloud was first studied by Nakano (1979)
adopting a method seeking for magnetohydrostatic configuration (Mouschovias 1976a, b). Nakano(1979) obtained that the mass to magnetic flux ratio (µ) near the center of the cloud does notincrease much before the central density reaches a critical density above which magnetohydrostatic
configurations no more exist. Paleologou & Mouschovias (1983), Shu (1983), and Mouschovias,
Paleologou, & Fiedler (1985) gave a completely different result. That is, the mass to magnetic flux
ratio near the center increases much when the density becomes larger than n ≳ 106 cm−3 . Critical
discussion has been done between them7 . However, we have now realized that the difference comesfrom the geometry assumed.
Relation of the mass-to-flux ratio to the central density for the axisymmetric realistic (not planeparallel disk) cloud is plotted in Figure 9 (right) (Tomisaka, Ikeuchi, & Nakamura 1990). The central
mass-to-flux ratio increases by the ambipolar diffusion. This increases the central density. Lines
with arrows indicate the quasistatic evolutionary paths driven by the ambipolar diffusion. It is tobe noticed that the relation is completely different from that of the plane-parallel disk. Increase of
the mass-flux ratio is small although the central density increases much. This indicates that in thequasistatic evolution driven by the ambipolar diffusion the mass-to-flux ratio does not increase muchin contrast to the plane-parallel case.
7

One of the reasons why large increase in the mass/ﬂux ratio was favored is understood as follows: There has
been a long-standing “magnetic ﬂux problem of stars” in which the magnetic ﬂux of, say, 1 M⊙ main-sequencestar Φ ∼ 5 × 1026 G cm2 (R∗ /2R⊙ )2 (B0 /3 × 104 G) is much smaller than that of the parent cloud Φ ∼ 3 ×1031 G cm2 (Rcl /0.1 pc)2 (B0 /100 µG). The magnetic ﬂux must be reduced in the star formation process. If the mass-ﬂux
ratio increases much in these density ranges, this might resolve the magnetic ﬂux problem of stars.

18

Figure 10: (Left:) Hall parameters for ions and electrons. Two boundaries come from the conditionof ions’ Hall parameter βi = 1 and that of electrons’ βe = 1, both of which are given by
B/ρ =const. Generally, magnetic field strength B increases in proportion to ρ1/2 in thecontraction. In consequence, contraction proceeds from ambipolar diffusion region to theohmic dissipation region. Taken from Wardle (2007). (Right:) Direction of drift velocity
for various Hall parameter. Taken from Wardle (2004).

In fully ionized plasma, the magnetic field is coupled to the matter. However, in the low ionizedgas the magnetic field is coupled only to the ions. Thus the induction equation of magnetic fieldshould be written∂B+ ∇ × (vi × B) = 0.
(4.17)∂tThis leads to the final expression as()
1∂B(4.18)+ ∇ × (vn × B) + ∇ ×[(∇ × B) × B] × B = 0,3/2∂t
4πCγρnwhere we used the expression of the drift velocity [vD = vi − vn ; eq. (4.8)]. In the isothermal cloud,there is a maximum mass supported by the magnetic field. Since in the subcritical cloud M < Mcr
the magnetic flux escapes from the center, the critical mass decreases with time. After the criticalmass becomes smaller than the actual mass (Mcl > Mcr ), there is no hydrostatic configuration. The
cloud evolves into the supercritical cloud region and experiences dynamical collapse.

4.3

Dissipation of Magnetic Field by Ambipolar Diﬀusion, Hall Eﬀect, and Ohmic
Dissipation∗

Beside the ambipolar diffusion process, there are several processes which induce dissipation of magnetic field in the low ionization plasma. These processes break the ideal MHD condition and flux
freezing character (see §1.3.2). Magnetic dissipation is by the resistivity of the plasma. Measuring
the resistivity, we begin with the drift of charged particles against the neutrals. Drift is controlled
by the balance between the electromagnetic force and momentum exchange with neutral molecules
via mutual collisions (Wardle 2007). Considering ionized species named j (ions and electrons), thebalance is expressed asvj
Zj eE ′ + Zj e × B − mj γj ρvj = 0,(4.19)cwhere E ′ and vj represent the electric field in the frame of neutral molecules, velocity of j-th charged
species (mass mj and charge Zj e) relative to the neutral molecules (drift velocity), respectively.
Momentum exchange rate between j-th charged species and neutral molecules are given using
γj = ⟨σv⟩j /(mj + m).
19(4.20)

′′
Figure 11: (Left:) Definition of E//
and E⊥. (Right:) Electric current J is divided into three parts:
current parallel to B is equal to σ0 E// , that flowing in the direction of electric field perpendicular to B as σ1 E⊥ , and that flowing perpendicular to both E⊥ and B as σ2 [(B/|B|) ×
E⊥ ]. Taken from http://wdc.kugi.kyoto-u.ac.jp/ionocond/exp/icexp.html.

In deriving equation (4.19), we assumed that the balance is achieved in much shorter timescale thanthe macroscopic evolution timescale. Drift direction is controlled by the ratio of the second to thethird terms, which is called Hall parameter asβj =|Zj |eB 1.
mj c γj ρ(4.21)It should noticed that the Hall parameter is dependent on magnetic field strength and the density,
B/ρ. When βj ≫ 1, for j-th ionized species, magnetic force is dominated over the collision withneutrals. These charged particles are coupled with the magnetic field. On the other hand, even forcharged particles, j-th species are decoupled from the magnetic field if βj ≪ 1. In Figure 10 (left),the Hall parameters for ions and electrons are plotted in the B − ρ plane (Wardle 2007).
Ambipolar diﬀusion region: When both βi and βe are much larger than unity (1 ≪ βi ≪ βe ),ions and electrons are coupled with magnetic field.
Hall region: When βi ≪ 1 ≪ βe , electrons are coupled with magnetic field, but ions are not.
Ohmic resistive region: When βi ≪ βe ≪ 1, not only ions but also electrons are decoupled fromthe magnetic field.
The boundary between ambipolar diffusion and Hall regions is given by the condition of ion’s Hallparameter equal to 1. The other boundary between the Hall region and the Ohmic resistive region isfrom βe = 1. From equation (4.21), both boundaries are expressed by B ∝ ρ. Since in the contraction
by the self-gravity magnetic field strength B is amplified approximately in proportion to ρ1/2 , themajor process of magnetic dissipation changes from the ambipolar diffusion region to the Ohmic
resistive region, during the contraction.
P
P
From equation (4.19), the current made by the drift velocity, J =j Jj =
j Zj enj vj , isobtained depending on E ′ as′+ σHJ = σO E//

B
′
′
× E⊥
+ σ P E⊥|B|

(4.22)′ and E ′ represent components of E ′ parallel and perpendicular(generalized Ohm’s law), where E//⊥
to the magnetic field [see Fig. 11 (left)]. The first term represents the current flowing parallel to B,
in which the magnetic field does not play a role. Thus, this conductivity is equal to that measured

20

for no magnetic field (the Ohmic conductivity, σO ). The second term represents the current flowing′ and B, which comes from the charged particles drifting in the direction of
perpendicular both to E⊥′′ is called the Hall conductivity (σ ). The±(E⊥ × B) (the Hall effect). Thus, the coefficient to the E⊥
H′ and its coefficient is called the Pedersen
third represents the current flowing in the direction of E⊥
conductivity (σP ). Although the induction equation is written as ∂B/∂t = ∇ × (v × B) in the ideal
MHD condition, in this case an additional term should be added as∂B= ∇ × (v × B) − c∇ × E ′ .
∂t(4.23)c
Assuming only the Ohmic dissipation, E ′ = σ1O J = 4πσ∇ × B,
O∂B= ∇ × (v × B) − ∇ × (ηO ∇ × B) ,∂t

(4.24)

where we used Ampere’s law J = (c/4π)∇ × B and the magnetic diffusibility is given by the conductivity as ηO = c2 /4πσO . The last term for the Hall resistivity and the ambipolar diffusion are written
ashi−∇ × ηH (∇ × B) × B̂ ,(4.25)n hio∇ × ηA (∇ × B) × B̂ × B̂
(4.26)(Wurster 2016, Tsukamoto et al. 2017), where B̂ ≡ B/|B|. As mentioned in section 4.2, in a
dynamically contracting cloud, major magnetic flux loss expects in relatively high density where Hall
resistivity or Ohmic resistivity mainly operates.

21

5

Dynamical Collapse of Spherical Isothermal Clouds

Figure 12: Evolution of isothermal clouds more massive than the Bonnor-Ebert mass. Density (left)
and velocity (right) distributions are illustrated. Solid lines show the cores of the preprotostellar phase (prestellar core) and dashed lines show those of protostellar phase (protostellarcore). The evolution of the protostellar phase is studied by the sink-cell method, where weassume the gas that entered in the sink-cells is removed from finite-difference grids and add
to the point mass sitting at the center of the sink-cells which corresponds to a protostar.
Taken from Ogino, Tomisaka, & Nakamura (1999).

5.1

Preprotostellar Phase

In 1969, Larson (1969) and Penston (1969) found a self-similar solution which is suited for the dynamical contraction. Figure 12 is a radial density distribution for a spherical collapse of an isothermal
cloud, where the cloud has a four-times larger mass than that of the Bonnor-Ebert mass. We assume the initial density distribution is proportional to that of the hydrostatic equilibrium solution(ρ(r) ∝ ρeq (r)) and the proportional coefficient is taken as 4. Although the figure is taken from a
recent numerical study by Ogino et al. (1999), a similar solution was obtained in Larson (1969). We
can see that the solution has several characteristic points as follows:
1. The cloud evolves in a self-similar way. That is, the spatial distribution of the density [Fig. 12(left)] at t = 0.326 is well fitted by that at t = 0.293 after shifting in the −x and the +ydirections. As for the infall velocity spatial distribution (right panel), only a shift in the −x
direction is needed.
2. The density distribution in the envelope, which is fitted by ∝ r−2 , is almost unchanged. Only
the central part of the cloud (high-density region) contracts.
3. The time before the core formation epoch (the core formation time t0 is defined as the time atwhich the central density increases greatly) is a good indicator to know how high the central
density is. That is, reading from the figure, at t0 − t = 0.003 (t = 0.326) the central density
reaches ρc ≃ 104 and at t0 − t = 0.036 (t = 0.293) the density is equal to ρc ≃ 102 . This shows
22the maximum (central) density is approximately proportional to (t0 − t)−2 , which is reasonablefrom the description of the free-fall time ∝ ρ−1/2 (see p.52 of Appendix E).
This kind of solutions is often called runway collapse solution. In the preprotostellar phase, isothermal
gas exhibits the runway collapse.

5.2

Protostellar Phase

After the central density exceeds a limiting density 106 ρeq (r = 0) = ρlim , we set sink-cells around thecenter and excess mass larger than the limiting density in the sink-cells is removed from the sink-cells
and the removed mass is added to the sink (proto-star’s) mass M∗ asX∆M∗ =(ρi − ρlim ) dVi .
(5.1)
i:sink−cells

Gravity due to the removed mass from the sink cells ϕ = −GM∗ /r is added to the self-gravitational
potential.
Dashed line in Figure 12 represent the density and velocity distributions after t > t0 . Solution
after the protostar formation M∗ > 0 has following several characters:
1. After the protostar is made, inflow velocity becomes large near the center due to the point-massgravity.
2. Inflow velocity is proportional to |vr (r)| ∝ r−1/2 and increases with time from t0 .
3. Density distribution is proportional to ρ(r) ∝ r−3/2 and decreases with time from t0 .
4. Density and velocity evolve in a self-similar fashion, also in this protostellar phase.
These characters are similar to inside-out self-similar solution found by Shu (1978). Although the
Shu’s inside-out solution starts from static singular isothermal solution at t = t0 , this inside-out
solution starts from the previous dynamic runaway contraction. Hydrodynamical contraction is well
expressed by two self-similar solutions which are connected at t = t0 (Whitworth & Summers 1985).

23

6

Example of Numerical Simulation

Figure 13: Explanation of nested grid method.

6.1

Evolution of Rotating Cloud with Magnetic Field

We have described the cloud runaway collapse (preprotostellar phase) and succeeding accretion process (protostellar phase). In the latter phase, there is a possibility that the outflow is driven by theeffect of magnetic Lorentz force. A restricted number of numerical simulation could simulate such
evolution throughout before 2000.
The size of the molecular cores ∼ 0.1 pc is more than 3 × 104 times as large as the typical size ofthe first core ∼ 1 AU. Therefore numerical scheme to calculate such process must have a large spatial
dynamic-range. In the finite-difference scheme, the spatial dynamic range is restricted by the numberof cells. Although at least 3 × 105 × 3 × 105 grid points (in 2D) seem necessary to resolve a factor of
3 × 104 , to solve the MHD equation with the Poisson equation of the self-gravity has not been done.
This is done by the Eulerean nested-grid simulation. This method uses a finite number of grid
systems with different spatial resolutions simultaneously. Coarse grid covers whole the cloud and is
to see a global structure. Fine grid covers only the central part of the cloud and is to see the fine-scale
structure appeared near the center. Figure 13 explains grid cells of the nested grid method. The sizeof the cells of the n-th level grid (Ln) is taken equal to a half of those of the Ln − 1 (Tomisaka 1998,
2000, 2002).
A cylindrical cloud with coaxial magnetic field in hydrostatic balance is studied. Rotation vectorΩ and the magnetic field B are parallel to the symmetric axis of the cylinder (z-axis). For theenergy equation, a double-polytropic equation of state is adopted, since the radiation hydrodynamical
calculation (for example Masunaga & Inutsuka 2000) predicts ρ − T relations which consists of anumber of power-laws for some density ranges.
 2(for ρ < ρA ), cis ρ  
7/5
p=
(6.1)
ρ
 c2is ρA(for ρ ≥ ρA ),
ρA
24

(a)

(b)

(c)

Figure 14: Ωc = 10−13 rad s−1 , Bc = 10 µG, ρc = 104 H2 cm−3 , and cis = 190 (T /10K)1/2 m s−1 .
Snapshots at the time of 0.6066τﬀ = 1.061 Myr represented in different levels are shown:
L1 (upper-left), L5 (upper-right), and L10 (bottom). The actual size of the frames of L5(upper-right) and L10 (bottom) are, respectively, 1/16 and 1/512 smaller than that of
L1 (upper-left). Magnetic field lines (red solid lines) and isodensity contours (black solidlines) are presented as well as the velocity vectors by arrows. Taken from Fig. 1 of Tomisaka(2002).

25

where we take ρA = 1010 H2 cm−3 . In the density range of ρ < ρA , gas behaves like an isothermalgas, while after ρ > ρA , gas behaves like an adiabatic gas. Figure 14 shows the model of Ωc =10−13 rad s−1 , Bc = 10 µG, ρc = 104 H2 cm−3 , cis = 190 (T /10K)1/2 m s−1 , and scale-length H =0.341(cis /190 m s−1 ) (ρc /104 H2 cm−3 )−1/2 pc. This is a snapshot at t = 0.6066 τﬀ = 1.061 Myr. The
cylindrical cloud breaks into pieces with a length in z-direction equal to the most unstable wavelengthof the gravitational instability (Matsumoto, Nakamura, & Hanawa 1994) as
λmax ≃

2π(1 + α/2 + β)1/2

cis
,
i1/2
(4πGρc )1/2
1/3
0.72 (1 + α/2 + β) − 0.6
h

(6.2)

where α ≡ Bz2 (r)/4πc2is ρ(r) and β ≡ −[ρ(r)vϕ (r)2 /r]/[c2is dρ(r)/dr] represent the magnetic-to-thermalpressure ratio and the ratio of the centrifugal force to the thermal pressure force, respectively.
6.1.1

Preprotostellar Phase

The magnetic fields prevent gas motions perpendicular to the field lines B. Further, the centrifugal
force works in the perpendicular direction to Ω. A disk forms which running in the perpendicular
direction to B and Ω and thus to the z-axis. Figure 14 illustrates a snapshot at the time of 0.6066 τﬀ =1.061 Myr. Level 1 grid (upper-left) captures the global structure of the contracting disk perpendicular
to the z-axis. Close up view is shown in L5, which has 16 times finer spatial resolution (upper-right).
Since speed, vz , of infalling gas is larger than the sound speed, a number of shock fronts are formednear z ≃ 0.02H and z ≲ 0.01H. Further inner structure is represented in L10 grid (bottom) which has
32 times finer resolution than L5. Another discontinuity in density is now forming near z ∼ 0.0002H.
Since gas inside of this discontinuity has the density larger than the critical density ρ > ρA , in thisphase an adiabatic core (corresponding to the first core) begins to be formed. In the isothermal
collapse phase, a characteristic power-law for the Larson-Penston self-similar solution is seen in theradial distributions of density as ρ(r, z = 0) ∝ r−2 , and Bz (r, z = 0) ∝ r−1 . A disk formed in
the preprotostellar phase (gas is still isothermal) looks disklike but continues to contract even withrotation, which is called a pseudo-disk.
6.1.2

Protostellar Phase

After t > 0.6067τﬀ the first core8 is formed (protostellar phase). Radial infall motion is accelerated
toward the surface of the first core and the rotational motion vϕ also increases. This amplification
promotes the toroidal magnetic field, Bϕ . Since the torque is proportional to (∇ × Bϕ eϕ ) × Bp , theamplitude of the torque increases after the first core formation.
Structure just after ∼500 yr has passed from Figure 14 is shown in Figure 15. In τ ∼ 500 yr
from the epoch of Figure 14, flow is completely changed and outflow is launched. Panel (b), whichcorresponds to a four-times close-up of panel (a), shows clearly that outflow is ejected from a region
near r ≃ 1.5 × 10−4 H and z ≃ 0.5 × 10−4 H. And the outflow vectors and the disk has an angle ofapproximately ≃ 45◦ . The magneto-centrifugal wind model favors a small angle between the magnetic
field and the disk, θmag < 60◦ . This model predicts that the angle between the flow and the disk is also
8

From radiation hydrodynamical simulation, Larson (1969) explored formation process of protostars; First, a hydrostatic core made with hydrogen molecules H2 is formed, which is called ﬁrst core. The ﬁrst core grows in mass due tothe accretion from the envelope. Along with the mass accretion, central density and thus central temperature increases.
Finally, hydrogen molecules begin to dissociate to hydrogen atoms, H. Since the dissociation process absorbs energy,
second contraction begins from the center. After the dissociation is completed, another hydrostatic core is made, whichis called second core. The second core grows and ﬁnally becomes a protostar.

26

(a)

(b)

(c)
P3

P2

P1

Figure 15: Continuation of the evolution shown in Fig. 14. Snapshots at the time of 0.6069 τﬀ =1.066 Myr just ∼ 500 yr after from Fig. 14. Panel (a) plots the structure of L10 which
is similar to Fig. 14 (c). Spherical region inside ∼ 1.2 × 10−3 H has been swept by theoutflow. Panel (b) is a snapshot of L12 grid which has 4-times finer spatial resolution than
L10 (a). In panel (c), we plotted magnetic force dominated region with ‘+’, centrifugal
force dominated with ‘*’, and thermal pressure force dominated with blank ‘ ’. Taken from
Fig. 3 of Tomisaka (2002).

27

smaller than ≲ 60◦ . Numerical results is not inconsistent with the prediction of magneto-centrifugal
wind model.
To see which force is driving such an outflow, amplitudes of three forces are compared: the Lorentz
force Fm = (∇ × B) × B/4π, the centrifugal force Fc = (ρvϕ2 /r)er , and the thermal pressure gradientFp = −∇p. The parallel components to the poloidal magnetic field of the three forces are as
Fm ·

Bp
|Bp |

(∇ × B) × B Bp
·
,
4π
|Bp |
1 Bp
= −
· ∇(rBϕ )2
8πr2 |Bp |
=

(6.3)

(Ustyugova et al. 1999)
Fc ·

Bp
Br
= ρvϕ2 r ·
,
|Bp |
|Bp |

(6.4)

Fp ·

Bp
Bp
= −∇p ·
.
|Bp |
|Bp |

(6.5)

and

Figure 15 illustrates the largest force at each grid point:
Region C: Centrifugal force dominated region (region C), which is shown by asterisk (*), overlaps
the outflow region around point P1 (z = 10−4 H, r = 2.5 × 10−4 H).
Region M: Just radially exterior to this region C, there is another region near point P2 (z =2 × 10−4 H, r = 3.5 × 10−4 H) filled with plus signs where the magnetic force is dominated(region M).
Comparing with the flow vectors in the panel (b), the strong outflow coincides with this region M andthe above region C. Therefore, the outflow is driven by the centrifugal force and the Lorentz force(pressure gradient by toroidal magnetic fields).
6.1.3

U-type and I-type Outﬂows

Difference between models with different magnetic field strength but the same rotation speed isillustrated in Figure 16. In this comparison, we assume virtually a hard equation of state such asp ∝ ρ2 for ρ ≥ ρA gas in equation (6.1). Magnetic to thermal pressure ratio α is taken 1 [panels(a), and (d )], 0.1 [panels (b), and (e)], and 0.01 [panels (c), and (f )]. Panels (a), (b), (c) in theleft column are snapshots when the central density reaches ρA while panels (d ), (e), (f ) in the rightcolumn show the structure after τ = 4.5 × 10−3 τﬀ = 8000 yr passed. In the model with α = 1, theoutflow gas flows through a region whose shape resembles a capital letter U. This is similar to thatobserved in the model with α = 1 and Ω = 5 shown in Figure 15. Decreasing the magnetic fieldstrength from α = 1 to 0.01, the shape of disk changes from a flat disk to a round ellipsoid. Thisis reasonable, since the model with neither magnetic field nor rotation shows a spherical contractingcore. Outflow is also affected by decrease in magnetic field strength. Globally folded magnetic fieldlines are seen in the outflow region in panel (e). Smaller-scale folding structure of the magnetic fieldlines develops in the outflow region in panel (f ). The structure seen in the poloidal magnetic field
is well correlated to the distribution of toroidal magnetic field strength. That is, a packet of strong
Bϕ pinch the gas and poloidal magnetic field lines radially inwardly. This means that the toroidal
magnetic field is amplified by the rotation motion and becomes predominant over the poloidal one
in the case with weak initial magnetic field. Thus, the whoop stress, which is expressed by theterm of −(Bϕ /4πr)[∂(rBϕ )/∂r] and comes from the tension of toroidal magnetic field, pinches the
28

(a)

(d )

(b)

(e)

(c)

(f )

Figure 16: Difference in models with the same rotation speed Ω = 2 × 10−14 rad s−1 but different
magnetic field strength ((a) and (d ): Bc = 10 µG, (b) and (e): 3 µG, and (c) and (f ):1 µG). Snapshots in the left column corresponds to the epochs when ρmax = ρA , whilepanels in the right column show snapshots after 8000 yr from the left. Taken from Fig. 7of Tomisaka (2002).

29

plasma contained inside of the loop of Bϕ . For this case, the magnetic force dominant region is widely
spread compared with the centrifugal force dominant region, which is far from the model with strong
magnetic field (α = 1). In this case, the magnetic pressure gradient plays an important role. Gas ofoutflow is ejected perpendicularly from the disk and forms flow similar to the capital letter I. In
conclusion, U-type and I-type outflows are generated for strong magnetic model and weak magnetic
model, respectively.

6.2

Angular Momentum Distribution

Figure 17: Evolution of specific angular momentum contained in accumulated mass from the center.
Open circles, filled squares, and triangles denote, respectively, the beginning and the endof the runaway collapse phase and the accretion phase. Taken from Tomisaka (2000).

During this evolution, magnetic torque is working to the outflowing gas and let the gas rotate. As
the reaction, the magnetic torque is working the disk gas to reduce the angular momentum. To see
the time evolution of angular momentum distribution, we plotted angular momentum spectrum 9 inFigure 17 (Tomisaka 2000). Vertical axis is the specific angular momentum and the horizontal axis is
the accumulated mass from the center. In a spherical symmetric calculation, the accumulated massis well defined as
Z r
ρ4πr2 dr,
(6.6)
Mr =
0

and is regarded as an independent variable. However, in a 2D axisymmetric calculation this is notapparently defined. Since the gas density increases monotonically reaching the center, the mass shellwith higher density is considered to be located nearer to the center. Total gas mass whose density is
Angular momentum spectrum means literally that mass ∆M with speciﬁc angular momentum of j − j + ∆j (y-axis)
is plotted against j (x-axis). However, to consider the anglar momentum transfer, i.e., extraction and deposition, we
regard the accumulated mass from the center as the independent variable and its speciﬁc angular momentum as thedependent variable.
9

30

Z

higher than ρ1 is calculated
M (> ρ1 ) =

ρdV,

(6.7)

which contains the total angular momentum (z-component) ofZ
J(> ρ1 ) =ρrvϕ dV.

(6.8)

ρ>ρ1

ρ>ρ1

Specific angular momentum contained in M (> ρ1 ) is calculated asj(> ρ1 ) =J(> ρ1 ).
M (> ρ1 )(6.9)

In Figure 17, we plotted j(> ρ1 ) as a function of M (> ρ1 ).
Before the first hydrostatic core is made, the angular momentum distribution evolves from the linewith open circles to that with filled squares. The central ≃ 0.5 (cis /190 m s−1 )3 (nc /104 H2 cm−3 )−1/2 M⊙
reduces its angular momentum to 1/3 before the core is formed.
Since we found no outflows (the outflow extracts angular momentum from the central part) in thisphase, this reduction seems to be explained by the magnetic braking (Appendix G). This is explained
also by considering that gas with smaller angular momentum is preferentially accumulated during
the contraction and constructs the high density region (Tomisaka 2000).
Angular momentum distribution changes greatly after the protostar formation. In ≃ 7 × 103 yr(nc /104 H2 cm−3 )−1/2 , the distribution evolves from the line with filled squares to that of filled triangles in Figure 17. The central ≃ 0.17 (cis /190 m s−1 )3 (nc /104 H2 cm−3 )−1/2 M⊙ has been reduced toa factor ≲ 10−4 of the initial value. The angular momentum transfer at this stage is mainly occurredby the wind. That is, outward flowing gas brings a finite angular momentum and as a result angularmomentum is transferred outwardly.
In conclusion, the outflow process in molecular cores both with magnetic field and rotation plays
a role in segregating angular momentum and mass. That is, about 10 % of the accreting mass isejected by the magnetic wind mechanism. At the same time, the outflow brings almost all theangular momentum of the originally accreting gas. From the rest 90 % of the gas mass, the angularmomentum is extracted and the gas continues to accrete to the central protostar. The specific angularmomentum of the protostar is reduced to a factor of 10−4 of the initial value. Today, this reduction
process of angular momentum is essential to form a star, avoiding the centrifugal barrier.

31

7

Aligned vs. Misaligned Rotation

Figure 18: Time evolution of the directions of (1) the central magnetic field Bc (thin solid line),(2) the central angular momentum Jc (dotted line), (3) normal vector of the central disk
structure nd (dashed line), along with (4) the magnetic field in the large scale (thick solidline). Diamonds denote the epoch when the adiabatic core formation ρc = ρcr . This shows
three vectors (1), (2), and (3) converge each other. Taken from Fig. 5 of Matsumoto &Tomisaka (2004).

In preceding sections, we studied the evolution of molecular cloud whose angular momentum J andmagnetic field B are initially parallel each other and their direction coincides with the axisymmetricaxis. As shown in Appendix G, in the parallel configuration, the torsional Alfvén wave brings theangular momentum away from the cloud and the wave propagates before the cloud rotation is reduced
to 1/2 of the initial one asρclzswept ≃ h
(7.1)
ρext
where h, ρcl , and ρext represent the height of the cloud, the cloud density, and that of the external
medium, respectively [eq. (G.7)]. In contrast, when J ⊥ Ω, the wave front travels in the same foldingtimeρcl 1/4
rswept ≃ R
,
(7.2)
ρext
where R represents the cloud radius [eq. (G.6)]. Considering high contrast between ρcl and ρext ,we expect rswept ≪ zswept . Thus, we expect the perpendicular rotator J ⊥ Ω looses its angularmomentum more rapidly than the aligned rotator J //Ω (see Mouschovias & Paleologou 1979, 1980).
Using MHD simulations, Matsumoto & Tomisaka (2004) studied the evolution of misaligned
rotating cloud (J //\ Ω). The initial density distribution is similar to the Bonnor-Ebert sphere (Ebert
1955, Bonnor 1956) but its density is enhanced by a factor f . The cloud has the mass of Mcl =32

Figure 19: Final stage of the simulation. (Left:) cross-section along the plane which contains Bc andz-axis (Plane A of the right panel). (Middle:) cross-section along the plane which contains
Bc and x-axis (Plane B of the right panel). Grey scale represents the density. Velocity
component lying in the plane is shown by the arrow and that out of the plane is shownwith red and blue contour lines. (Right:) explanation of Planes A and B. Taken from
Fig. 1 (d ) - (f ) of Matsumoto & Tomisaka (2004).

6.130 M⊙ , radius Rcl = 0.178 pc, the central density ρc = 1 × 10−19 g cm−3 , and the rotation rate
Ωcl = 7.11 × 10−7 yr−1 . In commonly used non-dimensional parameters, the thermal to gravitationalenergy ratio is α0 ≡ Eth /|Egrav | = 0.5 and the rotational to gravitational energy ratio is β0 ≡Emag /|Egrav | = 0.02. Magnetic field and rotation are initially uniform as




sin θ
0
(7.3)
B = B0  0  , Ω =  0  ,
cos θ
Ωcl
where θ represents the angle between B and the z-axis (and thus Ω). Respective models with θ = 0◦and θ = 90◦ correspond to the parallel and perpendicular rotators. Similar to equation (6.1), thedouble-polytropic equation of state is assumed for the gas pressure
 2(ρ < ρcr ), cis ρ  
5/3
(7.4)
p=
ρ
 c2is ρcr(ρ ≥ ρcr ),ρcrwhere ρcr = 2 × 10−11 g cm−3 and cis = 0.19 km s−1 .
Magnetic field and angular momentum at the center of the cloud are defined as
R
BdV
ρ>ρ
Bc = R 0.1max,
ρ>ρ0.1max dVZρr × vdV,
Jc =

(7.5)

ρ>ρ0.1max

where we average the quantities for the central gas with higher density than 1/10 of the maximum
density (ρ > 0.1ρmax ). From the moment of inertia Iij averaged for the same area, we obtain the
shortest principal axis, which gives the direction of the disk-like structure (pseudo-disk), nd .
We begin with the case B0 = 18.6 µG, which corresponds to 1/4 of critical magnetic flux density(mass to magnetic flux ratio µ = 4) and magnetic to gravitational energy ratio is equal to γ0 ≡33

Figure 20: Jc /Mc2 represents the rotation speed normalized by the central free-fall timescale. Solidline represents Jc,// /Mc2 , in which Jc,// is the component of angular momentum parallel tothe local magnetic field, Bc . Jc,⊥ /Mc2 (dotted line) is calculated from the perpendicularcomponent of Jc,⊥ . MF00, MF45, MF90 are the models with θ = 0◦ , 45◦ , and 90◦ . Other
parameters are the same. NM represents the result of nonmagnetic model B0 = 0. Taken
from Fig. 4 of Matsumoto & Tomisaka (2004).

Emag /Egrav = 0.721. Figure 18√shows the√ evolution of Bc , Jc and nd . Initially Bc is parallel◦ to◦ thet
global magnetic field B = (1/ 2, 0, 1/ 2)B0 of equation (7.3), which points to (x, y) = (45 , 0 ) in
the figure. The large-scale B (thick solid line) rotates counter-clockwise as the cloud rotates as awhole. Although Bc (thin solid line) rotates also counter-clockwise, its rotation speed is much fasterthan the large-scale B. Just after the simulation begins, a pseudo-disk is formed in the directionnd ≃ (30◦ , 0◦ ) (dashed line). This vector rotates also counter-clockwise. Although the angularmomentum Jc (dotted line) is set in the z-direction initially, the direction of Jc begins to approach
those of Bc and nd . After the core formation epoch (open diamond), the outflow begins to flow. It
is shown that directions of the three vectors converge to the point ≃ (16◦ , 33◦ ). It should be noticed
that this direction agrees neither with the initial direction of magnetic field nor that of the angularmomentum. In conclusion, the magnetic field, rotation, and the disk normal direction converge eachother finally.
Figure 19 shows the gas flow around the protostar at the final stage of this simulation or at theendpoint of Figure 18. Left and middle panels show 2D cross-cuts, in which density and velocity field
are displayed. Panel (d ; left) corresponds to the plane determined by Bc and the z-axis (plane A of
panel (f )) and panel (e; middle) corresponds to that by Bc and the x-axis (plane B of panel (f )).
Both panels capture the so-called U-type outflow (see p.30). Thus, we can conclude that even fromthe misaligned rotator θ ≲ 45◦ , outflow is ejected, in which flow is almost axisymmetric.
R From the angular momentum of2 the central part Jc [eq. (7.5)] and the corresponding mass Mc =
ρ>0.1ρmax ρdV , the quantity Jc /Mc indicates the normalized central rotation rate, which gives theangular rotation speed ωc normalized with the central free-fall timescale (4πGρc )−1/210 We divide Jc
into the components, of which J// and J⊥ are the angular momentum parallel and perpendicular to10

cis J 2 /(GM ) = cis ωc r2 /(Gσπr 2 ) = ωc /(2πGρc )1/2 , where we used σ = (2c2is ρc /πG)1/2 [eq. (1.33)].

34

the local magnetic field as

Jc,// = |Jc · Bc |/|Bc |,Jc,⊥ = |Jc × Bc |/|Bc |.

(7.6)

In Figure 20, we plotted Jc,// /Mc2 (solid line) and Jc,⊥ /Mc2 (dotted line) as a function of the maximum
density ρmax , which is substituted for the elapse time t, since ρmax is a monotonically increasing
function of t. Figure 20 shows following points:
1. Comparing the solid and dotted lines of MF45 model (θ = 45◦ ), Jc,⊥ decreases faster than Jc,// .
2. Comparing MF00 model (θ = 0◦ ) and MF90 model (θ = 90◦ ), the perpendicular rotator(θ = 90◦ ) loses its angular momentum more rapidly than the aligned rotator (θ = 0◦ ).
3. In NM model (dashed line), Jc /Mc2 increases with ρmax before the adiabatic core formation(open diamond). This is explained by the fact that gas which was originally located far fromthe center and had larger angular momentum falls to the center later. Comparing models NM
and MF00, magnetic braking works in MF00 and reduce the normalized rotation rate to ≃ 1/5.
Matsumoto & Tomisaka (2004) concluded that the perpendicular rotator (θ = 90◦ ) looses its angularmomentum faster than the aligned rotator (θ = 0◦ ).
Joos, Hennebelle, & Ciardi (2012) revisited the same problem. Their setup is a bit different from
Matsumoto & Tomisaka (2004): the density distribution is taken as ρ = ρ0 /[1 + (r/r0 )2 ], the thermal
to gravitational energy ratio α0 = 0.25, the rotational to gravitational energy ratio β0 = 0.03, mass
to flux ratio ratio are in the range from µ = 2 to µ = 17. They obtained the central specific angularmomentum by the formulaR
ρ>ρ1 ρr × vdV
R
j(ρ > ρ1 ) =(7.7)
ρ>ρ1 ρdV
where ρ1 is chosen as 108 cm−3 , 109 cm−3 , and 1010 cm−3 . From their Figure 4, for µ = 3 and5, and for ρ1 = 109 cm−3 and ρ1 = 1010 cm−3 , the specific angular momentum j(ρ > ρ1 ) of theperpendicular rotator (θ = 90◦ ) is larger than that of the aligned rotator (θ = 0◦ ). See Figure 21(A).
They concluded that the perpendicular rotator (θ = 90◦ ) looses its angular momentum slower than
the aligned rotator (θ = 0◦ ), which contradicts with the conclusion of Matsumoto & Tomisaka (2004).
The contradiction is now solved by paying attention to the value of the thermal to gravitationalenergy ratio, α0 (Tsukamoto et al. 2018, see also Hirano et al. 2020). Tsukamoto et al. (2018) compared
their models with α0 = 0.4 similar to Matsumoto & Tomisaka (2014) and that of α0 = 0.2 similar to
Joos et al. (2012). When α0 = 0.2, j(θ = 90◦ ) is larger than j(θ = 0◦ ) until ρc reaches 10−10 g cm−3
(Fig. 21 (left)). In contrast, in the model of α0 = 0.4, j(θ = 90◦ ) < j(θ = 0◦ ) until the end of thecalculation (Fig. 21 (right)). The cloud with small α0 contracts faster than that with large α0 . Thus,
the duration of the isothermal contraction in which the magnetic braking is working is longer forlarger α0 .
Beside the effect of α0 , they found a mechanism that weakens the magnetic braking in the modelwith θ = 90◦ and strengthens that in θ = 0◦ . This selective accretion mechanism is as follows:θ = 0◦ : In the model with θ = 0◦ , gas flows mainly in the z-direction and gas initially distributednear the z-axis accretes and form a protostar. Since such gas has slight angular momentum,
angular momentum of the pseudo-disk is thus small (strengthens the magnetic braking).
θ = 90◦ : In the model with θ = 90◦ , gas initially distributed far from the z-axis accretes. Since suchgas has relatively large angular momentum, angular momentum of the pseudo-disk is thus large(weakens the magnetic braking).

35

(A)

(B)

Figure 21: (A:) Evolution of specific angular momentum of the central portion. Top two rows aretaken from Joo et al. (2012). Their parameters are α0 = 0.25, β0 = 0.03, µ = 2 (top) andµ = 5 (middle). Left, middle, and right rows plot j(n > 1010 cm−3 ), j(n > 109 cm−3 ), andj(n > 108 cm−3 ) as jc . Perpendicular rotation model θ = 90◦ (red solid line) has larger
jc than the parallel model θ = 0◦ (blue dotted line). (B :) Bottom rows are taken from
Fig. 10 of Tsukamoto et al. (2018). jc (ρ > 10−12 g cm−3 ) is plotted against ρc . Model with
α0 = 0.2 (left) and that of α0 = 0.4 (right:) are compared. In each α0 , 3 models of θ areshown: θ = 0◦ (aligned rotator; solid line), θ = 45◦ (misaligned rotator; dashed line), andθ = 90◦ (perpendicular rotator; dotted line).

36

Figure 22: Outflow and disks at the epoch of 5000 yr after the protostar formation. Model parameters
are identical in all models. Initial B runs in z-direction in all models but the angularmomentum vector J is inclined at an angle θ from B. Panels (a), (b), (c), and (d )correspond to the model of θ = 0◦ , 10◦ , 45◦ , and 90◦ . Box sizes are taken as 3000 AU(left), 400 AU (middle), and 100 AU (right). A disk is displayed as face-on view in theright box. Taken from Fig. 4 of Hirano et al. (2020).

37

Figure 23: Outflows and disk with different (α0 , θ). Comparison with models different α0 clearlyshows that angular momentum brought to the central region increases with decreasing α0 .
This is understood the magnetic braking is effective in larger α0 . Taken from Fig. 10 of
Hirano et al. (2020).

As a result, the angular momentum jc is affected by the duration of magnetic braking and the selective
accretion mechanism.
The former mechanism mainly works during the isothermal phase before the protostar formation.
The latter works even in the adiabatic core stage after the protostar formation. How does themisalignment work in the later evolution? As accretion proceeds onto the protostar, the selective
accretion mechanism seems more important in the later phase. Figure 22 shows the outflow expectedfor θ = 0◦ , 10◦ , 45◦ , and 90◦ (Hirano et al. 2010). Initial condition of Figure 22 is relatively compact
compared with the previous studies: nc = 106 cm−3 , Mcl = 1 M⊙ , Rcl = 5900 AU, B0 = 430µG,Ωcl = 4.73 × 10−6 rad s−1 , which corresponds to α0 = 0.4, β0 = 0.02, and µ0 = 3. Owing to the sinkcell method, their simulation can reach ≃ 5000 yr after the core formation.
Figure 22 shows that a prominent outflow is observed in (a) and (b), that is, the angle between B
and J is θ ≲ 10◦ . Outflow becomes weaker with increasing θ. This seems to show that the componentof magnetic field parallel to the angular momentum is necessary to drive a prominent outflow. In
contrast, the rightmost box shows that the radius of the disk increases with increasing θ. This isunderstood owing to the selective accretion mechanism; specific angular momentum of accreting
gas increases as θ increases up to 90◦ .
Effect of α0 is explored by comparison of models with different α0 . Figure 23 plots the outflow and disk at the epoch of 5000 yr after the core formation. Models of Figure 23 have different(α0 , θ) but the same β0 = 0.02, µ = 3, and Rcl = 6200 AU. Thus B0 , Ωcl , and Mcl are chosen as: (B0 , Ωcl , Mcl ) = (970 µG, 6.63 × 10−6 yr, 2.5 M⊙ ) for models α0 = 0.2, (B0 , Ωcl , Mcl ) =(480 µG, 4.73×10−6 yr, 1.2 M⊙ ) for models α0 = 0.4, and (B0 , Ωcl , Mcl ) = (260 µG, 3.79×10−6 yr, 0.82 M⊙ )for models α0 = 0.6.
The effect of the angle θ is seen similar to α0 = 0.4 (Fig. 22) even for α0 = 0.2 and α0 = 0.6.

38

That is, the smaller θ, the more prominent the outflow; the larger θ, the larger the disk. Comparing
three columns, models α0 = 0.2 have a larger outflow in θ = 0 case and a larger disk in θ = 90◦
case. On the other hand, models α0 = 0.6 have a smaller outflow in θ = 0 case and a smaller disk inθ = 90◦ case. This variation is understood if we pay attention to the duration of magnetic braking
in the isothermal collapse phase is longer in α0 = 0.6 models compared with those of α0 = 0.2. This
indicates that a long timescale evolution is affected by (α0 , θ) and that observation of outflows anddisks can specify the initial condition of the star formation.

39

AppendixA

Maxwell Equation

Appearance of Maxwell equation depends on the unit adopted. In SI (MKSA) unit, the Maxwell
equation in the vacuum is written asρ,ϵ0∇ · B = 0,∇·E =∇×E =−

∂B
,
∂t


(A.1)

∂E
J + ϵ0∂t∇ × B = µ0,where we use notations of E (the electric field), B (the magnetic field), ρ (the electric charge density),j (the electric current density), ϵ0 (the permittivity), µ0 (the permeability), and c (the speed of light).
While, using the cgs Gauss unit, the Maxwell equation is written as∇ · E = 4πρ,∇ · B = 0,
∇×E =−
∇×B =

1
c

1 ∂B
,
c ∂t
4πJ +

In this document, we use the latter cgs Gauss unit.

40

(A.2)
∂E
∂t


.

B

Fluid Dynamics

B.1

What is ﬂuid?

m
V
ρ

m
V’
ρ’

t

t+∆t

Figure 24: Movement of fluid element

Gas and liquid change their shape according to the shape of the container. This is a definition of
ﬂuid.

B.2

Equation of Motion

Consider a fluid element with mass m as in Figure 24. From the Newton’s second law of mechanics,
the acceleration of an element is written as
dv
F= fm ≡ ,
dt
m

(B.1)

where m, v, F , and f m represent the mass, the velocity of the element, and force working on theelement and the force per mass working on the element, respectively.
In fluid dynamics, using the mass density ρ and the force working on the unit volume f v equation(B.1) is rewritten as
dv
(B.2)
ρ
= f v.
dt
Which kind of force works in a fluid? Pressure force does work in any fluids. Beside this, if thereis the gravity, ρg should be included in f v . If the electric currents is running in the fluid and themagnetic field exists, the Lorentz force j × B should be added.
To write down the expression of the pressure force, consider a fluid element between x and x + ∆x
as shown in Figure 25. Pressure force exerting on the surface S at x is p(x)S, while that on the
opposite side is −p(x + ∆x)S. The net pressure force working on the volume of S∆x is equal to(p(x) − p(x + ∆x))S, which is approximated as −(∂p/∂x)∆xS + O(∆x2 ). Thus, the pressure force
working on the unit volume is written −∂p/∂x. Equation (B.2) can be rewritten asρdv= −∇p + ρg,
dt

(B.3)

when gravity g is working.

B.3

Lagrangian and Eulerian Equations

The time derivative appearing in equation (B.1) expresses how the velocity of a speciﬁc particle
changes. Therefore, what appears in equation (B.3) means the same, that is, the position of the gas
element concerning in equation (B.3) moves and the positions at t and t + ∆t are generally different.
41

S
p(x)

p(x+∆x)x

x+∆x

Figure 25: Pressure force exerted on an element with volume S∆x

However, considering the velocity field in the space, the time derivative of the velocity should be
calculated staying at a fixed point x.
These two time derivative are different each other and should be distinguished. The former time
derivative is called Lagrangian time derivative and is expressed using d/dt. On the other hand,
the latter is called Eulerian time derivative and is expressed using ∂/∂t. These two are relatedwith each other. Consider a function F whose independent variables are time t and position x, that is
F (x, t). The difference dF
dt ∆t, using the Lagrangian time derivative of F , represents the the differenceof F (t + ∆t) from F (t) focusing on a specific fluid element, whose positions are different owing to its
motion. The element at the position of x0 at the epoch t0 moves to x0 + v 0 ∆t in time span of ∆t.
Thus the difference is expressed asdF∆t = F (x0 + v 0 ∆t, t0 + ∆t) − F (x0 , t0 ),
dt




∂F
∂F
≃
· v 0 ∆t +∆t,
∂x t x0 ,t0∂t x x0 ,t0(B.4)

where we used the Taylor expansion of F . The difference corresponding to the Eulerian derivative iswritten down as∂F∆t = F (x0 , t0 + ∆t) − F (x0 , t0 ),(B.5)∂t
and this is equal to the second term of the rhs of equation (B.4). Comparing equations (B.4) and(B.5), the Lagrangian derivative contains an extra term besides the Eulerian derivative. That is, theLagrangian derivative is expressed by the Eulerian derivative as
d
∂
∂
∂
=
+v·
=
+ v · grad.
dt
∂t
∂x
∂t

(B.6)

Applying the above expression on equation of motion based on the Lagrangian derivative (B.3),
we obtain the Eulerian equation motion:
∂v+ v · ∇v = −∇p + ρg.
(B.7)
ρ
∂t

B.4

Continuity Equation

Another basic equation comes from the mass conservation. This is often called the continuity equation,
which relates the change of the volume to its density. Consider a fluid element whose volume is equal42to ∆V . The mass contained in the volume is constant. Thus
dρ∆V
dρ
d∆V
=
∆V +
ρ = 0.
dt
dt
dt
The variation of the volume d∆Vdt is rewritten as
Z
Z
d∆V
=
v · dS =divvdV,
dt
∂∆V
∆V

(B.8)

(B.9)

where ∂∆V represents the surface of the fluid element ∆V . From equations (B.8) and (B.9), we
obtain the mass continuity equation for Lagrangian time derivative asdρ+ ρdivv = 0.
dt

(B.10)

Using equation (B.6) this is rewritten to Eulerian form as∂ρ+ div(ρv) = 0.
∂t

(B.11)

Basic equations using the Lagrangian derivative are equations (B.3) and (B.10), while those of
the Eulerian derivative are equations (B.7) and (B.11).

B.5

Expression for Momentum Density

It is useful to describe the equation for momentum ρv. Eulerian derivative of momentum density ρv
is rewritten as
∂ρv
∂t

∂ρ
∂v
v+ρ ,
∂t
∂t
= −∇(ρv)v − ρv · ∇v − ∇p + ρg,=(B.12)

where we used equations (B.7) and (B.11). In the Cartesian coordinate (x, y, z) = (xi , x2 , x3 ), thisbecomes

 

 




ρvx
∂/∂x
ρvx vx ρvx vy ρvx vz
∂/∂x
p 0 0
gx
∂ 
ρvy  +  ∂/∂y   ρvy vx ρvy vy ρvy vz  +  ∂/∂y   0 p 0  = ρ  gy  .
∂t
ρvz
∂/∂z
ρvz vx ρvz vy ρvz vz
∂/∂z
0 0 p
gz
(B.13)
Using the momentum stress tensor
πi j = ρvi vj + δi j p,
(B.14)
the above equation is written as∂ρvi ∂πi j= ρgi ,+
∂t∂xj

(B.15)

where δi j represents the Kronecker’s delta as δi j = 0 for i ̸= j and δi j = 1 for i = j.
If there is no external force, rhs of equation (B.15) is equal to zero. If the momentum density ρvwere to obey the continuity equation, equation (B.15) would be∂ρvi ∂ρvi vj
+
= 0.
∂t
∂xj

(B.16)

However, this is incorrect, because there exists the pressure force in the fluid and thus the momentumof the fluid element is not conserved.
43

B.6

Energy Equation

The above basic equations (B.3) and (B.10) or equations (B.7) and (B.11) contain three dependent
variables ρ, p, and v. The number of the variables, 3, is larger than the number of equations, 2.
Therefore, an extra equation is needed to close the basic equations.
B.6.1

Polytropic Relation

If the pressure of the fluid, p, is expressed only by the density, ρ,
p = P (ρ),(B.17)

the number of dependent variables is reduced to two and the above equations are sufficient to describe
the dynamics of the fluid. Occasionally, the pressure is assumed proportional to the power of ρ as
p = KρΓ ,(B.18)

where Γ is a constant. This assumption is called polytropic relation.
The fact should be reminded that the validity of the assumption comes from the physical conditionof the system. In the case that the temperature of the gas is kept constant owing to the cooling andheating process, the gas pressure is proportional to the density
p = c2is ρ,(B.19)

where cis = (kT /µmp )1/2 (k: Boltzmann constant, µ average molecular weight, and mp is the protonmass) represents the isothermal sound speed and is constant.
Another example is the isentropic fluid, in which the entropy is kept constant. In this case thepressure is proportional to ργ , as γ
ρp = p0,(B.20)
ρ0where γ is the specific heat ratio = cp /cv . In these cases the polytropic relations such as equations(B.19) and (B.20) play a role as the third equation of basic equations of hydrodynamics.
B.6.2

Energy Equation from the First Law of Theromodynamics

In more general cases, the last equation comes from the first law of the thermal physics,∆Q = dU + pdV.

(B.21)

where U and V represent the internal energy and the volume and ∆Q is the heat added to or
subtracted from the system. Using the total energy per unit volume,1U
ϵ = ρv 2 + ,2V(B.22)

the above equation gives the equation for the total energy per unit volume:∂ϵ+ div(ϵ + p)v = ρv · g.
∂t

(B.23)

Since the total energy per unit volume is expressed using the basic physical quantities as11
ϵ = ρv 2 +
p
2
γ−1
44

(B.24)

for the ideal gas, equation (B.23) is the final basic equation for hydrodynamics. Equations (B.7),(B.11), and (B.23) are basic equations of hydrodynamics using the Eulerian time derivative.
Problem 2
From the continuity equation (B.11) and the equation of motion (B.7) or the equation for momentum density (B.15), show that the equation for the kinetic energy density is expressed asfollows∂ 12 ρv 2
1 2+∇·
ρv v = ρg · v.
(B.25)∂t2
Then, obtain equation for the total energy (B.23).




45

C

Polytrope
1
0.8
0.6
0.4
0.2

1

2

4

3

5

7

6

Figure 26: Solutions of eq. (C.6) are plotted for indices of n = 3 or Γ = 4/3 (the most extended one),n = 5/2 or Γ = 7/5 (the middle), and n = 3/2 or Γ = 5/3 (the most compact one).
Although, these three solutions have zero-point at finite ξ, for n ≥ 5 gas extends infinitely
and such a solution has no zero-point.

If we choose the polytropic equation of state11 ,p = Kρ1+1/n ,

(C.1)

 2 
r dp
= −4πGρ,
ρ dr

(C.2)

the hydrostatic balance is expressed by
1 d
r2 dr

where we used equations (D.1) and (D.2). A hydrostatic gaseous star composed with a polytropic
gas is called polytrope. Normalizing the density, pressure and radius asρ ≡ ρc θ n ,
p ≡

(C.3)

Kρ1+1/nθn+1 ,c(1−n)/nr ≡ Hξ =(n + 1)Kρc
4πG

(C.4)

!1/2
ξ,

(C.5)

we obtain a normalized equation as
1 d
ξ 2 dξ


ξ

2 dθ

dξ


= −θn ,

(C.6)

which is called Lane-Emden equation of index n. The boundary condition at the center of polytrope
which is regular at the center (r = 0) should beθ = 1,dθ= 0,dξat ξ = 0. Solution of this equation is plotted for several n in Figure 26.
11

Reference book for this appendix is chapter 4 of Chandrasekhar (1939).

46

(C.7)
(C.8)

Mass of the polytrope is written down asZ r
M =4πρr2 dr,
0
Z ξ
3
= 4πH
θn ξ 2 dξ,0Z ξd 2 dθ
= 4πH 3
ξ
dξ,
dξ
0 dξ
dθ
= −4πH 3 ξ 2
dξ ξ1


(n + 1)K 3/2 (3−n)/2n 2 dθ
= −4π
ρc
ξ
,
4πG
dξ ξ1

(C.9)

where ξ1 represents the zero point of θ or the surface radius normalized by H. For n = 3 or
Γ ≡ 1 + 1/n = 4/3, equation (C.9) reduces to
M = −4π

K
πG

3/2
ξ2

dθ
.
dξ ξ1

(C.10)

Thus, the mass does not depend on the central density for Γ = 4/3 polytrope. For n = 5/2 orΓ = 1 + 1/n = 7/5, M is written down as
M = −4π

7K
8πG

3/2
ρc1/10 ξ 2

dθ
,
dξ ξ1

(C.11)

where ξ1 = 5.35528 and −ξ12 dθ/dξ|ξ1 = 2.18720. Polytrope with Γ = 7/5 gas, the mass-density
relation becomes ρc ∝ M 10 . While, for n = 3/2 of Γ = 1 + 1/n = 5/3, M is written down as
M = −4π

5K
8πG

3/2

where ξ1 = 3.65375 and −ξ12 dθ/dξ|ξ1 = 2.71406.

47

ρc1/2 ξ 2

dθ
,
dξ ξ1

(C.12)

D

Bonnor Ebert Mass

D.1

Hydrostatic Equilibrium in Tenuous Gas

Consider hydrostatic balance of a spherical isothermal cloud. By the gas density, ρ, the isothermalsound speed, cis , and the gravitational potential, Φ, the force balance is written as−

c2is dρ dΦ
−
= 0,
ρ dr
dr

and the gravity is calculated from a density distribution asZdΦGMr
4πG r 2−=− 2 =− 2
ρr dr,
dr
r
r
0

(D.1)

(D.2)

for a spherical symmetric cloud, where Mr represents the mass contained inside the radius r. Gravity
for a cylindrical cloud with infinite length is expressed
Z
dΦ
2Gλr
4πG r
−
=−
=−
ρrdr,
(D.3)
dr
r
r
0
where λr represents the mass per unit length within a cylinder of radius being r (line mass).
For the spherical symmetric case, the equation becomes the Lane-Emden equation with the polytropic index of ∞. This has no analytic solutions. However, the numerical integration gives us
a solution shown in Figure 27 (left). Only in a limiting case with the infinite central density, thesolution is analytically given asc2
ρ(r) = is r−2 .
(D.4)
2πG
Increasing the central density, the solution reaches the above singular isothermal solution (SIS).
On the other hand, a cylindrical cloud has an analytic solution (Stodólkiewicz 1963; Ostriker1964) as−2r2
ρ(r) = ρc 1 +
,
(D.5)
8H 2
where
H 2 = c2is /4πGρc .

(D.6)

Far from the cloud symmetric axis, the distribution of equation (D.5) gives

while the spherical symmetric cloud has

ρ(r) ∝ r−4 ,(D.7)

ρ(r) ∝ r−2(D.8)

distribution.
Problem 3Show that the SIS is a solution of equations (D.1) and (D.2).

Problem 4Show that the density distribution of equation (D.5) is a solution of equations (D.1) and (D.3).




48

1

1
’fort.4’’1.1.d’’1.d’’0.9.d’

0.1

0.1

0.01

0.01

0.001

0.001

0.00010.0001

1e-05
0.01

0.1

1

10

100

1e-05
0.01

1000

0.1

1

10

100

1000

Figure 27: Radial density distribution of a spherical cloud (left) and a cylindrical cloud (right). In theright panel, solutions for polytropic gases with Γ ≡ 1 + 1/n = 1.1 or n = 10 (a relatively
compact curve plotted with “+”) and Γ = 0.9 or n = −10 (relatively extended curve with“*”) are plotted as well as the isothermal one (middle curve with “×”).

D.2

Bonnor-Ebert Mass

In the preceding section [Fig. 27 (left)], we have seen the radial density distribution of a hydrostatic
configuration of the isothermal gas. Consider a circumstance that such kind of cloud is immersed
in a low-density medium with a pressure p0 . To establish the pressure equilibrium, the pressure
at the surface c2is ρ(R) must equal to p0 . This means that the density at the surface is constantρs = ρ(R) = p0 /c2is .
Figure 28 (left) shows three models of density distribution, ρc = ρ(r = 0) = 10ρs , 102 ρs , and3
10 ρs . Comparing these three models, it should be noticed that the cloud size (radius) decreases withincreasing the central density ρc . The mass of the cloud is obtained by integrating the distribution,
which is illustrated against the central-to-surface density√ ratio ρc /ρs in Figure 28 (right). The y-axis
represents a normalized mass as m = M/[4πρs (cis / 4πGρs )3 ]. The maximum value of m = 4.026
means
c2is
Mmax ≃ 1.14
.
(D.9)
1/2
G3/2 p0
This is the maximum mass which is supported against the self-gravity by the thermal pressure with anisothermal sound speed of cis , when the cloud is immersed in the pressure p0 . This is called BonnorEbert mass [Bonnor (1956), Ebert (1955)]. It is to be noticed that the critical state Mcl = Mmax isachieved when the density contrast is rather low ρc ≃ 16ρs ≡ ρcr .
Problem 5Show that the normalized form of equations (D.1) and (D.2) are represented as
1 df
f dξ

=

m
,
ξ2
Z ξ
f ξ 2 dξ,m =03
where nondimensional variables√ are defined as ξ ≡ r/H, f ≡ ρ/ρs , and m ≡ M/(4πρs H ) using
the scale-length of H = cis / 4πGρs .





49


Problem 6Numerically solve the above simultaneous differential equations and obtain the critical mass [eq.
(D.9)] and the central density of the critical state, ρcr . Boundary conditions of the above equations
are: at the center (ξ = 0), m = 0 and f = 1.
It should be noted that a gas cloud with Mcl < Mmax with isothermal sound speed cis andimmersed in the external pressure pext will have a mechanical balance between the pressure and theself-gravity (mechanical equilibrium). However, that with Mcl > Mmax is not supported against the
self-gravity. The former is called “subcritical” and the latter “supercritical.” Later, we will show
the critical mass is increased by rotation (centrifugal force) and magnetic field (Lorentz force).
Another important result from Figure 28 (right) is the stability of an isothermal cloud. Even fora subcritical cloud with Mcl < Mmax , any clouds on the part of (∂Mcl /∂ρc )p0 < 0 are unstable, whoseclouds are distributed on the branch with ρc > ρcr . This is understood as follows: For a hydrostatic
cloud, the mass should be expressed with the external pressure and the central density [Fig.28 (right)]as
Mcl = Mcl (p0 , ρc ).
(D.10)In this case, a relation between the partial derivatives such as






∂p0
∂ρc
∂Mcl
·
·
= −1,
∂p0 ρc∂ρc Mcl∂Mcl p0(D.11)

is satisfied, unless each term is equal to zero. Figure 27 (left) shows that the cloud mass Mcl is adecreasing function of the external pressure p0 = ρs c2is , if the central density is fixed. Since this means


∂Mcl
< 0,
(D.12)∂p0 ρcequation (D.11) gives us



∂p0
∂ρc




·
Mcl

∂ρc
∂Mcl


> 0.

(D.13)

p0

For a cloud with ρc < ρcr ≃ 16p0 /c2is , mass is an increasing function of the central density as


∂Mcl
> 0.
(D.14)∂ρc p0Thus, equation (D.13) leads to the relation


∂ρc
∂p0


> 0,

(D.15)

Mclfor ρ < ρcr . This means that gas cloud contracts (the central density and thus pressure increase),
when the external pressure increases. This is an ordinary reaction of a stable gas.
On the other hand, the cloud on the part of (∂Mcl /∂ρc )p0 < 0 (for 16 ≲ ρc /(p0 /c2is ) ≲ 2000)
behaves


∂ρc
< 0,
(D.16)∂p0 Mcland this represents that an extra external pressure decreases the central density and thus the pressure.
This means the pressure encourages expansion. This reaction is unstable.
50

1000

5
’F1000.d’’F100.d’’F10.d’’rhoc-M.d’

4

100

M

rho

3

2
10

1

1
0.001

0
0.01

0.1
x

1

10

1

10

100
rho

1000

10000

Figure 28: (Left:) radial density distribution. Each curve has√different ρc . The x-axis denotes theradial distance normalized by a scale-length as c2is / 4πGp0 . It is shown that the radiusincreases with decreasing ρc in this range (ρc = (10, 102 , 103 ) × p0 /c2is ). (Right:) the
relation between mass (Mcl ) and the central density (ρc ) is plotted, under the conditionof constant external pressure. The x-axis represents the central density √normalized byρs ≡ p0 /c2is . The y-axis represents the cloud mass normalized by 4πρs (cis / 4πGρs )3 .

51

E

Free-fall TimeEven if the pressure force can be neglected in the equation of motion (B.1), the gravitational force
remains. Assuming the spherical symmetry, consider the gravity gr (r) at the position of radial
distance from the center being equal to r. Using the Gauss’ theorem, gr is related to the mass
inside of r, which is expressed by the equation
Z r
Mr =ρ4πr2 dr,(E.1)
0

and gr is written as

GMr
.
r2
This leads to the equation motion for a cold gas under a control of the self-gravity is written
gr (r) = −

(E.2)

d2 r
GMr
=− 2 .
dt2
r

(E.3)

Analyzing the equation is straightforward, multiplying dr/dt gives
dv 2 /2
d GMr
=+
,
dt
dt
r

(E.4)

which leads to the conservation of mechanical energy as 
1 dr 2 GMr
−
= E,
2 dt
r

(E.5)

in which E represents the total energy of the pressureless gas element and it is obtained from theinitial condition. If the gas is static initially at the distance R, the total energy is negative asE=−

GMr (R),R(E.6)

because at t = 0, r = R and dr/dt = 0.
The solutions of equation (E.5) are well known as follows:
1. the case of negative energy E < 0. Considering the case that the gas sphere is inflowing v < 0,equation (E.5) becomes


dr
1 1/2
1/2 1
= − [2GMr (R)]
−
,
(E.7)
dt
r R
where we assumed initially dr/dt = 0 at r = R. Using a parameter η(t), the radius of the gas
element at some epoch t is writtenr = R cos2 η.
(E.8)
In this case, equation (E.7) reduces tocos2 η
This gives us the expression of t as

t=

dη
=
dt



R3
2GMr (R)

GMr (R)
2R3

1/2 
η+

1/2.

sin 2η
2

(E.9)


.

This corresponds to the closed universe in the cosmic expansion (Ω0 > 1).
52(E.10)

1
0.8
0.6
0.4
0.2
0.25

0.5

1

0.75

1.25

1.5

Figure 29: Free-fall. The x- and y-axes represent η + sin 2η/2 and cos2 η, respectively.

2. if the energy is equal to zero, the solution of equation (E.5) is written as2/3  9GM (R) 1/3
r
3/2
3/2
r −R
=
t2/3 ,
2(E.11)

where R = r(t = 0). This corresponds to the Einstein - de Sitter Universe with Ω0 = 1.
3. If E > 0, the expansion law is given by!−1 
 2 1/2
RṘ2 R
sinh 2ηt=−η
(E.12)
2E
2GMr (R)
2


and
x=

E
GMr (R)r = sinh2 η,(E.13)

where E represents the total energy
E=Ṙ2 GMr (R)
−
> 0.
2
R

(E.14)


Problem 7Solve equation (E.5) for E = 0 and obtain (E.11).
Problem 8Solve equation (E.5) for E > 0 and obtain equations (E.12) and (E.13).


In the present case, at t = 0, since dr/dt = 0 the energy is negative. Equation (E.8) shows us
r becomes equal to zero (the gas collapses) if η = π/2 as well as η = 0 at t = 0. Equation (E.10)
indicates it occurs at the epoch of
1/2

R3
π
t = tﬀ =,2GMr (R)2
3π 1/2=,(E.15)
32Gρ̄where ρ̄ represents the average density inside of r, that is Mr /(4πr3 /3). This is called “free-falltime” of the gas. This gives the time-scale for the gas with density ρ̄ to collapse. In the actualinterstellar space, the gas pressure is not negligible. However, tﬀ gives a typical time-scale for a gas
cloud to collapse and to form stars in it.
53

F

Alfvén Wave

Figure 30: Alfvén wave.

Consider uniform (ρ = ρ0 = const) gas contains uniform B-field B = t (B0 , 0, 0). We assume gas
is cold p = 0 and the gas has initially no motion v0 = 0. Basic equation of magnetohydrodynamics
is as follows:∂ρ+ ∇ · (ρv) = 0,
(F.1)
∂t


∂v
1
1
ρ
+ v · ∇v = − ∇B 2 +(B · ∇)B,
(F.2)
∂t
8π
4π
∂B
= ∇ × (v × B).
(F.3)∂t
Consider an incompressible wave δρ = 0 which coexists with transverse component velocity δv =t (0, δv, 0) and B-field δB = t (0, δB , 0). Equations for perturbed variables are from (F.2) and (F.3)
y
as follows:
∂δBy
∂δv
1
ρ0
=
B0
,
(F.4)
∂t
4π
∂x
∂δBy
∂
=
δvy B0 .
(F.5)
∂t
∂x
Eliminating δv with use of ∂/∂x of equation (F.4) and ∂/∂t of equation (F.5),∂ 2 δBy
B02 ∂ 2 δBy=.
2∂t
4πρ0 ∂x2(F.6)

This is satisfied by a wave equation as δBy = A exp(iωt − ikx). This mode is called Alfvén wave.
The phase velocity of the wave isvphase ≡

ω
B0
=±
.
k
(4πρ0 )1/2

This phase speed is called as Alfvén speed and is expressed with VA .

54

(F.7)

G

Magnetic Braking
(b)

(a)
ρext
B0
ρext
B0

z

Ωcl
Σcl

Ωcl

2h

z=0
R

Figure 31: Schematic model of magnetic braking. (a:) Aligned rotator (Ωcl //B0 ). A rotating cloudwith the column density Σcl and angular speed of Ωcl is considered. Torsional Alfvénwave propagating in the intercloud medium with the magnetic field strength of B0 andthe density of ρext . (b:) Perpendicular rotator (Ωcl ⊥ B0 ). A cloud with radius R andhalf height h is immersed in the intercloud medium. Torsional Alfvén wave is assumed topropagate radially only in a restricted region of −h < z < h. Cloud-to-intercloud density
ratio ρcl /ρext controls which rotation is braked more rapidly.

Angular momentum transfer or redistribution is possible in several ways: (for simplicity, we
consider 2D axisymmetric system, here)
Advection: Gas with angular momentum ρrvϕ moves with velocity v = (vr , vz ) brings angularmomentum flux density as ρrvϕ v.
Magnetic tension: Since the torque N = r × F changes angular momentum, even in the axisymmetric geometry, the ϕ-component of Lorentz force transfers the angular momentum.

G.1

Aligned Rotator

A simple mechanism to reduce the angular momentum of the cloud by magnetic tension is called
magnetic braking (Mouschovias & Paleologou 1980). Shown in Figure 31 (a), consider a cloud withthe column density Σcl being immersed in the intercloud medium with density ρext and magnetic fieldof B0 . The cloud is rotating with Ωcl but the intercloud gas has no angular momentum. Torsional
Alfvén wave propagates in the both vertical directions (z > 0 and z < 0). Wave fronts propagate with
the Alfvén speed as ±VA = B0 /(4πρext )1/2 . After passing the wave front |z| < VA t, the intercloud
medium begins corotation with the same angular speed of the cloud Ωcl . Since the inertia increases
with a rate of ρext VA , the angular momentum conservation is expressed asΣcl Ωcl (t = 0) = [2ρext VA t + Σcl ] Ωcl (t),(G.1)

which gives the time derivative of Ωcl as
1 dΩclΩcl dtΣcl 2ρext VA(2ρext VA t + Σcl )22ρext VA
∼ −
,
Σcl
= −

55

(G.2)

where we assumed 2ρext VA t ≪ Σcl to fit the early phase evolution. Thus, the rotation stopping
timescale is given as
Σcl
τ∼
.
(G.3)
2ρext VAIt should be noticed that the timescale is longer for lower intercloud gas density, although the prop−1/2
agating speed is increased as VA ∝ ρext .

G.2

Perpendicular Rotator

In contrast, when the cloud is rotating in the perpendicular direction to the magnetic field, Ω ⊥ B(Fig. 31 (b)), the moment of inertia of intercloud gas swept by the wave increases
Z
Iext =
r2 dm,
Z r=
r2 2hρext 2πrdr,0= ρext πhr4 ,(G.4)

where h represents the half thickness of the cloud and we assumed the Alfvén wave propagates in theintercloud restricted region of −h < z < h. After sweeping the external medium of r < rswept , Iext
becomes equal to that of the cloud
Icl = ρcl πhR4 ,(G.5)
where ρcl and R represent the cloud density and its radius. This gives
rswept = R

ρcl
ρext

1/4
.

(G.6)

In the aligned rotator, equation (G.3) gives the height where the Alfvén wave has reached asρclzswept = h
.
(G.7)
ρext
Comparing equations (G.6) and (G.7), rswept /zswept = (R/h)/(ρcl /ρext )3/4 ≪ 1, since ρcl ≫ ρext . The
perpendicular rotator seems to stop its rotation much faster than the aligned rotator (Mouschovias& Paleologou 1979).

56

H

Self-Similar Solution of Dynamical Contraction

The basic equations of spherically symmetric isothermal flow are

∂ρ
1 ∂
+ 2
r2 ρv = 0,
∂t
r ∂r
∂v
∂v c2is ∂ρ GM
+v
+
+ 2
= 0,
∂t
∂r
ρ ∂r
r
Z r
M (r, t) = M (0, t) +4πr′2 ρ(r′ , t)dr′ ,

(H.1)
(H.2)
(H.3)

0

where M (r, t) represents the mass included in the radius r and M (0, t) denotes the mass of the
protostar. A self-similar solution which has a formΩ(ξ),
4πG(t − t0 )2v(r, t) = cis V (ξ),
c3 |t − t0 |M (r, t) = is
m(ξ),
G
r
ξ =
,
cis |t − t0 |ρ(r, t) =

(H.4)
(H.5)
(H.6)
(H.7)

should be found, where Ω and V are functions only dependent on ξ. For example, equation (H.4)
asks the shape of the density distribution is the same after resizing of equation (H.7) r → r/cis |t − t0 |and re-normalizing ρ → ρ · 4πG(t − t0 )2 as Ω(ξ).

H.1

Derivation of Equation for Self-Similar Variables∗

Equations (H.4)-(H.7) give the relation between self-similar variables (independent variable, ξ and dependent
variables, Ω(ξ), V (ξ), m(ξ)). For t < t0 , equation (H.7) is written as ξ = −r/cs (t − t0 ), which gives
 
−1
∂ξ
=
(H.8)
∂r t
cs (t − t0 )
 
∂ξ
−r
ξ
=
=−
(H.9)
.
2
∂t r
cs (t − t0 )t − t0
When t − t0 < 0, ∂ξ/∂t > 0 means that a typical point of the solution (ξ = const) shrinks with time. The firstterm of the continuity equation (H.1) is rewritten as


∂
Ω(ξ)
∂ρ
=
,
∂t
∂t 4πG(t − t0 )2
1
dΩ ∂ξ
−2Ω
+
,
=
4πG(t − t0 )34πG(t − t0 )2 dξ ∂t


−1
dΩ
(H.10)=
2Ω + ξ
,3
4πG(t − t0 )dξand the second term is

1 ∂
r2 ρv
2
r ∂r

=
=
=
=

∂ρv
2ρv
+
,
∂r 
r

2Ωcs V∂Ωcs V +,∂r 4πG(t − t0 )2−cs (t − t0 )ξ4πG(t − t0 )2


d
−1
Ω
2Ωcs Vcs V,+
dξ 4πG(t − t0 )2
cs (t − t0 ) −cs (t − t0 )ξ4πG(t − t0 )2


dΩV
−1
2ΩV
+
.
4πG(t − t0 )3
dξ
ξ

57

(H.11)

Thus, we obtain the continuity equation for similarity variables as
dΩV
2ΩV
dΩ
+
+ξ
+ 2Ω = 0.
dξξdξ(H.12)

Equation (H.3) or ∂Mr /∂r = 4πr2 ρ is rewritten using self-similar variables as 3dcs (t − t0 )
dξ
Ω
2
−
m
= 4π [cs (t − t0 )ξ]
,
dξ
G
dr
4πG(t − t0 )2
−1
dm
c2s 2
c3 (t − t0 )
=
ξ Ω,
∴− s
G
cs (t − t0 ) dξGdm∴
= ξ 2 Ω.
dξOn the other hand, ∂Mr /∂t = −4πr2 ρv leads to 3
∂cs (t − t0 )
−
m
∂t
G
c3
c3 (t − t0 ) dm ∂ξ∴ − sm− s
G
G
dξ ∂t
dm
∴ −m + ξdξ2= 4π [cs (t − t0 )ξ] cs V
=

c3s 2
ξ VΩ
G

=

ξ 2 V Ω,

(H.13)Ω
4πG(t − t0 )2(H.14)

where we used equation (H.9). With help of equation (H.13), this is reduced tom = ξ 2 (ξ − V )Ω,(H.15)

by which equation (H.16) is derived.
Problem 9Referring to the above, derive the equation of motion using the self-similar variables from equation (H.2).
Then, show the other two equations (H.17) and (H.18).



H.2

Self-Similar Solution

After calculation as shown in the previous subsection, the basic equations for the spherical symmetricmodel yieldm = (ξ − V )ξ 2 Ω,



 dV
2
2
(ξ − V ) − 1= Ω(ξ − V ) −(ξ − V ),
dξ
ξ


2
(ξ − V )2 − 1 dΩ= Ω − (ξ − V ) (ξ − V ),Ωdξ
ξ

(H.16)(H.17)(H.18)[eqs. (2.9)-(2.14) of Whitworth & Summers 1985; Hanawa 1999]. Equations (H.17) and (H.18) havea singular point at which (ξ − V )2 − 1 = 0 or V = ξ ± 1. Since the point of ξ =const moves with cis ,the flow velocity relative to this ξ =const is equal to v − cis ξ = (V − ξ)cis . Thus the singular pointξ∗ at which V = ξ ± 1 corresponds to a sonic point. Therefore, since the flow has to pass the sonicpoint smoothly, the rhs of equations (H.17) and (H.18) have to be equal to zero at the singular pointξ∗ . This gives at the sonic point ±ξ∗ > 0,Ω−2(ξ − V )2
= Ω(ξ − V ) − = 0,ξξ58

(H.19)

which leads toV∗ = ξ∗ ∓ 1,

(H.20)

2.
ξ∗(H.21)Ω∗ = ±

These equations (H.16), (H.17) and (H.18) have an analytic solution
V (ξ) = 0, Ω =2, m = 2ξ, −∞ < ξ < ∞.
ξ2

(H.22)

This is a solution which agrees with the Chandrasekhar’s SIS [eq. (D.4)]. Generally, solutions are
obtained only by numerical integration. For ξ → ∞, the solution has to converge to an asymptotic
form ofA − 2 V∞ 4V∞ + (A − 2)(A − 6)+ 2 ++ O(ξ −4 ),ξξ6ξ 3Ω∞ (A − 2)A−
+ O(ξ −6 ).
2
ξ
2ξ 4

V (ξ) = V∞ −(H.23)Ω(ξ) =(H.24)

This shows that for sufficiently large radius the gas flows with a constant inflow velocity V∞ cis .
This has a solution in which the density and the infall velocity should be regular with reaching
the center (ξ ≪ 1). Such kind of solution is plotted in Figure 32 (solid lines represent the solutionsof τ = t − t0 < 0). This time evolution is expected from the self-similar solution. This shows that≃ ρc(in the central region),ρ(H.25)∝ r−2 (in the outer envelope),∝r(in the central region),v(H.26)≃ 3.28cis (in the outer envelope).
Reaching the outer boundary, the numerical solution (Fig. 12) differs from the self-similar solution(Fig. 32). For example, v is reduced to zero in the numerical simulations, while it reaches a finite value
3.28 in the self-similar solution. And as for the density distribution, ρ drops near the outer boundary
in the numerical simulations while it decreases proportional to ∝ r−2 . However, in the region except
for the vicinity of the outer boundary the self-similar solution expresses well the dynamical collapseof the spherical isothermal cloud. This solution gives the evolution of a pre-protostellar core formed
in a supercritical cloud/cloud core.

H.3

Inside-out Collapse Solution

In 1977, Shu found another self-similar solution which is realized after a central protostar with
infinitesimal mass is formed in the singular isothermal sphere (SIS) solution. The gas begins to accrete
to the protostar. Outside the region where the accretion occurs, the initial SIS is kept unchanged,
since the SIS is a hydrostatic solution. And the front of accretion expands radially outward in time.
Since the inflow region expand outwardly, he called it the inside-out collapse solution. In Figure
33, the evolution is shown for density and inflow velocity. This solution gives the evolution of a
protostellar core formed in a subcritical cloud/cloud core.
Using equations (H.23) and (H.24) and assuming the inflow velocity should reduce at large radius,
we obtain
A − 2 (A − 2)(A − 6)
V (ξ) = −−+ ··· ,
(H.27)
ξ
6ξ 3
59

16

3000

t=0y

14

+104
12 -10vr [m s-1]

log n [m-3]

2000
4

+105
10

-105

1000

+106
8

+104 +105 +106
t=0y

-106

-104
6
13

14

15
16
log r [m]

0
13

17

14

-105

-106

15
16
log r [m]17

Figure 32: A self-similar solution indicating a dynamical collapse of isothermal spherical cloud(Larson-Penston solution). Spatial distributions of the density and inflow velocity which
are expected from the self-similar solution are plotted. Dashed lines show the evolutionof prestellar core phase and solid lines show that of protostellar core phase. Attachednumbers represents time τ counted from the epoch (t0 ) when the central density reaches
infinity as τ = t − t0 . Taken from Hanawa (1999).

60

A(A − 2)A−+ ··· .
ξ2
2ξ 4
Since ξ → ∞ means t → t0 (if r is finite), Ω(ξ) → A/ξ 2 means thatΩ(ξ) =(H.28)

Ac2is.
(H.29)
4πGr2Comparing with the SIS, when A = 2 this gives the SIS and when A > 2 this gives a density
distribution in which the pressure is inefficient and the cloud is contracting. The solution with A > 2
is obtained by a procedure as (1) at a sufficiently large radius ξ1 , calculate V (ξ1 ) and Ω(ξ1 ). (2) fromthese values, integrate equations (H.17) and (H.18) inwardly. Figure 33 show the solution of thistype. The solution with A > 2 inflow speed is accelerated towards the center. Decreasing A (A → 2),
it is shown that an outer part ξ ≳ 1 reaches V → 0. For A = 2+ , the solution reaches the singularline V = −ξ + 1 at ξ = 1 (V = 0) 12 . Since V = 0 and Ω = 2 at ξ = 1, this solution with A = 2+
converges to the SIS at ξ = 1. This means that if there is an infinitesimally small amount of excess
mass at the center of SIS, the accretion begins from the center while outside a radius the cloud is left
static. The inner part of the solution ξ ≲ 1, V and Ω are well expressed as V ∝ ξ 1/2 and Ω ∝ ξ 3/2 .
ρ(r, t0 ) =Figure 33: Self-similar solution which shows the inside-out collapse. (Left:) Infall velocity (−V =−v/cis ) is plotted against the similarity variable ξ ≡ r/cis (t − t0 ). Three curves correspondto models A = 2.4, 2.1, and 2.01. (Right:) Density (Ω = 4πGρ(t − t0 )2 ) is plotted against
the similarity variable ξ ≡ r/cis (t − t0 ).

The power-law distributions of V ∝ ξ 1/2 and Ω ∝ ξ 3/2 are explained as follows: Conservation of
the total energy of an inflowing gas shell is expressed as
v 2 GMr
GMr0,(H.30)−=−
2rr0
where r0 denotes the initial radius of a gas shell and Mr0 = Mr represents the mass inside the gas
shell. Neglecting GMr0 /r0 compared with the term GMr /r (r0 ≫ r), this gives
2GMr 1/2∝ r−1/2 ,(H.31)v≃r
12This means A = 2 + ϵ where ϵ > 0 and ϵ → 061

where we assumed a major part of Mr comes from the mass of a protostar M∗ , that is,
Z r
M r = M∗ +
ρ4πr2 dr ≃ M∗ .

(H.32)

0

Since the average density inside the radius r is written asRrρ4πr2 drρ̄(< r) = R0 r
= 3ρ(r),
2
0 4πr dr

(H.33)

for SIS distribution, the time necessary for a gas shell to reach the center is proportional to ∝ ρ̄(<r)−1/2 ∝ ρ(r)−1/2 ∝ r1 ∝ (t − t0 )1 , where we used a fact that the front of accretion expands witha constant sound speed cis . This means that the time necessary for the gas shell to travel from theradius of the accretion wave front to the center is proportional to t − t0 . Since the mass of the shell
which begins accretion in a unit time is equal to ρ(r)4πr2 cis and is constant irrespective of t − t0 .
These two facts indicate that the mass accretion rate is constant in time. That is,Ṁ = 4πr2 ρv = const.

(H.34)

Since the mass swept by the sound wave per unit time is equal to Ṁsw = 4πr2 ρSIS cis = 2c3s /G, theaccretion rate is proportional to c3is /G. Shu (1977) obtained the accretion ratec3
Ṁ = 0.975 is ,G

(H.35)

for SIS A = 2+ . Using this equation, equations (H.31) and (H.34) indicate that the spatial density
distribution is expressed by a power-law as
ρ(r) ∝ r−3/2 ,which is valid for the region except for the vicinity of the front of the accretion.

62(H.36)

References[1] Alfvén, H. 1942, Existence of Electromagnetic-Hydrodynamic Waves, Nature, 150, 405[2] Blandford, R. D., & Payne, D. G. 1982, Hydromagnetic Flows from Accretion Discs and theProduction of Radio Jets, MNRAS, 199, 883[3] Bonnor, W. B. 1956, Boyle’s Law and Gravitational Instability, MNRAS, 116, 351[4] Chandrasekhar, S. & Fermi, E. 1953, Magnetic Fields in Spiral Arms, ApJ, 118, 113[5] Crutcher, R. M. 2004, Observations of Magnetic Fields in Molecular Clouds, in “The MagnetizedInterstellar Medium,” Proceedings of the conference, held in Antalya, Turkey, September 8 - 12,
2003, Eds: B. Uyaniker, W. Reich, and R. Wielebinski, Copernicus GmbH, Katlenburg-Lindau.,
p. 123-132.
[6] Crutcher, R. M. 2012, Magnetic Fields in Molecular Clouds, Annual Review Astronomy & Astrophys., 50, 29[7] Davis, L. 1951, The Strength of Interstellar Magnetic Fields, Physical Review, 81, 890[8] Ebert R. 1955, Über die Verdichtung von H i-Gebieten, Z. Astrophys., 37, 217[9] Hanawa, T. 1999, Star Formation, in “Active Universe” in Japanese (Shokabo, Tokyo) ed. by K.
Shibata, J. Fukue, R. Matsumoto, & S. Mineshige[10] Hirano, S., Tsukamoto, Y., Basu, S., & Machida, M. N. 2020, The Eﬀect of Misalignment between
the Rotation Axis and Magnetic Field on the Circumstellar Disk, ApJ, 898,118[11] Joos, M., Hennebelle, P., & Ciardi, A. 2012, Protostellar Disk Formation and Transport ofAngular Momentum during Magnetized Core Collapse, AAp, 543, 128[12] Larson, R. B. 1969, Numerical calculations of the dynamics of collapsing proto-star, MNRAS,
145, 271[13] Matsumoto, T. & Tomisaka, K. 2004, Directions of Outﬂows, Disks, Magnetic Fields, and Rotation of Young Stellar Objects in Collapsing Molecular Cloud Cores, ApJ, 616, 266[14] Mouschovias, T. Ch. & Paleologou, E. V. 1979, The Angular Momentum Problem and Magnetic
Braking: an Exact, Time-dependent Solution, ApJ, 230, 204[15] Mouschovias, T. Ch. & Paleologou, E. V. 1980, Magnetic Braking of an Aligned Rotator during
Star Formation - an Exact, Time-dependent Solution, ApJ, 237, 877[16] Mouschovias, T. Ch. & Spitzer, L., Jr. 1976, Note on the Collapse of Magnetic Interstellar
Clouds, ApJ, 210, 326[17] Masunaga, H. & Inutsuka, S. 2000, A Radiation Hydrodynamic Model for Protostellar Collapse.
II. The Second Collapse and the Birth of a Protostar, ApJ, 531, 350[18] Nakano, T. 1979, Quasistatic Contraction of Magnetic Protostars due to Magnetic Flux Leakage- Part One - Formulation and an Example, PASJ, 31, 713

63

[19] Nakano, T. Nishi, R., & Umebayashi, T. 2002, Mechanism of Magnetic Flux Loss in Molecular
Clouds, ApJ, 572, 199[20] Ogino, S., Tomisaka, K., & Nakamura, F. 1999, Gravitational Collapse of Spherical Interstellar
Clouds, PASJ, 51, 637.
[21] Paleologou, E. V., & Mouschovias, T. Ch. 1983, The Magnetic Flux Problem and Ambipolar
Diﬀusion during Star Formation - One-Dimensional Collapse. I - Formulation of the Problem andMethod of Solution, ApJ, 275, 838[22] Shu, F. H. 1978, Self-similar Collapse of Isothermal Spheres and Star Formation, ApJ, 214, 488[23] Shu, F. H. 1983, Ambipolar Diﬀusion in Self-Gravitating Isothermal Layers, ApJ, 273, 202[24] Stodólkiewicz, J. S. 1963, On the Gravitational Instability of Some Magneto-Hydrodynamical
Systems of Astrophysical Interest. Part III., AcA, 13, 30[25] Tomisaka, K. 1998, Collapse-Driven Outﬂow in Star-Forming Molecular Cores, ApJL, 502, L163
[26] Tomisaka, K. 2000, The Evolution of the Angular Momentum Distribution during Star Formation, ApJL, 528, L41[27] Tomisaka, K. 2002, Collapse of Rotating Magnetized Molecular Cloud Cores and Mass Outﬂows,
ApJ, 575, 306[28] Tomisaka, K., Ikeuchi, S., & Nakamura, T. 1990, The Equilibria and Evolutions of Magnetized,
Rotating, Isothermal Clouds. IV - Quasi-Static Evolution, ApJ, 362, 202[29] Tsukamoto, Y. et al. 2015, Eﬀects of Ohmic and Ambipolar Diﬀusion on Formation and Evolution of First Cores, Protostars, and Circumstellar Discs, MNRAS, 452, 278[30] Ustyugova, G. V., Koldoba, A. V., Romanova, M. M., Chechetkin, V. M., & Lovelace, R. V.
E. 1999 Magnetocentrifugally Driven Winds: Comparison of MHD Simulations with Theory, ApJ,
516, 221[31] Wardle, M. 2004, Star Formation and the Hall Eﬀect, ApSS, 292, 317[32] Wardle, M. 2007, Magnetic Fields in Protoplanetary Disks, ApSS, 311, 35[33] Whitworth, A. & Summers, D. 1985, Self-Similar Condensation of Spherically Symmetric SelfGravitating Isothermal Gas Clouds, MNRAS, 214, 1[34] Wurster, J. 2016, NICIL: A Stand Alone Library to Self-Consistently Calculate Non-Ideal Magnetohydrodynamic Coeﬃcients in Molecular Cloud Cores, PASP, 33, e041

64Index
Alfvén speed, 55
Alfvén wave, 11, 54, 55
aligned rotator, 32
ambipolar diffusion, 16, 19
angular momentum (specific), 31

magneto-centrifugal wind mechanism, 14
magnetocentrifugal wind, 13
Maxwell stress tensor, 2
misaligned rotator, 32
momentum flux density tensor, 2Bonnor-Ebert mass, 8, 49

nested grid method, 24
non-ideal magnetohydrodynamics , 4continuity equation, 42
critical mass, 50
Davis–Chandrasekhar-Fermi Method, 11
double-polytropic equation of state, 24, 33
drift velocity, 16
dynamical collapse, 22
equation of motion, 41
Eulerian time derivative, 42
Faraday’s induction equation, 3
first law of thermodynamics, 44
flux freezing, 4free-fall time, 52

Ohm’s law, 3
ohmic dissipation, 19
polytrope, 46
polytropic relation, 44
preprotostellar phase, 22
protostellar phase, 23
pseudo-disk, 26
rotational to gravitational energy ratio (β0 ), 33,
35
runaway collapse, 24
runway collapse, 23

I-tyle outflow, 30
ideal magnetohydrodynamics , 4
inside-out collapse, 59

selective accretion mechanism, 35
self-similar solution, 22
singular isothermal solution (SIS), 48, 59
singular isothermal sphere (SIS), 59
sink cells, 23
subcritical, 12
subcritical cloud, 50
supercritical, 12supercritical cloud, 50

Kepler speed, 14

thermal to gravitational energy ratio (α0 ), 33, 35

Lagrangian time derivative, 42
Lane-Emden equation, 46
Larson-Penston solution, 22
Lorentz force, 1, 8, 16

U-type outflow, 30

Gauss’s theorem, 4gravity flux tensor, 2Hall effect, 19
Hall parameter, 20

virial relation, 7

magnetic braking, 31, 35, 55
magnetic pressure, 1, 8
magnetic Reynolds number, 4magnetic tension, 1
magnetic to gravitational energy ratio (γ0 ), 34
magneto-centrifugal wind, 13, 28
65