Star Formation1
Kohji Tomisaka
National Astronomical Observatory Japan
Thursday 22nd November, 20121

http://th.nao.ac.jp/MEMBER/tomisaka/Lecture Notes/StarFormation/6.pdf

cCopyright 2001-2012
by Tomisaka, K.
version 0.5: August 5, 2001
version 1.0: August 30, 2001
version 2.0: November 8, 2003
version 3.0: June 1, 2004
version 4.0: Oct. 29, 2007
version 5.0: Dec. 9, 2009version 6.0: Oct. 3 – , 2012
Since this is currently under updating for a lecture of autumn semester of 2012, Dept. of Astronomical
Science, GUAS, be careful if you download and use. This was compiled at 11:25 of Thursday 22nd
November, 2012.

Contents
1 Introduction
1.1 Interstellar Matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Case Study — Taurus Molecular Clouds . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 T Tauri Stars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Spectral Energy Distribution (SED) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Protostars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5.1 B335 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5.2 L1551 IRS 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 L 1544: Pre-protostellar Cores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7 Rotation of Outﬂow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8 Magnetic Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8.1 Serpens-South Cloud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8.2 Prestellar Core . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8.3 Cores with Protostars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8.4 Magnetic Field ∼1000AU Scale . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.9 Density Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.10 Mass Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.11 Line Width - Size Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
3
5
7
10
10
14
17
17
20
21
21
23
23
27
30
32

2 Physical Background37
2.1 Basic Equations of Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.2 The Poisson Equation of the Self-Gravity . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.3 Free-fall Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.1 Accretion Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.4 Gravitational Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.4.1 Linear Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.4.2 Sound Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.5 Jeans Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.6 Gravitational Instability of Thin Disk . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.6.1 Rotating Disk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.7 Convective Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.8 Super- and Subsonic Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.8.1 Flow in the Laval Nozzle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.8.2 Steady State Flow under an Inﬂuence of External Fields . . . . . . . . . . . . . 49
2.8.3 Stellar Wind — Parker Wind Theory . . . . . . . . . . . . . . . . . . . . . . . 51
2.9 Virial Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
i

ii

CONTENTS
2.9.1 Virial Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.10 Radiative Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.10.1 Radiative Transfer Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.10.2 Einstein’s Coeﬃcients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.10.3 Relation of Einstein’s Coeﬃcients to Absorption and Emissivity . . . . . . . . .
2.10.4 Excitation Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.10.5 Critical Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56
57
57
58
59
61
61

3 Galactic Scale Star Formation
63
3.1 Schmidt Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.1.1 Global Star Formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.1.2 Local Star Formation Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2 Distributions of Atomic and Molecular Hydrogen . . . . . . . . . . . . . . . . . . . . . 67
3.3 Gravitational Instability of Rotating Thin Disk . . . . . . . . . . . . . . . . . . . . . . 69
3.3.1 Tightly Wound Spirals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.3.2 Toomre’s Q Value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.4 Spiral Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.5 Density Wave Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.5.1 Group Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.6 Galactic Shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4 Local Star Formation Process81
4.1 Hydrostatic Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.1.1 Bonnor-Ebert Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.1.2 Equilibria of Cylindrical Cloud . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Virial Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.2.1 Magnatohydrostatic Clouds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3 Subcritical Cloud vs Supercritical Cloud . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.3.1 Observation of Magnetic Field Strength . . . . . . . . . . . . . . . . . . . . . . 88
4.3.2 Observation of Mass-to-Flux Ratio . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.4 Ambipolar Diﬀusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.4.1 Ionization Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.4.2 Ambipolar Diﬀusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.5 Dynamical Collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.5.1 Inside-out Collapse Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.5.2 Protostellar Evolution of Supercritical Clouds . . . . . . . . . . . . . . . . . . . 99
4.6 Star Formation Time Scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.7 Accretion Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.8 Outﬂow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.8.1 Magneto-driven Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.8.2 Entrainment Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.9 Evolution to Star . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.10 Example of Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.11 Evolution in the H-R diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.11.1 Main Accretion Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.11.2 Premain-sequence Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

CONTENTS
Bibliography
A Basic Astrophysical Quantities

iii
121
129

B Basic Equation of Fluid Dynamics
131
B.1 What is ﬂuid? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
B.2 Equation of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
B.3 Lagrangian and Euler Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
B.4 Continuity Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
B.4.1 Expression for Momentum Density . . . . . . . . . . . . . . . . . . . . . . . . . 133
B.5 Energy Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
B.5.1 Polytropic Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
B.5.2 Energy Equation from the First Law of Theromodynamics . . . . . . . . . . . . 135
B.6 Shock Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
B.6.1 Rankine-Hugoniot Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
C Basic Equations of Magnetohydrodynamics
137
C.1 Magnetohydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
C.1.1 Flux Freezing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
C.1.2 Basic Equations of Ideal MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
C.1.3 Axisymmetric Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
D Hydrostatic Equilibrium
141
D.1 Polytrope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
D.2 Magnetohydrostatic Conﬁguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
E Basic Equations for Radiative Hydrodynamics
145
E.1 Radiative Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
F Random Velocity

147

G X-factor

149

H Abundance
151
H.1 Solar system abundance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
H.2 Isotope ratio (Isotopic abundance) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
H.3 Molecular Weight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

iv

CONTENTS

Chapter 1

Introduction
1.1

Interstellar Matter
8.0
7.0"Coronal" Gas

Stellar Wind

6.0

Supernovalog T (K)Star

Interaction with Supernova

5.0

Dissociation of H2

HII Regions

Second Collapse

4.0
Mass Loss from Red GiantIntercloud Gas
3.0
2.0

Photo-dissociation

Diffuse Cloud

1.0
0.0
-4.0

Warm or Hot Core

Quasi-static Contraction

Isothermal Collapse

Molecular Cloud
-2.0

0.0

2.0

First Stellar Core

Globule
4.0

6.0

8.0

10.0 12.0
log n (cm-3)

14.0

16.0

Second Stellar Core
18.0

20.0

22.0

Figure 1.1: Multiphases of the interstellar medium. The temperature and number density of gaseous
objects of the interstellar medium in our Galaxy are summarized. Originally made by Myers (1978),reconstructed by Saigo (2000).
Figure 1.1 shows the temperature and number density of gaseous objects in our Galaxy. Coldinterstellar medium forms molecular clouds (T ∼ 10K) and diﬀuse clouds (T ∼ 100K). Warm inter4
stellar medium 103 K<∼T <∼10 K are thought to be pervasive (wide-spread). HII regions are ionizedby the Ly continuum photons from the nearby early-type stars. There are coronal (hot but tenuous)gases with T ∼ 106 K in the Galaxy, which are heated by the shock fronts of supernova remnants.
4−3
Pressures of these gases are in the range of 102 K cm−3 <∼ p<∼10 K cm , except for the HII regions.
This may suggest that gases are in the pressure equilibrium.
1

2

CHAPTER 1. INTRODUCTION

In this ﬁgure, a theoretical path from the molecular cloud core to the star is also shown. We will
see the evolution more closely in Chapter 4.

Figure 1.2: (left) Temperature-density relation of interstellar gas in which thermal equilibrium isachieved between two eﬀects of heating and cooling. Above the line, cooling due to radiative loss isdominated, while below it the heating is dominated over the cooling. (right) The same as the left one,but for pressure-density relation. Figures are taken from Figs. 1.15 and 1.16 of Hayashi et al.(1978)(see also Draine 2011).
From theoretical point of view, in Figure 1.2, we plot the temperature-density (left) and pressuredensity (right) relations for a gas in which the thermal equilibrium is achieved between two eﬀects−13−2
to heat the gas and to cool the gas. If the interstellar pressure p is higher than p>∼10 dyn cm ,−22
−3
Figure 1.2 (right) indicates that interstellar gas density is higher than ρ>∼10 cm . On the other−14 dyn cm−2 , low-density interstellar matter
hand, if the interstellar pressure p is lower than p<
10∼−26 cm−3 . Between these two pressures, there are three states coexisting in the interstellar
exists ρ<
10∼
space. Since the middle state is unstable (because increasing the temperature leads to increasing theheating rate or decreasing the cooling rate), there exist two stable phases in pressure equilibrium.
This is a picture of “two-phase interstellar medium” (Field et al. 1969; Spitzer 1978). Intercloud gas
and diﬀuse cloud gas are explained by this mechanism. Historically, the discovery of hot coronal gas
5of T >∼several × 10 K stimulated an idea that the interstellar gas is globally divided into three phases(Cox & Smith 1974; McKee & Ostriker 1977).
Globally, the molecular form of Hydrogen H2 is much abundant inside the Solar circle, while theatomic hydrogen HI is more abundant than molecular H2 in the outer Galaxy. In Figure 1.3 (left),the radial distributions of molecular and atomic gases are shown. The right panel shows similar
distributions for four typical external galaxies (M51, M101, NGC6946, and IC342). This indicates1.2. CASE STUDY — TAURUS MOLECULAR CLOUDS

3

these distributions are similar with each other. HI is distributed uniformly, while H2 density increases
greatly reaching the galaxy center. In other words, only in the region where the total (HI+H2 ) density
exceeds some critical value, H2 molecules are distributed.

Figure 1.3: Radial distribution of H2 (solid line) and HI (dashed line) gas density. (Left:) our Galaxy.
Converting from CO antenna temperature to H2 column deity, n(CO)/n(H2 ) = 6 × 10−5 is assumed.
Taken from Gordon & Burton (1976). (Right:) Radial distribution of H2 and HI gas for external
galaxies. The conversion factor is assumed constant X(H2 /CO) = 3 × 1020 H2 /K km s−1 . Taken from
Honma et al (1995). About the X-factor refer to appendix G.

1.2

Case Study — Taurus Molecular Clouds

Figure 1.4 (left) shows the 13 CO total column density map of the Taurus molecular cloud (Mizuno
et al 1995) whose distance is 140 pc far from the Sun. Since 13 CO contains 13 C, a rare isotope ofC, the abundance of 13 CO is much smaller than that of 12 CO (for isotopic abundances, see Table
H.3). Owing to the low abundance, the emission lines of 13 CO are relatively optically thiner thanthat of 12 CO. Using 13 CO line, we can see deep inside of the molecular cloud. The distributions of
T Tauri stars and 13 CO column density coincide with each other. Since T Tauri stars are young
pre-main-sequence stars with M ∼ 1M , which are in the Kelvin-Helmholtz contraction stage anddo not reach the main-sequence Hydrogen burning stage, it is shown that stars are newly formed inmolecular clouds.
Since 18 O is much more rare isotope (18 O/16 O  13 C/12 C, see Table H.3), the distribution ofmuch higher-density gases is explored using C18 O lines. Figure 1.4 (right) shows C18 O map of Heiles
cloud 2 in the Taurus molecular cloud by Onishi et al (1996). This shows us that there are many
molecular cloud cores which have much higher density than the average. Many of these molecular
cloud cores are associated with IRAS infrared sources and T Tauri stars. It is shown that starCHAPTER 1. INTRODUCTION

4

Figure 1.4: (Left) 13 CO total column density map of the Taurus molecular cloud (Mizuno et al 1995).
Taken from their home page with URL of http://www.a.phys.nagoya-u.ac.jp/nanten/taurus.html (inJapanese). T Tauri stars, which are thought to be pre-main-sequence stars in the Kelvin-Helmholtz
contraction stage, are indicated by bright spots. (Right) C18 O map of Heiles cloud 2 region in the
Taurus molecular cloud (Onishi et al. 1996). This shows clearly that the cloud is composed of anumber of high-density regions.
formation occurs in the molecular cloud cores in the molecular cloud. They found 40 such cores in
the Taurus molecular cloud. Typical size of the core is ∼ 0.1 pc and the average density of the coreis as large as ∼ 104 cm−3 . The mass of the C18 O cores is estimated as ∼ 1 − 80M .
H13 CO+ ions are excited only after the density is much higher than the density at which CO
molecules are excited1 . H13 CO+ ions are used to explore the region with higher density than thatobserved by C18 O. Figure 1.5 shows the map of cores observed by H13 CO+ ions. The cores shown in
the lower panels are accompanied with infrared sources. The energy source of the stellar IR radiationis thought to be maintained by the accretion energy. Here we estimate the accretion luminosity Lacc(energy released by accretion in a unit time). Since the gravitational potential energy per unit massat the surface of a protostar with a radius R∗ and a mass M∗ is equal toΦ  −GM∗ /R∗ ,(1.1)

the kinetic energy of the gas accreting on the stellar surface is approximately equal to ∼ GM∗ /R∗ .
The energy inﬂow rate owing to the accretion is
∼

GM∗
c3
GM∗
× Ṁ ∼
×A s,
R∗
R∗
G

(1.2)

where Ṁ = Ac3s /G is the mass accretion rate (cs represents the sound speed of accreting gas). Thisrate was estimated as follows: (1) A mass of gravitationally bounded object (the Jeans mass) is given
as MJ ∼ ρ × c3s t3ﬀ , where tﬀ represents the gravitational free-fall time tﬀ ∼ (Gρ)−1/2 (eq.[2.26]).
(2) Thus, the accretion rate is estimated as Ṁ ∼ MJ /tﬀ ∼ ρ × c3s t2ﬀ ∼ c3s /G. Thus, the “natural”
accretion rate is proportional to the sound speed cs cubed and the inverse of gravitational constant
G. (See also §§ 4.5.1 and 4.5.2.). Accretion rate of−6Ṁ ∼ 1.6 × 10
1

M yr

−1



cs
190 m s−1

3

(1.3)

In section 2.10.5 of the next chapter, we will touch the critical density, above which the upper state level are well
excited by collision with H2 . In Table 2.1 we show the values of critical densities of CO and HCO+ .

1.3. T TAURI STARS

5

Figure 1.5: Pre-protostellar vs protostellar cores (H13 CO+ map). Upper panel shows the C18 O coreswithout associated IRAS sources. Lower panel shows the cores with IRAS sources. Taken from Fig.1
of Mizuno et al. (1994)
gives the accretion luminosity of

Lacc ∼ 5L

M∗
1M



10R
R∗



cs
190 m s−1

3

.

(1.4)

In the upper panels of Figure 1.5, the cores without IR sources are shown. This core does not show
accretion but collapse. That is, before a protostar is formed, the core itself contract owing to thegravity, which is explained closely in chapter 4.
In Figure 1.5, H13 CO+ total column density maps of the C18 O cores are shown. Cores in thelower panels have associated IRAS sources, while the cores in the upper panels have no IRAS sources.
Since the IRAS sources are thought to be protostars or objects in later stage, the core seems to evolvefrom that without an IRAS sources to that with an IRAS source. From this, the core with an IRAS
source is called protostellar core, which means that the cores contain protostars. On the other
hand, the core without IRAS source is called pre-protostellar core or, in short, pre-stellar core.
Figure 1.5 shows that the prestellar cores are less dense and more extended than the protostellar
core. This seems to suggest the density distribution around the density peak changes between before
and after the protostar formation.

1.3

T Tauri Stars

T Tauri stars are observationally late-type stars with strong emission lines and irregular light variations associated with dark or bright nebulosities. T Tauri stars are thought to be low-mass pre-mainsequence stars, which are younger than the main-sequence stars. Since these stars are connecting
between protostars and main-sequence stars, they attract attention today. More massive counterparts are called Herbig Ae-Be stars. They are doing the Kelvin-Helmholtz contraction in which theown gravitational energies released as it contracts gradually and this is the energy source of theluminosity.
The gravitational energy of a star with mass of M∗ and radius of R∗ is given as
W = −α

GM∗2
,
R∗

6

CHAPTER 1. INTRODUCTION

Figure 1.6: HR diagram of T Tauri stars. Many emission lines are found in the spectra of T Tauristars. WTTS (Weak Emission T Tauri Stars) and CTTS (Classical T Tauri Stars) are classiﬁed by
the equivalent widths of emission lines. That is, the objects with an EW of Hα emission < 1nm is
usually termed a WTTS. Taken from Fig.1.2 of Hartmann (1998).

1.4. SPECTRAL ENERGY DISTRIBUTION (SED)

7

where α is a numerical factor of order of unity which depends on the radial density distribution ρ(r).
From the Virial theorem [eq.(2.118)], total energy of E = W/2 = −U for monoatomic gas γ = 5/3.
When the star loses its energy by radiation in a rate of L = −dE/dt, the gravitational energy
decreases as dW/dt = −2L∗ (in other words, absolute value of the gravitational energy increases).
Since dW/dt = αGM∗2 /R∗2 · dR∗ /dt = −2L∗ , the star is contracting in a rate of
2L∗ R∗2
dR∗
=−
.
dt
αGM∗2
This process is called Kelvin-Helmholtz contraction. In short, the gravitational energy is used tocompensate the luminosity. When the system loses its energy, although the potential energy W
decreases, its thermal energy U increases. Thus, the central temperature increases and such a starreaches to the main-sequence star in which the nuclear fusion reaction of 1 H to create 4 He produces
energy.
ProblemObtain the Kelvin-Helmholtz time of our Sun. The Kelvin-Helmholtz time is given
τKH ∼

GM∗2
.
R∗ L∗

Many emission lines are observed in the spectra of T Tauri stars. WTTS (Weak emission T Tauri
Stars or Weak line T Tauri Stars) and CTTS (Classical T Tauri Stars) are classiﬁed by their equivalent
widths (EWs) of emission lines. That is, the objects with an EW of Hα emission < 1nm = 10Å is
usually termed a WTTS. Figure 1.6 is the HR diagram (Teﬀ − Lbol ) of T Tauri stars in Taurus-Auriga
region (Kenyon & Hartman 1995). WTTSs distribute near the main-sequence while CTTSs are foundeven far from the main-sequence. A number of theoretical evolution tracks for pre-main-sequence stars
with M ∼ 0.1M − 2.5M are shown in a solid line, while the isochorones for ages of 105 yr, 106 yr,
and 107 yr are plotted in a dashed line. Vertical evolutionary paths are the Hayashi convective track.
Since D=2 H has a much lower critical temperature (and density) for a fusion nuclear reaction to makeHe than 1 H, Deuterium begins to burn before reaching the zero-age-main sequence. This occurs near
the isochrone for the age of 105 yr and some activities related to the ignition of Deuterium seem tomake the central star visible (Stahler 1983).
Disk Frequency
Infrared studies of T Tauri stars in star-forming regions have suggested that initial disk frequency
is rather high and that the disk lifetimes are relatively short 3 − 15Myr. From JHKL photometry,
Haish, Lada, & Lada (2001) obtained the fraction of disk-bearing stars for 6 star formation regions.
L band excess emission indicates an accompanied disk. The fraction is a decreasing function of ageof the cluster as in Figure 1.7. This ﬁgure shows clearly that the disk fraction is initially very high(>∼80%) and rapidly decreases with increasing cluster age. In 3 Myr a half of the disk stars lose theirdisks. Overall disk lifetime is estimated as ∼ 6Myr.

1.4

Spectral Energy Distribution (SED)

A tool to know the process of star formation is provided by the spectral energy distribution (SED)
mainly in the near- and mid-infrared light. T Tauri stars and protostars have typical respective SEDs.

8

CHAPTER 1. INTRODUCTION

Figure 1.7: Fractions of IR excess sources in respective clusters are plotted against the age of theclusters.
IR SEDs of T Tauri stars were classiﬁed into three as Class I, Class II, and Class III, from a stand-pointof relative importance of the radiation from a dust disk to the stellar black-body radiation. Today,
the classiﬁcation is extended to the protostars, which is precedence of the T Tauri stars, and they
are called Class 0 objects. (Unfortunately, there is no zero in Roman numerals.) The classiﬁcation
uses the slope in near- and mid-IR spectrum as α ≡ |d log λFλ /d log λ|λ=2.2μm−25μm . In Figure 1.8,
typical SEDs and models of emission regions are shown.
1. Class III is well ﬁtted by a black-body spectrum, which shows the energy mainly comes froma central star. This SED is observed in the weak-line T Tauri stars. Although T Tauri stars
show emission lines of such as the Hydrogen Balmar sequence, the weak-line T Tauri stars do
not show prominent emission lines, which indicates the amount of gas just outside the star
(this seems to be supplied by the accretion process) is small. In this stage, a disk has been
disappeared or only an extremely less-massive disk is still alive. α < −2.
2. Class II SED is ﬁtted by a single-temperature black-body plus excess IR emission. This showsthat there is a dust disk around a pre-main sequence star and it is heated by the radiationfrom a central star. The width of the spectrum of the disk component is much wider than that
expected from a single-temperature black-body radiation. Thus, the disk has a temperature
gradient which decreases with increasing the distance from the central star. In this stage, thedust disk is more massive than that of Class III sources. Classical T Tauri stars have such
SEDs. −2 < α < 0.
3. In Class I SED, the mid infrared radiation which seems to come from the dust envelope ispredominant over the stellar black-body radiation. Since the stellar black-body radiation seems
to escape at least partially from the dust envelope, a relatively large solid angle is expected fora region where the dust envelope does not intervene. α > 0.
4. Class 0 SED seems to be emitted by an isothermal dust with ∼ 30K. The protostar seems to
be completely covered by gas and dust and is obscured with a large optical depth by the dust

1.4. SPECTRAL ENERGY DISTRIBUTION (SED)

9

Figure 1.8: Spectral Energy Distribution (SED) of young stellar objects (YSOs) and their models.
(Left:) ν − νFν plot taken from Lada (1999). (Right:) λ − λFλ plot taken from André (1994)
envelope. No contribution can be reached from the stellar-black body radiation.
The reason why the emission from the disk becomes wide in the spectral range is understood(Fig.1.9) as follows: Temperature of the disk is determined by a balance of heating and cooling.
Assuming the disk is geometrically thin but optically thick, the cooling per unit area is given by
the equation of the black-body Planck radiation. Therefore, the temperature is determined by theheating predominantly by viscous heating and extra heating by the radiation from the central star.
The ﬂux density emitted by the disk is given by
νFν ∼



νπBν [T (r)]2πrdr ∼ r(T or ν)2 νBν .

(1.5)

Assuming the radial distribution of temperature as

T = T0

R
R0

−q

,

(1.6)

CHAPTER 1. INTRODUCTION

10

Figure 1.9: Explanation for the spectral index of the emission from a geometrically thin but opticallythick disk. Taken from Fig.16 of Lada (1999).
(q = 3/4 for the standard accretion disk) and taking notice that each temperature in the disk radiates
at a characteristic frequency ν ∝ T (the Wien’s law for black-body radiation)
νFν ∼ r 2 νBν ∝ ν 4 T −2/q ∝ ν 4−2/q ,(1.7)

where we used the fact that the peak value of Bν ∝ ν 3 . Therefore, it is shown thatνFν ∝ ν n ;2n =4− .
q

(1.8)

As shown in the previous section, we have no young stellar objects found by IR before a protostar
is formed. These kinds of objects (pre-protostellar core) are sometimes called Class −1. The classiﬁcation was originally based on the SED and did not exactly mean an evolution sequence. However,today YSOs are considered to evolve as the sequence of the classes: Class −1 → Class 0 → Class I→ Class II → Class III → main-sequence star.

1.5

Protostars

1.5.1

B335

Infall Motion
B335 is a dark cloud (Fig.1.10) with a distance of D  250pc. Inside the dark cloud, a Class 0 IR
source is found. The object is famous for the discovery of gas infall motion. In Figure 1.11, theline proﬁles of CS J = 2 − 1 line emissions are shown (Zhou et al. 1993). The relative position of
the proﬁles correspond to the position of the beam. (9,9) represents the oﬀset of (9”,9”) from thecenter. At the center (0,0), the spectrum shows two peaks and the blue-shifted peak is brighter than
the red-shifted one. This is believed to be a sign of gas infall motion. The blue-red asymmetry is
explained as follows:

1.5. PROTOSTARS

11

Figure 1.10: Near infrared images of B335, which is Class 0 source. Left is H band image and rightis K band image.

Figure 1.11: Line proﬁle of CS J = 2 − 1 line radio emission. Model spectra illustrated in a dashedline (Zhou 1995) are overlaid on to the observed spectra in a solid line (Zhoug et al 1993).

12

CHAPTER 1. INTRODUCTION

Figure 1.12: Explanation of blue-red asymmetry when we observe a spherical symmetric inﬂow motion. An isovelocity curve for the red-shifted gas is plotted in a solid line. That for the blue-shifted
gas is plotted in a dashed line. Taken from Fig.14 of Lada (1999).
1. Considering a spherical symmetric inﬂow of gas, whose inﬂow velocity vr increases with reaching
the center (a decreasing function of r)
2. Considering a gas element at r moving at a speed of vr (r) < 0, the velocity projected on aline-of-sight is equal to(1.9)
vline−of−sight = vsystemic + vr cos θ,
where vsystemic represents the systemic velocity of the cloud (line-of-sight velocity of the cloudcenter) and θ is the angle between the line-of-sight and the position vector of the gas element.
The isovector lines, the line which connect the positions whose procession/recession velocities
are the same, become like an ellipse shown in Fig.1.12.
3. An isovelocity curve for the red-shifted gas is plotted in a solid line. That for the blue-shifted
gas is plotted in a dashed line. If the gas is optically thin, the blue-shifted and red-shifted gases
contribute equally to the observed spectrum and the blue- and red-shifted peaks of the emission
line should be the same.
4. In the case that the gas has a ﬁnite optical depth, for the red-shifted emission line a cold gas in
the fore side absorbs eﬀectively the emission coming from the hot interior. On the other hand,for the blue-shifted emission line, the emission made by the hot interior gas escapes from thecloud without absorbed by the cold gas (there is no cold blue-shifted gas).
5. As a result, the blue peak of the emission line becomes more prominent than that of the redshifted emission. This is the explanation of the blue-red asymmetry.
In Figure 1.11, model spectra calculated with the Sobolev approximation or LVG (large-velocity
gradient; Zhou 1995) are shown. These show the blue-red asymmetry (the blue line > the red line).

1.5. PROTOSTARS

13

Figure 1.13: C18 O total column density map (left) and H13 CO+ channel map (right) of B335 alongwith the position-velocity maps along the major and minor axes. Taken from Fig.3 of Saito et al(1999).

Figure 1.14: Column density distribution NH (r) derived from the H13 CO+ and C18 O data taken by
the Nobeyama 45 m telescope. Taken from Fig.9 of Saito et al (1999).

CHAPTER 1. INTRODUCTION

14
Outﬂow Motion and a Disk

Early observation of star forming regions have revealed fast molecular outﬂow are often ejected−1from the central protostar with >∼10km s . B335 is also a typical outﬂow source. In Figure 1.13,
distributions of high density gases traced by the C18 O and H13 CO+ lines are shown as well as the
bipolar outﬂow whose outline is indicated by a shadow (Hirano et al 1988). Comparing left and rightpanels, it is shown that the distribution of C18 O gas is more extended than that of H13 CO+ which
traces higher-density gas. And the distribution of the H13 CO+ is more compact and the projectedsurface density seems to show the the actual distribution is spherical. And the molecular outﬂow
seems to be ejected in the direction of the minor axis of the high-density gas. It may suggest that (1) amolecular outﬂow is focused or collimated by the eﬀect of density distribution or that (2) collimation
is made by the magnetic ﬁelds which run preferentially perpendicularly to the gas disk. This gas disk
is observed by these high density tracers.
Combining the C18 O and H13 CO+ distributions, the surface density distribution along the major
axis is obtained by Saito et al (1999). From the lower panel of Figure 1.14, the column density
distribution is well ﬁtted in the range from 7,000 to 42,000 AU in radius,Σ(r) = 6.3 × 1021 cm−2



r
4
10 AU

−0.95,(1.10)

where they omitted the data of r <∼7000 since the beam size is not negligible. Similar power-law
density distributions are found by the far IR thermal dust emission.

1.5.2

L1551 IRS 5

Figure 1.15: (Left:) 13 CO column density distribution. The contour lines represent the distributionof 13 CO column density. 2.2 μm infra-red reﬂection nebula is shown in grey scale which was observedby Hoddap (1994). (Right:) Schematic view of L1551 IRS5 region.
L1551 IRS 5 is one of the most well studied protostellar objects. This has an infra-red emission
nebulosity (Fig.1.15). It is believed that there is a hole perpendicular to the high-density disk and theemission from the central star escapes through the hole and irradiate the nebulosity. In this sense this1.5. PROTOSTARS

15

Figure 1.16: Isovelocity contours measured by the 13 CO J = 1 − 0 line. It should be noticed that theisovelocity lines run parallelly to the major axis. The north-eastern side shows a red-shift and thesouth-western side shows a blue-shift.
is a reﬂection nebula. L1551 IRS 5 has an elongated structure of dense gas similar to that observedin B335. The gas is extending in the direction from north-west to south-east [Fig.1.15 (left)]. Since
the opposite side of the nebulosity is not observed, the opposite side of nebulosity seems to be locatedbeyond the high-density disk and be obscured by the disk. This is possible if we see the south surfaceof the high-density disk as in Figure 1.15 (right).
Infall Motion
The inﬂow motion is measured. Figure 1.16 shows the isovelocity contours measured by the 13 CO
J = 1−0 observation (Ohashi et al 1996). It should be noticed that the isovelocity lines run parallelly
to the major axis. The north-eastern side shows a red-shift and the south-western side shows a blueshift. Considering the conﬁguration of the gas disk shown in Fig.1.15 (right), this pattern of isovelocity
contours indicates not outﬂow but inﬂow. That is, the north-east side is a near side of the disk and thesouth-west side is a far side. Since a red-shifted motion is observed in the near side and a blue-shiftedmotion is observed in the far side, it should be concluded that the gas disk of the L1551 IRS5 is now
infalling.
Optical Jet
HST found two optical jets emanating from L1551 IRS5. This has been observed by SUBARU
telescope, which found jet emission is dominated by [FeII] lines in the J- and H-bands. The jet
extents to the south-western direction and disappears at ∼ 10”  1400AU from the IRS5. The
width-to-length ratio is very small <∼1/10 or less, while the bipolar molecular outﬂow shows a lesscollimated ﬂow. As for the origin of the two jets, these two jets might be ejected from a single source.
However, since there are at least two radio continuum sources in IRS5 within the mutual separationof ∼ 0.”5 [see Fig.1.17 (right)], these jets seem to be ejected from the two sources independently.

16

CHAPTER 1. INTRODUCTION

Figure 1.17: (Left:) Infrared image (J- and K-band) of the IR reﬂection nebula around L1551 IRS5
by SUBARU telescope. Taken from Fig.1 of Itoh et al. (2000). (A jpeg ﬁle is available from thefollowing url: http://SubaruTelescope.org/Science/press release/9908/L1551.jpg). (Right:)
Central 100 AU region map of L1551 IRS5. This is taken by the λ = 2.7cm radio continuum
observation. Deconvolved map (lower-left) shows clearly that IRS5 consists of two sources. Taken
from Looney et al. (1997).

1.6. L 1544: PRE-PROTOSTELLAR CORES

17

Figure 1.18: A mosaic image of HH 111 based on HST NICMOS images (left) and WFPC2 images(right). Taken from Fig.1 of Reipurth et al (1999) and rotated clockwise.
Although the lengths of these jets are restricted to 10”, Herbig-Haro jets, which are much larger
than the jets in L1551 IRS5, have been found. HH30 has a ∼ 500 AU-scale jet whose emission ismainly from the shock-excited emission lines. One of the largest ones is HH111, which is a member of
the Orion star forming region and whose distance is as large as D ∼ 400pc, and a jet with a length of∼ 4pc is observed. Source of HH111 system is thought to consist of at least binary stars or possibly
triple stars [Reipurth et al (1999)]. Star A, which coincides with a λ = 3.6 cm radio continuum
source (VLA 1), shows an elongation in the VLA map whose direction is parallel to the axis of thejet. Therefore, star A is considered to be a source of the jet. Since the VLA map of star A shows
another elongated structure perpendicular to the jet axis, star A may be a binary composed by two
outﬂow sources.

1.6

L 1544: Pre-protostellar Cores

L1544 is known as a pre-protostellar core (Taﬀala et al 1998). That shows an infall motion but thiscontains no IR protostars. In Figure 1.19(left), CCS total column density map is shown, which shows
an elongated structure. Ohashi et al (1999) have found both rotation and infall motion in the cloud.
Position-velocity (PV) diagram along the minor axis shows the infall motion. That along the major
axis indicates a rotational motion, which is shown by a velocity gradient. Due to a ﬁnite size of thebeam, a contraction motion is also shown in the PV diagram along the major axis.

1.7

Rotation of Outﬂow

Figures 1.20 and 1.21 show that the outﬂow found in globule CB 26 is rotating. CB 26 is a Bok globulewith bipolar NIR refrection nebulae between which a very young T Tauri star and ∼ 200AU-scale
high-density disk are observed (Fig.1.20). Figure 1.21 shows the map of intensity-weighted velocity(1st moment)
T (v)vdv,(1.11)
Vave = T (v)dv
where T (v) represents the brightness temperature for line-of-sight velocity v. Left panel indicatesclearly that the right part gas is departing from us while the left part is approaching. This means

CHAPTER 1. INTRODUCTION

18

Figure 1.19: CCS image of prestellar core L1544. (Left:) Total intensity map. (Right:) PV diagrams
along the minor axis (left) and along the major axis (right).
there exists a global velocity gradient perpendicular to the ﬂow axis, or in other words, the rotation
motion in which the rotation axis coincides with the symmetric axis. The rotation is toward the same
direction of the high-density disk observed by 13 CO (Launhardt & Sargent 2001). Right panel is an
expecting intensity-weighted velocity distribution for a simple model 12 CO (J = 2 − 1),Vr (r) = v0



r
,
r0


(1.12)
−1

Vrot = VKep


,

r
,
r0
 −q
r
,
TK (r) = T0
r0n(r) = n0(1.13)0−3(1.14)(1.15)

in which we assume (i) the outﬂow is conical, (ii) gas element at ( , z) in cylindrical coordinate waslaunched from a Keplerlian disk at 0 conserving angular momentum, (iii) the radial expanding speed
Vr (r) is simply proportional to the distance from the central star r (eq.[1.12]). The assumption (ii)leads to (eq.[1.13) as the rotation speed is inversely proportional to the distance from the rotation axis
= (r 2 − z 2 )1/2 . (iv) Density and temperature distributions are assumed as the density decreases
with the distance r in proportion to r −3 (eq.[1.14]), and the kinetic temperature decreases with the
distance r in proportion to r −q (eq.[1.14]).
Comparing with a simple model (right panel), such global rotation motion is seen evidently only
when the outﬂow is observed from edge-on.

1.7. ROTATION OF OUTFLOW

19

Figure 1.20: Globule CB 26 has an outﬂow which seems to be seen edge-on. Grey-scale indicates the
K-band image of the bipolar refrection nebula. Red contours which show SMA 1.1 mm dust continuum
emission indicate a high-density disk exists between two lobes of the bipolar nebula. Green contours
represent the 12 CO(J = 2 − 1) integrated intensity.

Figure 1.21: Intensity weighted velocity is shown (left). Right one is a model.

CHAPTER 1. INTRODUCTION

20

Figure 1.22: Explanation how the polarized radiation forms. Taken from Weintraub et al.(2000).

1.8

Magnetic Fields

Directions of magnetic ﬁeld are studied by (1) measuring the polarization of light which is suﬀeredfrom interstellar absorption. In this case the direction of magnetic ﬁeld is parallel to the polarization
vector. The reason is explained in Figure 1.22. In the magnetic ﬁelds, the dusts are aligned in a
way that their major axes are perpendicular to the magnetic ﬁeld lines. Such aligned dusts absorb
selectively the radiation whose E-vector is parallel to their major axes. As a result, the detected light
has a linear polarization parallel to the magnetic ﬁeld lines.
However, the polarization measurement in the near IR wavelength limited to the region withlow gas density, because background stars suﬀer severe absorption and becomes hard to be observed
if we want to measure the polarization of the high-density region. More direct method is (2) the
measurement of the linear polarization of the thermal emission from dusts in the mm (or sub-mm)
wavelengths; in this case the direction of magnetic ﬁeld is perpendicular to the polarization vector.
The mechanism is explained in Figure 1.22b. The aligned dusts, whose major axes are perpendicular
to the magnetic ﬁeld lines, emit the radiation whose E-vector is parallel to the major axes. Since theabsorption does not have a severe eﬀect in this mm wavelengths, this gives information about themagnetic ﬁelds deep inside the clouds.

1.8. MAGNETIC FIELDS

21

Figure 1.23: H band (left) and Ks band (right) polarizations shown by bars and λ = 1.1 mm dust
continuum image by ASTE/AzTEX. Distance is 260pc and this image covers ∼ 1.5 pc × 1.5 pc.

1.8.1

Serpens-South Cloud

Figure 1.23 shows the structure of the magnetic ﬁeld measured with H and Ks band polarimetory(Sugitani et al 2010). Serpens south cluster is located at the center of an elongated ﬁlament. At thesame time the location corresponds to the intersection of several minor ﬁlaments. This ﬁgure shows
clearly that magnetic ﬁeld runs in a perpendicular direction to the axis of the ﬁlamentary cloud.
Although the conﬁguration of magnetic ﬁeld lines in the ﬁlament is not explored suﬃciently by thenear IR polarimetory, in the envelope the magnetic ﬁeld direction is perpendicular to the ﬁlament aswell as outside of the ﬁlament.

1.8.2

Prestellar Core

Figure 1.24 illustrates the polarization maps of three prestellar cores (L1544, L183, and L43) donein the 850 μm band by JCMT-SCUBA. In L1544 and L183 the mean magnetic ﬁelds are at an angleof 30 deg to the minor axes of the cores. L43 is not a simple object; there is a T Tauri star locatedin the second core which extends to south-western side of the core (an edge of this core is seen near
the western SCUBA frame boundary). And a molecular outﬂow from the source seems to aﬀect
the core. The magnetic ﬁeld as well as the gas are swept by the molecular outﬂow. L43 seems an
exception. The fact that the mean magnetic ﬁelds are parallel to the minor axis of the high-density
gas distribution seems to mean gas contracts preferentially in the direction parallel to the magnetic
ﬁelds.
Figure 1.25(left) shows the optical image of dark cloud B68. Emission from background stars
suﬀers extinction due to the dust in the dark cloud. Figure 1.32(right) gives the distribution ofvisual extinction AV against the distance from the center. The distribution is well ﬁtted to that of
the isothermal cloud in a hydrostatic equilibrium. Thus, the amount of extinction gives an estimateof the column density of interstellar dusts in the dark cloud. Since the absorption by the aligned
dust grains induces the polarization, the direction of magnetic ﬁeld is obtained from the E vector of
the polarized emission. Figure 1.25(right) indicates the polarization vector and thus the direction of

22

CHAPTER 1. INTRODUCTION

Figure 1.24: Directions of B-Field are shown from the linear polarization observation of 850 μm thermal emission from dusts by JCMT-SCUBA. L1544 (upper-left) and L183 (upper-right), the magnetic
ﬁeld and the minor axis of the molecular gas distribution coincide with each other within ∼ 30deg.
Conﬁguration of magnetic ﬁeld in L43 (lower) is not simple, which may be aﬀected by the molecular
outﬂow. Taken from Ward-Thompson et al (2000).

1.8. MAGNETIC FIELDS

23

Figure 1.25: (left) Optical image of B68 (Alves et al. 2001) and (right) H band image and polarization
E-vector derived for background stars (Kandori et al. 2009). The region shown in (b) correspends to
the region shown by the square in (a).
interstellar magnetic ﬁeld in B68. Figure indicates a magnetic ﬁeld line like a “hour-glass.”

1.8.3

Cores with Protostars

Bok Globules
Bok globules are isolated dark clouds. Wolf, Launhardt & Henning (2003) have studied the relation
between the magnetic ﬁeld and outﬂow directions. Figure 1.26 displays the direction of the magnetic
ﬁeld obtained with JCMT-SCUBA polarization observation in 850μm as well as the outﬂow foundby 13 CO observation of Chandler & Sargent (1993) and 12 CO observation by Yun & Clemens (1994).
In B 335, CB 230 and CB 244 outﬂows are oriented in the direction perpendicular to the major
axis of the globules. The magnetic ﬁeld is running parallel to the outﬂow but perpendicular to themajor axis. This means that a disk is formed with a gas ﬂowing along the magnetic ﬁeld. Further, theoutﬂow is generated by a twisted magnetic ﬁeld which is an outcome of a rotating disk. In globule CB
54, the alignment is not perfect, that is, magnetic ﬁeld direction is slightly aligned with the outﬂow(Wolf, Launhardt & Henning 2003). This seems to strengthen the magnetic origin of the outﬂow andmagnetically guided disk formation.
It is clearly shown that the polarization anti-correlates with the intensity of the thermal emission.
This might be due to the observational error in the low intensity region (or low S/N region). Thismay be related to physical processes such as (1) the alignment owing to the magnetic ﬁeld becomesineﬃcient in high density region or (2) the magnetic ﬁeld is tangled in the dense region and thepolarization due to the aligned dust is canceled.

1.8.4

Magnetic Field ∼1000AU Scale

Single dish observation with JCMT can reveal a structure larger than 5000-10000AU (Fig.1.26).
Smaller-scale conﬁguration of magnetic ﬁeld is only able to be reached with interferometric observationsuch as SMA. Figure 1.27 indicates a 1000 AU-scale conﬁguration of magnetic ﬁeld for a low-massstar forming region NGC 1333 IRAS 4A (Girart, Rao & Marrone 2006). Although this is a complex
region, a core extends from the upper-left to the lower-right. The magnetic ﬁeld (Fig.1.27(right))
runs from the upper-right to the lower-left and this clearly shows ’hour-glass’ shape. This means thatthe magnetic ﬁeld is squeezed by the eﬀect of contraction in the star formation process.

24

CHAPTER 1. INTRODUCTION

Figure 1.26: Directions of magnetic ﬁeld are shown from the linear polarization observation of 850
μm thermal emission from dusts by JCMT-SCUBA. These four objects are known as Bok globules.
The bars with γ represent the average directions of polarization (E-vector), which are perpendicular
to the magnetic ﬁeld.

Figure 1.27: SMA observation low-mass (NGC1333 IRAS 4A) star formation region. (middle) The
total intensity of 879 μm thermal dust emission (contour lines), that of polarized intensity (false color)and polarization E-vector (red bars) are shown. (right) The polarization E-vector of the dust thermal
emissions is perpendicular to the direction of magnetic ﬁeld (Fig.1.22). The direction of magnetic
ﬁeld (red bars) is overlaid on the total intensity distribution (false color). 1 = 300AU. This clearlyshows that the magnetic ﬁeld looks like a hour-glass within ∼ 1000 AU scale.

1.8. MAGNETIC FIELDS

25

Figure 1.28: SMA observation high-mass (G31.41+0.31) star formation region. (left) The total
thermal dust emission at 879 μm is shown with contour lines, while the intensity of polarized intensity
is shown with false-color. 1 = 7000AU. Directions of magnetic ﬁeld are illustrated with solid thickbars. (middle) Directions of magnetic ﬁeld is overlaid on the ﬂux-weighted velocity map of the CH3 OH
147 -156 line, which indicates the core has a systematic velocity gradient.
They ﬁtted the magnetic ﬁeld line by a polynomial as
y = g(1 + Cx2 ),(1.16)

where x and y, distances along the symmetic axis and along the major axis of the disk, respectively.
They obtained C = 0.12 ± 0.06sec−1 .
SMA observation for G31.41+0.31 (Fig. 1.28) also reveals that the magnetic ﬁeld conﬁguration of
10000AU scale for high-mass star forming region is “hour-glass” shape (Girart et al. 2009). This may
indicate the magnetic ﬁeld is draged by the infall motion. However, this can not be ﬁtted with the
conﬁguration of the magnetic ﬁeld line obtained with single dish observations, since Fourier component
with the small wave-number k is not suﬃciently picked up in the interferometric observation. In otherwords, large-scale (more or less uniform) magnetic ﬁeld is underestimated in the interferometricobservations.
T Tauri Disks
Magnetic ﬁeld at the position of protostars and T Tauri stars are measured for IRAS 16293-2422,
L1551 IRS5, NGC1333 IRAS 4A, and HL Tau (Tamura et al. 1995). Although HL Tau is a T Tauri
star, it has a gas disk. Thus this is a Class I source. The others are believed to be in protostellar
phase (Class 0 sources). It is known that IRAS 16293-2422, L1551 IRS5, and HL Tau have diskswith the radii of 1500-4000 AU from radio observations of molecular lines. Further, near infraredobservations have shown that these objects have 300-1000 AU dust disks. Figure 1.29 shows the
E-vector of polarized light. If this is the dust thermal radiation, the direction of magnetic ﬁelds isperpendicular to the polarization E-vector. Figure shows the magnetic ﬁelds run almost perpendicular

26

CHAPTER 1. INTRODUCTION

Figure 1.29: Polarization of the radio continuum λ = 1mm, λ = 0.8mm. IRAS 16293-2422 (upper-left:
Class 0), L1551 IRS5 (upper-right: Class 0), NGC1333 IRAS 4A (lower-left: Class 0), and HL Tau(lower-right: Class I). Taken from Tamura et al (1995).
to the elongation of the gas disk. Global directions of magnetic ﬁeld outside the gas disk and the
direction of CO outﬂows are also shown in the ﬁgure by arrows. It is noteworthy that the directionsof local magnetic ﬁelds, global magnetic ﬁelds, and outﬂows coincide with each other within ∼ 30
deg.

1.9. DENSITY DISTRIBUTION

1.9

27

Density Distribution

Motte & André (2001) made λ =1.3 mm continuum mapping survey of the embedded young stellar
objects (YSOs) in the Taurus molecular cloud. Their maps include several isolated Bok globules,
as well as protostellar objects in the Perseus cluster. For the protostellar envelopes mapped in
Taurus, the results are roughly consistent with the predictions of the self-similar inside-out collapse
model of Shu and collaborators (§4.5.1). The envelopes observed in Bok globules are also qualitatively
consistent with these predictions, providing the eﬀects of magnetic pressure are included in the model.
By contrast, the envelopes of Class 0 protostars in Perseus have ﬁnite radii <∼10000 AU and are afactor of 3 to 12 denser than is predicted by the standard model.
Another method to measure the density distribution is to use the near IR extinction. From
(H − K) colors of background stars, the local value of AV in a dark cloud can be obtained using astandard reddening law,
AV = 15.87E(H − K)(1.17)
if the intrinsic colors of background stars are known. Here, the color excess is deﬁned as the diﬀerence
between the observed color and the intrinsic color: E(H − K) ≡ (H − K)obs − (H − K)intrinsic . We
can convert the extinction to the column density assuming the gas/dust ratio is constant
N (H + H2 ) = 2 × 1021 cm−2 mag−1 AV .

(1.18)

This is a standard method to obtain the local column density of the dark cloud using the near IR
photometry.
See Figure 1.30. If the density distribution is expressed asρ(r) = ρ0

r
r0

−α

,

(1.19)

where r is a physical distance from the center. The column density distribution against the projecteddistance of the line-of-sight from the center of the cloud is given as (R2 −p2 )1/2

Nρ (p) = 2ρ (s2 + p2 )1/2 ds,(1.20)

0where R represents the outer radius of the cloud. Using equation (1.18), this yields AV distribution−23

AV (p) = 10

ρ0 r0α (R2 −p2 )1/2

(s2 + p2 )−α/2 ds.

(1.21)

0If background stars are uniformly distributed, the number of stars with AV |obs is proportional to thearea which satisﬁes AV |obs = AV (p). That is, if we plot AV (p) against 2πpdp, this gives the number
distribution of background stars with AV . Figure 1.32 shows the result of L977 dark cloud by Alveset al (1998).
Recently, Alves et al (2001) derived directly the radial distribution of NH by comparing the NH (p)model distribution for B68 (see also Kandori et al 2009). They obtained a distribution is well ﬁtted
by the Bonner-Ebert sphere in which a hydrostatic balance between the self-gravity and the pressure
force is achieved (lower panels of Fig.1.32) (see section 4.1.1).
In this ﬁeld, we should pay attention to the density distribution in cylindrical clouds. As seenin the Taurus molecular cloud, there are many ﬁlamentary structures in a molecular cloud. In §4.1,
we will give the distribution for a hydrostatic spherically symmetric and that of a cylindrical cloud.
The former is proportional to ρ ∝ r −2 and the latter is ρ ∝ r −4 . Therefore, the distribution ρ ∝ r −4

28

CHAPTER 1. INTRODUCTION

Figure 1.30: Schematic view to explain an AV distribution.

Figure 1.31: (Left:) Radial intensity proﬁles of the environment of 7 embedded YSOs (a-g) and 1starless core (h). (Right:) Same as left panel but for 4 isolated globules (a-d) and 4 Perseus protostars(e-h). Taken from Motte & Andre (2001).

1.9. DENSITY DISTRIBUTION

29

Figure 1.32: Density distribution of L977 (top) and B68 (bottom) dark clouds. (Top-left:) L977
dark cloud dust extinction map derived from the infrared (H-K) observations. (Top-right:) Observed
frequency distribution of extinction measurements for L977 and the predictions from clouds models
with density structures ρ(r) ∝ r −α having α = 1 (dashed line), 2 (solid line), 3 (dotted line), and4 (dash-dotted line). (Bottom-left:) B68 images (false color images made from B, V, and I images

CHAPTER 1. INTRODUCTION

30

Figure 1.33: A structure of magnetic ﬁelds in the L1641 region. Polarization of light from embedded
stars (Vrba et al. 1988) is shown by a bar. The direction of magnetic ﬁelds in the line-of-sight isobserved using the HI Zeeman splitting, which is shown by a circle and cross (Heiles 1989).
was expected for cylindrical cloud. From near IR extinctions observation (Alves et al 1998), even if
a cloud is rather elongated [Fig.1.32 (top-left)], the power of the density distribution is equal to not−4 but  −2. Fiege, & Pudritz (2001) proposed an idea that a toroidal component of the magnetic
ﬁeld, Bφ , plays an important role in the hydrostatic balance of the cylindrical cloud (Fig.1.33).

1.10

Mass Spectrum

We have seen that a molecular cloud consists in many molecular cloud cores. For many years, there
are attempts to determine the mass spectrum of the cores.
From a radio molecular line survey, a mass of each cloud core is determined. Plotting a histogramof number of cores against the mass, we have found that a mass spectrum can be ﬁtted by a powerlaw asdN= Mn(1.22)
dM
where dN/dM represents the number of cores per unit mass interval. Many observation indicate thatn ∼ −1.5 as Table 1.1.
Figure 1.34(left) (André, Ward-Thomson, & Barsony 2000) shows a mass spectrum function
dN/dM for 59 ρ Oph cloud prestellar cores obtained at the IRAM 30-m telescope with the MPIfR
19-channel bolometer array. Presetellar cores with a mass ∼ M has a spectrum of dN/dM ∝ M −2.5
(n = −2.5) for M >∼0.5M . While the right panel (Motte et al 2001) shows the cumulative mass
spectrum (N (> M ) vs. M ) of the 70 starless condensations identiﬁed in NGC 2068/2071. The mass
spectrum for the 30 condensations of the NGC 2068 sub-region is very similar in shape. The best-ﬁt
power-law is N (> M ) ∝ M n+1 ∝ M −1.1 above M >∼0.8M . That is, n = −2.1. This power derived

1.10. MASS SPECTRUM

31

Table 1.1: Mass spectrum indicies derived with molecular line surveys.
Papern
RegionObservation
Mass range
Loren (1989)
−1.1
ρ Oph
10M <∼M <∼300M
18
3<
Stutzki& −1.7 ± 0.15
M17 SW
C O (J = 2 − a few M ∼M <∼a few 10 M
34
Guesten (1990)1), S S (J =2−1, J = 3−2)
Lada et al. (1991) −1.6
L1630
CS (J = 2 − 1) M >∼20M
13
Nozawa et al. −1.7
ρ Oph North
CO (J = 1 − 3M <∼M <∼160M
(1991)
0)
Tatematsu et al. −1.6 ± 0.3
Orion A
CS (J = 1 − 0) M >∼50M
(1993)
13 CO (J = 1 − M >100M
Dobashi et al. −1.6
Cygnus∼(1996)
0)
Onish et al (1996) −0.9 ± 0.2
TaurusC18 O (J = 1 − 3M <
∼M <
∼80M
0)
Motte, André, & −1.5ρOph
λ = 1.3mm
M<∼0.5M
Neri (1998)−2.5
M>∼0.5M
Testi & Sargent −2.1
Sarpens
λ = 3mm
1M <
∼M <
∼30M
(1998)
12 CO,
13 CO, 10−4 M <M <104 M
Krameret −1.6 ∼ −1.8 L1457 etc∗∼ ∼
al.(1998)
C18 O (J = 1 −0, 2 − 1)
Heithausen et al. −1.84
MCLD 123.5 12 CO(J = 1−0 MJ <
∼M <
∼10M
(1998)+ 24.9, Polaris and 2 − 1)
FlareJohnstone et al. −2 ∼ −2.5
ρOph
λ = 850μm
M>∼0.6M
(2000)−1.5
M<∼0.6M
>Johnstone et al. −2.5 ∼ −3Ori B
λ = 850μm
M ∼M
(2001)
Onishietal. −2.5
Taurus
H13 CO+ (J = 3.5M <∼M <∼20.1M
(2002)
1 − 0)
Ikeda, Sunada, & −2.3 ± 0.1
Ori A
H13 CO+ (J = M >∼9M
Kitamura (2007)
1 − 0)−0.4 ± 0.1
M<∼9M
Ikeda & Kitamura −2.3 ± 0.3
Ori A
C18 O (J = 1 − M >∼5M
(2009)0)∗ MCLD126.6+24.5, NGC 1499 SW, Orion B South, S140, M17 SW, and NGC 7538

CHAPTER 1. INTRODUCTION

32

Figure 1.34: Cumulative mass distribution of the 70 pre-stellar condensations of NGC 2068/2071. The
dotted and dashed lines are power-laws corresponding to the mass spectrum of CO clumps (Kramer
et al. 1996) and to the IMF of Salpeter (1955), respectively. Taken from Fig.3 of Motte et al (2001).

from the dust thermal emission is diﬀerent from that derived by the radio molecular emission lines.
Table 1.1 summarizes the observations to calculate the power index of core mass function.
The mass function of newborn stars is called as initial mass function (IMF). IMF for ﬁeld stars inthe solar neighborhood has been obtained as shown in Figure 1.36. The most famous one is Salpeter’s
IMF as dN∗ /dM∗ ∝ M∗−2.35 (Salpeter 1955). The low-mass end is ﬂatter than that of M∗ >∼M as−2.7 for 1M <M <10M (Meyer
<
while
dN
/dM
∝
M
∗
M
1M
dN∗ /dM∗ ∝ M ∗−1.2 for 0.1M <
∗
∗
∼ ∗∼

∼ ∗∼ 
et al. 2000). The powers of stellar (−2.7 ∼ −2.35) and prestellar (−2.5 ∼ −2.1) mass functions are
similar. If one prestellar core forms one star, the stellar mass M∗ is proportional to the prestellar
core mass as M∗ = f Mcore and f <∼1, the mass spectrum of prestellar cores completely determines the
IMF.
Figure 1.35 plots the power-law indices against the typical gas densities of respective observations,in which the critical density is taken as a typical density for molecular line studies. The powerlaw index is an increasing function of the typical gas density ngas (Ikeda 2007) and the cores with
5
−3
ngas >∼10 cm have the same power-law index as the IMF. This may indicates these cores are thesites of star formation or direct parents of new-born stars2 .

1.11

Line Width - Size Relation

Larson (1981) compiled the observations for molecular cloud complexes, molecular cloud and molecular clumps published in 1974-1979. He obtained an empirical relation that the size of a structure iswell correlated to the random velocity in the structure which is measured by the width of the emission
line (see Appendix F). Figure 1.37(left) shows this correlation and this is well expressed asσ  1.10km s
2

−1



L
1pc

0.38

,

(1.23)

However, the same authors (Ikeda and Kitamura 2009) obtained n  −2.3 ± 0.3 even for C18 O cores (its criticaldensity ∼ 103 cm−3 ) in OMC-1 region in the Orion A cloud. Thus, the relationship between the core mass function and
IMF has not been cleared yet.

1.11. LINE WIDTH - SIZE RELATION

33

Figure 1.35: The power-law index of n = γ of core mass function dN (M )/dM ∝ M −n . n is plottedagainst typical densities of observed molecular cores. In molecular line observations, the typicaldensities are assumed equal to the critical densities of the transitions of respective molecular line
observations. For mm and sub-mm dust emission, the density is assumed ∼ 106 cm−3 . Taken from
Fig. 8.1 of PhD Thesis of Ikeda (2007).

CHAPTER 1. INTRODUCTION

34

Figure 1.36: Comparison of initial mass functions for ﬁeld stars in the solar neighbourhood. Respective symbols represent S55: Salpeter (1955), MS: Miller & Scalo (1979), Scalo 86: Scalo (1986) and
KTG93: Kroupa, Tout, & Gilmore 1993. Taken from Fig. 1 of Meyer et al (2000).
where σ and L represent respectively the three-dimensional random speed of gas and the size of thestructure. A similar correlation is found only for giant molecular clouds (Sanders, Scoville, & Solomon1985) as

1/2


3
L 0.62−1
ΔVFWHM = (0.535 ± 0.16)km s,(1.24)σ=23 ln 2
1pc
for GMCs whose sizes are larger than 10pc (be careful the typos in their abstract: power was -0.62).
Larson (1981) also found another correlation between the mass M and the random velocity like Figure
1.37 (right), which is expressed asσ  0.43km s−1



M
1M

0.20

.

(1.25)In the next chapter (§2.9), we will see the virial relation, that is, for an isolated system to achievea mechanical equilibrium the gravitational to thermal energy ratio has to be equal to 2:1 for γ = 5/3
gas. The ratio of the gravitational energy ∼ (3/5)GM 2 /(L/2) to the thermal energy M σ 2 /2 is also
ﬁtted as

L 0.14
12 GM

1.1
,
(1.26)
5 σ2 L
1pc
which is weakly dependent of the size or the mass. This seems to mean the ratio is nearly constantirrespective of the mass or size of the clouds.
Since there is a mutual relation between mass, size, and the velocity dispersion to achieve amechanical equilibrium (the Virial relation), there is only one independent correlation in the above twocorrelations (eqs.[1.23] and [1.25]). Although several reasons to explain the correlation are proposed,
we have no consensus yet.

1.11. LINE WIDTH - SIZE RELATION

35

Figure 1.37: The left shows the relation between cloud size (horizontal axis) and the three-dimensional
internal velocity (vertical axis). The right shows a similar correlation between mass and the random
velocity.

36

CHAPTER 1. INTRODUCTION

Chapter 2

Physical Background2.1

Basic Equations of Hydrodynamics

The basic equation of hydrodynamics are (1) the continuity equation of the density [equation (B.11)],∂ρ+ div(ρv) = 0,∂t(2) the equation of motion [equation (B.7)]ρ∂v+ (v · ∇) v = −∇p + ρg,
∂t

(2.1)

(2.2)

and (3) the equation of energy [equation (B.23)]∂+ div( + p)v = ρv · g.
(2.3)∂t
Occasionally barotropic relation p = P (ρ) substitutes the energy equation (2.3). Especially polytropic relation p = KρΓ is often used on behalf of the energy equation. In the case that the gas is
isothermal or isentropic, the polytropic relations of
p = c2is ρ (isothermal)(2.4)

p = c2s ργ (isentropic)

(2.5)

and
are substitution to equation (2.3). [We can replace equation (2.3) with equations (2.4) and (2.5).]

2.2

The Poisson Equation of the Self-Gravity

In this section, we will show the basic equation describing how the gravity works. First, compare
the gravity and the static electric force. Consider the electric ﬁeld formed by a point charge Q at a
distance r from the point source as
1 Q
,
(2.6)
E=
4π 0 r 2
where 0 is the electric permittivity of the vacuum. On the other hand, the gravitational accelerationby the point mass of M at the distance r from the point mass is written down asM
g = −G 2 ,
r
37

(2.7)

CHAPTER 2. PHYSICAL BACKGROUND

38

where G = 6.67 × 10−8 dyn cm2 g−2 is the gravitational constant. Comparing these two, replacing Q
with M and at the same time 1/4π 0 to −G these equations (2.6) and (2.7) are identical with eachother.
The Gauss theorem for electrostatic ﬁeld asρe(2.8)
div E = ,0

and another expression using the electrostatic potential φe as∇2 φe = −

ρe

,

(2.9)

div g = −4πGρ,(2.10)∇2 φ = 4πGρ,(2.11)

0lead to the equations for the gravity as

andwhere ρe and ρ represent the electric charge density and the mass density. Equation (2.11) is called
the Poisson equation for the gravitational potential and describes how the potential φ is determinedfrom the mass density distribution ρ, after the boundary condition for φ is given.
Problem
Consider a spherically symmetric density distribution. Using the Poisson equation, obtain the potential (φ) and the gravitational acceleration (g) for a density distribution shown below.
ρ= ρ0 for r < R= 0 for r ≥ R

Hint: The Poisson equation (2.11) for the spherically symmetric system is


∂φ
1 ∂
r2
r 2 ∂r∂r

2.3

= 4πGρ.

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
 r

Mr =

ρ4πr 2 dr,(2.12)

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

GMr
d2 r
=− 2 .
2dt
r(2.14)

2.3. FREE-FALL TIME

39

Analyzing the equation is straightforward, multiplying dr/dt givesddv 2 /2
=+
dt
dt



GMr
r



,

(2.15)

GMr= E,
r(2.16)

which leads to the conservation of mechanical energy as
1
2



dr
dt

2

−

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(2.17)

because at t = 0, r = R and dr/dt = 0.
The solutions of equation (2.16) are well known as follows:
1. the case of negative energy E < 0. Considering the case that the gas sphere is inﬂowing v < 0,equation (2.16) becomes


1 1/2
dr
1/2 1
= − [2GMr (R)]−,
(2.18)
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 written(2.19)
r = R cos2 η.
In this case, equation (2.18) reduces tocos2 ηdη=dt

GMr (R)
2R3

1/2

.

(2.20)

This gives us the expression of t ast=R3
2GMr (R)1/2 

sin 2η
η+.
2(2.21)

This corresponds to the closed universe in the cosmic expansion (Ω0 > 1).
2. if the energy is equal to zero, the solution of equation (2.16) is written as
r 3/2 − R3/2

2/3

where R = r(t = 0).
Problem
Solve equation (2.16) and obtain (2.22).

=

9GMr (R)2

1/3

t2/3 ,(2.22)

CHAPTER 2. PHYSICAL BACKGROUND

40
1
0.8
0.6
0.4
0.2
0.25

0.5

0.75

1

1.25

1.5

Figure 2.1: Free-fall. x-axis and y-axis represent cos2 η and η + sin 2η/2.
3. If E > 0, the expansion law is given by
t=

R2
2E

1/2



and
x=

Ṙ2 R
2GMr (R)−1 sinh 2η−η
2

(2.23)

E
r = sinh2 η,GMr (R)(2.24)

Ṙ2 GMr (R)
−
> 0.
2
R

(2.25)

where E represents the total energy
E=
Problem
Solve equation (2.16) and obtain equations (2.23) and (2.24).
In the present case, at t = 0, since dr/dt = 0 the energy is negative. Equation (2.19) shows us r
becomes equal to zero (the gas collapses) if η = π/2 as well as η = 0 at t = 0. Equation (2.21)
indicates it occurs at the epoch oft = tﬀ

=


=

1/2

R3
2GMr (R)
3π
32Gρ̄

1/2

,

π
,
2
(2.26)

where ρ̄ represents the average density inside of r, that is Mr /(4πr 3 /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.

2.3.1

Accretion Rate

Equation (2.26) indicates that the gas shell with a large ρ̄ reaches the center earlier than that witha small ρ̄. Immediately, this means a spherical cloud with a uniform density ρ0 contracts uniformly2.4. GRAVITATIONAL INSTABILITY

41

and all the mass reaches the center at t = tﬀ = (3π/32Gρ0 )1/2 . In this case, the mass accretion
rate to a central source becomes inﬁnity at an epoch t = tﬀ . In contrast, consider a cloud whose
density gradually decreases outwardly. In this case, the outer mass shell has smaller ρ̄ than the innermass shell. Therefore even when the inner mass shell collapses and reaches the center, the outermass shell are contracting and does not reach the center. This gives a smaller mass accretion rate
than a uniform cloud. If the gas pressure is neglected, the accretion rate is determined by the initial
spatial distribution of the density. We will compare the accretion rate derived here with results ofhydrodynamical calculation in §4.7

2.4

Gravitational Instability

Here, we will study a typical size where the self-gravity plays an important role and forms density
inhomogeneities with the Jean wavelength.

2.4.1

Linear Analysis

Consider a uniform gas with density ρ0 and pressure p0 without motion u0 = 0. In this uniform gas
distribution, we assume small perturbations. As a result, the distributions of the density, the pressure
and the velocity are perturbed from the uniform ones asρ = ρ0 + δρ,(2.27)

p = p0 + δp,(2.28)

u = u0 + δu = δu,(2.29)

andwhere the amplitudes of perturbations are assumed much small, that is, |δρ|/ρ0  1, |δp|/p0  1
and |δu|/cs  1. We assume the variables changes only in the x-direction. In this case the basic
equations for isothermal gas are∂ρ ∂ρu+= 0,(2.30)
∂t
∂x


∂u
∂p
∂u
+u
+ ρgx ,=−
(2.31)
ρ
∂t
∂x
∂x
and
p = c2is ρ,(2.32)

where u and gx represent the x-component of the velocity and that of the gravity, respectively.
Using equations (2.27), (2.28), and (2.29), equation (2.30) becomes∂ρ0 + δρ ∂(ρ0 + δρ)(u0 + δu)
+
= 0.
∂t
∂x

(2.33)

Noticing that the amplitudes of variables with and without δ are completely diﬀerent, two equations
are obtained from equation(2.33) as∂ρ0 ∂ρ0 u0+∂t∂x∂δρ ∂ρ0 δu + δρu0
+
∂t
∂x

= 0,

(2.34)= 0,(2.35)

CHAPTER 2. PHYSICAL BACKGROUND

42

where the above is the equation for unperturbed state and the lower describes the relation between thequantities with δ. Equation (2.34) is automatically satisﬁed by the assumption of uniform distribution.
In equation (2.35) the last term is equal to zero. Equation of motion∂u0 + δu∂u0 + δu+ (u0 + δu)(ρ0 + δρ)
∂t
∂x



=−

∂p0 + δp
∂φ0 + δφ+ (ρ0 + δρ),∂x∂x(2.36)

gives the relationship between the terms containing only one variable with δ as follows:
ρ0

∂δp
∂δu
∂δφ
=−
− ρ0
.
∂t
∂x
∂x

(2.37)The perturbations of pressure and density are related with each other as follows: for the isothermal
gas
 
δp
p0
∂p
=
=
= c2is ,(2.38)
δρ∂ρ T
ρ0and for isentropic gas
 
p0
∂p
δp
=
=γ
= c2s .
(2.39)
δρ
∂ρ ad
ρ0

2.4.2

Sound Wave

If the self-gravity is ignorable, equation (2.35)
∂δu
∂δρ
+ ρ0
= 0,
∂t
∂x

(2.40)and equation (2.37)∂δu∂δρ= −c2is,(2.41)∂t∂x
where we assumed the gas is isothermal, these two equations describe the propagation and growth of
perturbations. If the gas acts adiabatically, replace cis with cs .
Making ∂/∂x×(2.40) and ∂/∂t×(2.41) vanishes δρ and we obtainρ02∂ 2 δu2 ∂ δu
−
c
= 0.
is
∂t2
∂x2

(2.42)

Since this leads to
∂δu
∂δu
− cis
∂t
∂x
∂δu
∂δu
+ cis
∂t
∂x

= 0,

(2.43)= 0,(2.44)

equation (2.42) has a solution that a wave propagates with a phase velocity of ±cs . Since thedisplacement (∝ δu) is parallel to the propagation direction x, and the restoring force comes fromthe pressure, this seems the sound wave.
ProblemInterstellar gas contains mainly Hydrogen and Helium, whose number ratio is approximately 10:1.
Obtain the value of average molecular weight for the fully ionized interstellar gas with temperature
T = 104 K (components are ionized H+ (HII) and He+2 (HeIII) and electron e−1 ). How about themolecular gas (T = 10K) containing molecular H2 , neutral He (HeI) and no electron?

2.5. JEANS INSTABILITY

2.5

43

Jeans Instability

Sound wave seems to be modiﬁed in the medium where the self-gravity is important. Beside thecontinuity equation (2.35)
∂δu
∂δρ
+ ρ0
= 0,
(2.45)∂t∂x
and the equation of motion (2.37)
ρ0

∂δu
∂δρ
∂δφ
= −c2is
− ρ0
,
∂t
∂x
∂x

(2.46)∂ 2 δφ= 4πGδρ,∂x2(2.47)

the linearized Poisson equation

should be included. ∂/∂x× (2.46) gives
ρ0∂ 2 δu∂x∂t= −c2is∂ 2 δρ
− 4πGρ0 δρ.
∂x2(2.48)

where we used equation (2.47) to eliminate δφ. This yields
2∂ 2 δρ
2 ∂ δρ
=c+ 4πGρ0 δρ.
is
∂t2
∂x2

(2.49)

where we used ∂/∂t× (2.45).
This is the equation which characterizes the growth of density perturbation owing to the selfgravity. Here we consider the perturbation are expressed by the linear combination of plane wavesas
A(ω, k) exp(iωt − ikx),(2.50)
δρ(x, t) =where k and ω represent the wavenumber and the angular frequency of the wave. Picking up a plane
wave and putting into equation (2.49), we obtain the dispersion relation for the gravitational
instability as(2.51)
ω 2 = c2is k2 − 4πGρ0 .
Reducing the density to zero, the equation gives us the same dispersion relation as that of the soundkJ = (4πGρ0 )1/2 /cis ), since ω 2 > 0 the wave is ordinary
wave as ω/k = cis . For short waves (koscillatory wave. Increasing the wavelength (decreasing the wavenumber), ω 2 becomes negative fork < kJ = (4πGρ0 )1/2 /cis . For negative ω 2 , ω can be written ω = ±iα using a positive real α. In thiscase, since exp(iωt) = exp(∓αt), the wave which has ω = −iα increases its amplitude exponentially.
This means that even if there are density inhomogeneities only with small amplitudes, they grow ina time scale of 1/α and form density inhomogeneities with large amplitudes.
The critical wavenumber(2.52)
kJ = (4πGρ0 )1/2 /cis
corresponds to the wavelength
2π
=
λJ =
kJ

πc2is
Gρ0

1/2

,

(2.53)

which is called the Jeans wavelength. Ignoring a numerical factor of the order of unity, it is shownthat the Jeans wavelength is approximately equal to the free-fall time scale (eq.[2.26]) times theCHAPTER 2. PHYSICAL BACKGROUND

44

Figure 2.2: Dispersion Relation

Figure 2.3: Thin disk.
sound speed. The short wave with λ  λJ does not suﬀer from the self-gravity. For such a scale, theanalysis we did in the preceding section is valid.
Typical values in molecular clouds, such as cis = 200m s−1 , ρ0 = 2 × 10−20 g cm−3 , give us theJeans wavelength as λJ = 1.7 × 1018 cm = 0.56pc. The mass contained in a sphere with a radiusr = λJ /2 is often called Jeans mass, which gives a typical mass scale above which the gas collapses.
Typical value of the Jeans mass is as follows


MJ 

4π
λJ
ρ0
3
2

3

=

π
6



π
Gρ0

3/2

c3is ρ0 .

(2.54)

Using again the above typical values in the molecular clouds, cis = 200m s−1 , ρ0 = 2 × 10−20 g cm−3 ,the Jeans mass of this gas is equal to MJ  27M .

2.6

Gravitational Instability of Thin Disk

Disks are common in the Universe. Spiral and barred spiral galaxies have disks where stars areformed. In more small scale, gas and dust disks are often found around protostars. Moreover, such

2.6. GRAVITATIONAL INSTABILITY OF THIN DISK

45

a disk may become a proto-planetary disk. It is valuable to study how the self-gravity works in suchthin structures. Here, we assume a thin disk extending in x- and y-directions whose surface density∞ρdz, in other word the density is written using the Dirac’s delta function δ(z) as
is equal to σ = −∞ρ(x, y, z) = σ(x, y)δ(z).

(2.55)

Integrating along the z-direction basic equations (2.45) and (2.46),the linearized basic equations forthe thin disk are as follows:
∂δu
∂δσ
+ σ0
= 0,
(2.56)
∂t
∂x
∂δu
∂δσ
∂δφ
= −c2is− σ0,(2.57)σ0∂t∂x∂x∂ 2 δφ ∂ 2 δφ+
= 4πGδσδ(z),(2.58)∂x2∂z 2
where we assumed σ = σ0 + δσ, u = δu, φ = φ0 + δφ and took the ﬁrst order terms (those contain
only one δ).
Outside the disk, the rhs of equation (2.58) is equal to zero. It reduces to the Laplace equation∂ 2 δφ ∂ 2 δφ
+
= 0.
∂x2
∂z 2

(2.59)

δX(x, t) = δA exp(iωt − ikx),(2.60)Taking a plane wave of
equation (2.59) is reduced to∂ 2 δφ− k2 δφ = 0.
∂z 2This has a solution which does not diverge at the inﬁnity z = ±∞ asδφ = δφ(z = 0) exp(−k|z|).

(2.61)(2.62)On the other hand, integrating equation (2.58) from z = −0 to z = +0 or in other word, applying theGauss’ theorem to the region containing the z = 0 surface, it is shown that the gravity δgz = −∂δφ/∂z
has a jump crossing the z = 0 surface as∂δφ 
∂δφ 
−= 4πGδσ.
∂z z=+0∂z z=−0(2.63)

Equations (2.62) and (2.63) lead a ﬁnal form of the potential asδφ = −2πGδσexp(−k|z|).
k(2.64)

Putting this to equations (2.57), and using equations (2.56) and (2.57), we obtain the dispersion
relation for the gravitational instability in a thin disk asω 2 = c2is k2 − 2πGσ0 k.

(2.65)

This reduces to the dispersion relation of the sound wave for the short wave k
2πGσ0 /c2is . While2
for a longer wave than λcr = cis /Gσ0 , an exponential growth of δσ is expected. The dispersion
relation is shown in Fig.2.2.

CHAPTER 2. PHYSICAL BACKGROUND

46

2.6.1

Rotating Disk

In the case of a rotating disk with an angular rotation speed Ω, basic equations are

σ0

∂δu
∂δσ
+ σ0
= 0,
∂t
∂x

(2.66)∂δu= −c2is ∇δσ − σ0 ∇δφ − 2σ0 Ω × δu,∂t(2.67)∂ 2 δφ− k2 δφ = 4πGδσδ(z),(2.68)∂z 2
where the third term of r.h.s. of equation (2.67) represents the Colioris force. Equation (2.68) gives
a solution identical to equation (2.64) as
δφ = −

2πGδσ
exp [i(ωt − kx) − k|z|] .
k(2.69)

Choosing a direction in which the wavenumber vector can be expressed as k = (k, 0, 0), we can reduce
equations (2.66) and (2.67) to(2.70)iωδσ − ikσ0 δux = 0,iωδux = −c2is (−ik)δσ− 2πiGδσ + 2Ωδuy ,σ0(2.71)iωδuy = −2Ωδux .

(2.72)

These three equations together with equation (2.69) bring us a dispersion relation asω 2 = c2is k2 − 2πGσ0 k + 4Ω2 .

(2.73)

Comparing with equation(2.65), this indicates rotation works to stabilize the system.
Equation (2.73) is rewritten asω2= c2is

πGσ0
k− 2
cis

2





(πGσ0 )2 4Ω2 c2is+
−1 .
(πGσ0 )2c2is

(2.74)

DeﬁningQ ≡ 2Ωcis /πGσ0 ,(2.75)

we can see that ω 2 > 0 for all wavenumbers if Q > 1 and that if Q < 1 for some range of wavenumberω 2 becomes negative. A rotating disk with Q < 1 is unstable for some range of wavenumber. This
number is called Toomre’s Q. This is useful to see whether a galactic disk is stable or not. For thegalactic disk, Toomre’s Q must be modiﬁed asQ ≡ κcis /πGσ0(2.76)

where κ represents the epicyclic frequency as
dΩ2
+ 4Ω2
κ≡ R
dR
See section 3.5 for the galactic disk.

1/2.

(2.77)

2.7. CONVECTIVE INSTABILITY

47zρ∗(z+Δz)p*(z+Δz)

Δz

ρ∗(z)
p*(z)

Figure 2.4: Convection.

2.7

Convective Instability

If water is heated from the bottom and temperature diﬀerence between the top and the bottom
exceeds a limit, convection is driven. Water heated from the bottom climbs and cool water on the top
descends. This transfers the thermal energy from the bottom to the top. In this section, we describe
the condition in which convection is driven.
Consider a hydrostatic balanced atmosphere in which the hydrostatic balance equation is satisﬁed:
−

dφ
dp
dp
−ρ
= − + ρg = 0,
dz
dz
dz

(2.78)

where we assumed v = 0 in equation (2.2) and the gravity is working downwards in z-direction (g < 0)
The pressure and density of the atmosphere are p(z) and ρ(z). We consider a gas element (hatchedregion in Fig.2.4), whose density ρ∗ and p∗ are equal to those of the atmosphere ρ(z) and p(z) asρ∗ (z) = ρ(z),p∗ (z) = p(z).

(2.79)

Further, we assume this gas element to move from z to z + Δz adiabatically, that is,p∗ (z)p∗ (z + Δz).
γ =[ρ∗ (z + Δz)][ρ∗ (z)]γ(2.80)

Pressure balance is required between the pressures of the gas element at z + Δz, p∗ (z + Δz), andthe atmosphere p(z + Δz). If the density of the gas element at z + Δz, ρ∗ (z + Δz), is smaller thanthat of the atmosphere, the gravity force in equation (2.78) is weaker than the pressure force and theelement keeps climbing. Thus the condition for the convective instability is written asρ∗ (z + Δz) < ρ(z + Δz).

(2.81)Usingp(z + Δz) = p(z) +

dp
Δz
dz

and ρ(z + Δz) = ρ(z) +

dρ
Δz,
dz

(2.82)

CHAPTER 2. PHYSICAL BACKGROUND

48

Figure 2.5: Left: Explanation of Laval nozzle. Right: The relation between the cross-section S(x)and the ﬂow velocity vx .
we can rewrite equation (2.81) into the relation in the variables of the atmosphere (variables without*) asd ln ρ
d ln p<γ(2.83)dzdzord ln(p/ργ )<0(2.84)dz
This means the speciﬁc entropy s = ln(p/ργ ) + K decreases upwardly. Thus, if we consider the
adiabatic process, the atmosphere in which a speciﬁc entropy decreases upwardly is unstable for theconvection.
ProblemObtain equation (2.83) from equation (2.81) using equation (2.82).

2.8

Super- and Subsonic Flow

Flow whose velocity is faster than the sound speed is called supersonic, while that slower than thesound speed is called subsonic. The subsonic and supersonic ﬂows are completely diﬀerent.

2.8.1

Flow in the Laval Nozzle

Consider a tube whose cross-section, S(x), changes gradually, which is called Laval nozzle. Assuming
the ﬂow is steady ∂/∂t = 0 and is essentially one-dimensional, the continuity equation (2.1) isrewritten as
ρuS = constant,(2.85)
or

1 ∂u
1 ∂S
1 ∂ρ
+
+
= 0.
ρ ∂x u ∂x S ∂x(2.86)

1 ∂p
c2 ∂ρ
∂u
=−
=− s
,
∂x
ρ ∂x
ρ ∂x

(2.87)

Equation of motion (2.2) becomesu2.8. SUPER- AND SUBSONIC FLOW

49

where we used the relationship of
∂p
=
∂x



∂p
∂ρ



∂ρ
∂ρ
= c2s .
∂xs ∂x

(2.88)

When the ﬂow is isothermal, use the isothermal sound speed c2is instead of the adiabatic one. Fromequations (2.86) and (2.87), we obtain
1 ∂S
u2
1 ∂u
=
,
−1
2
cs
u ∂x
S ∂x

(2.89)

where the factor M = u/cs is called the Mach number. For supersonic ﬂow M > 1, while M < 1
for subsonic ﬂow.
In the supersonic regime M > 1, the factor in the parenthesis of the lhs of equation (2.88) is
positive. This leads to the fact that the velocity increases (du/dx > 0) as long as the cross-section
increases (dS/dx > 0). On the other hand, in the subsonic regime, the velocity decreases (du/dx < 0)while the cross-section increases (dS/dx > 0). See right panel of 2.5.
If M = 1 at the point of minimum cross-section (throat), two curves for M < 1 and M > 1 havean intersection. In this case, gas can be accelerated through the Laval nozzle. First, a subsonic ﬂow
is accelerated to the sonic speed at the throat of the nozzle. After passing the throat, the gas follows
the path of a supersonic ﬂow, where the velocity is accelerated as long as the cross-section increases.

2.8.2

Steady State Flow under an Inﬂuence of External Fields

Consider a ﬂow under a force exerted on the gas whose strength varies spatially. Let g(x) represent
the force working per unit mass. Assuming the cross-section is constantρu = constant,(2.90)

immediately we have

1 ∂u
1 ∂ρ
+
= 0.
ρ ∂x u ∂x
On the other hand, the equation motion is
u

c2 ∂ρ
∂u
=− s
+ g(x),∂xρ ∂x(2.91)(2.92)

From equations (2.91) and (2.92), we obtain
g(x)
u2
1 ∂u
= 2 .
−1
2
cs
u ∂x
cs

(2.93)

Consider an external ﬁeld whose potential is shown in Figure 2.6(Left). (1) For subsonic ﬂow, thefactor in the parenthesis is negative. Before the potential minimum, since g(x) > 0, u is decelerated.
On th other hand, after the potential minimum, u is accelerated owing to g(x) < 0. Using equation(2.90), this leads to a density distribution in which density peaks near the potential minimum. (2)For supersonic ﬂow, the factor is positive. In the region of g(x) > 0, u is accelerated. After passing
the potential minimum, u is decelerated. The velocity and the density distribution is shown in Figure
2.6(right-lower panel).
The density distribution of the subsonic ﬂow in an external potential is similar to that of hydrostatic one. That is, considering the hydrostatic state in an external potential, the gas density peaks at

50

CHAPTER 2. PHYSICAL BACKGROUND

Figure 2.6: Left: External potential ﬁeld. Right: Velocity and density variations. Gas ﬂows in theexternal ﬁeld whose potential is shown in the left panel. The upper panel represents a subsonic ﬂow.
The lower panel does a supersonic ﬂow.

2.8. SUPER- AND SUBSONIC FLOW

51

the potential minimum. On the other hand, The density distribution of the supersonic ﬂow looks likethat made by ballistic particles which are moving freely in the potential. Owing to the conservationof the total energy (kinetic + potential energies), the velocity peaks at the potential minimum. And
the condition of mass conservation leads to the distribution in which the density decreases near the
potential minimum.

2.8.3

Stellar Wind — Parker Wind Theory

Stellar winds are observed around various type of stars. Early type (massive) stars have large luminosities; the photon absorbed by the bound-bound transition transfers its outward momentum tothe gas. This line-driven mechanism seems to work around the early type stars. On the other hand,acceleration mechanism of less massive stars are thought to be related to the coronal activity or dust
driven mechanism (dusts absorb the emission and obtain outward momentum from the emission).
Here, we will see the identical mechanism in §2.8.1 and §2.8.2 works to accelerate the wind froma star. Consider a steady state and ignore ∂/∂t = 0. The continuity equation (2.1) givesr 2 ρu = const,
where we useddivρv =

(2.94)

1 ∂r 2 ρu .
2r ∂r(2.95)

This leads to

2 1 ∂ρ 1 ∂u++
= 0.
r r ∂r u ∂r
The equation of motion (2.2) is as follows:u(2.96)c2 ∂ρ GM∗
∂u
=− s
− 2 ,
∂r
ρ ∂r
r

(2.97)

where we used g = −GM∗ /r 2 (M∗ means the mass of the central star). From these two equations(2.96) and (2.97) we obtain 1 ∂u2 GM∗ 1= − 2 2,M2 − 1(2.98)u ∂rr
cs r
where M ≡ u/cs represents the Mach number of the radial velocity. Take notice that this is similar
to equations (2.89) and (2.93). That is, the fact that the rhs of equation (2.98) is positive corresponds
to increasing the cross-section dS/dx > 0. On the contrary, when the rhs is negative, the ﬂuid acts
as the cross-section S is decreasing.
For simplicity, we assume the gas is isothermal. The rhs of equation (2.98) varies shown in Figure
2.7. Therefore, near to the star, the ﬂow acts as the cross-section of nozzle is decreasing and far
from the star it does as the cross-section is increasing. This is the same situation that the gas ﬂows
through the Laval nozzle.
Using a normalized distance x ≡ r/(GM∗ /2c2is ), equation (2.98) becomesM2 − 1This is rewritten asddx

 1 ∂M

M ∂x

=

2
2
− 2.
x x

2
M2
− log M − 2 log x −2x(2.99)= 0,(2.100)CHAPTER 2. PHYSICAL BACKGROUND

52

Figure 2.7: Right-hand side of equation (2.98) is plotted against the distance from the center r.
we obtain the solution of equation (2.99) as4+ C.
(2.101)xThis gives how the Mach number M varies changing x. To explore this, we deﬁne two functions:
M2 − 2 log M = 4 log x +f (M) = 2M2 − 2 log M(2.102)is a function only depending on M andg(x) = 4 log x +4+Cx(2.103)is a function only depending on x.
Since the minima of f (M) and g(x) are respectively 1 and C + 4, the permitted region in (x, M)changes for the values of C.
1. If C = −3, for all values of x > 0 there exist M which satisﬁes f (M) = g(x). This corresponds
to the two curves which pass through a critical X-point of (x, M) = (1, 1) in Figure 2.9.
2. If C < −3, the minimum of g(x) is smaller than that of f (M). In this case, for x whereg(x) < 1 = min(f (M)), there is no solution. Thus, f (M) = g(x) has solutions for x < x1
and x > x2 , where x1 < x2 , and g(x1 ) = g(x2 ) = 1. This corresponds to the curves running
perpendicularly in Figure 2.9.
3. If C > −3, the minimum of f (M) is smaller than that of g(x). In this case, f (M) = g(x) has
solutions for M < M1 and M > M2 , where M1 < M2 , and f (M1 ) = f (M2 ) = C + 4. This
corresponds to curves running horizontally in Figure 2.9.
Out of the two solutions of C = −3, a wind solution is one with increasing M while departingfrom the star. This shows us the outﬂow speed is slow near the star but it is accelerated and asupersonic wind blows after passing the critical point. Since the equations are unchanged even if we
replace u with −u, the above solution is valid for an accreting ﬂow u < 0. Considering such a ﬂow,
the solution running from the lower-right corner to the upper-left corner represents the accretion ﬂow,
in which the inﬂow velocity is accelerated reaching the star and ﬁnally accretes on the star surface
with a supersonic speed.

2.8. SUPER- AND SUBSONIC FLOW

53

fHML
2

gHxL
2

1.75

1.75

1.5

1.5

1.25

1.25

1

1

0.75

0.75

0.5

0.5

0.25

0.25
0.5

1

1.5

2

M

0.5

1

1.5

2

x

Figure 2.8: f (M) (left) and g(x) for C = −3 (right)

Figure 2.9: M vs x (the x-axis is the normalized distance from the center x ≡ r/(GM∗ /2c2is ) and they-axis is the Mach number M.)

CHAPTER 2. PHYSICAL BACKGROUND

54

2.9

Virial Analysis

For a system to achieve a mechanical equilibrium, a relation must be satisﬁed between energies
such as potential, thermal and kinetic energies. This is called the Virial relation. For example,
a harmonic oscillator md2 x/dt2 = −kx has a potential energy of WH = (1/2)kx2 and a kinetic
energy KH = (1/2)mẋ2 . Averaging these two energies over one period, both energies give thesame absolute value proportional to the oscillation amplitude a squared as < WH >= (1/4)ka2 and< KH >= (1/4)ka2 . Another example is a Kepler problem. For simplicity, consider mass m runningon a circular orbit with a radius a from a body with a mass M . The gravitational and kinetic energies
are equal to WK = −GmM/a and KK = (1/2)mv 2 = (1/2)GmM/a, where we used the centrifugalbalance as mv 2 /a = GM m/a2 . As a result, for the harmonic oscillator < WH >:< KH >= 1 : 1
while for the circular Kepler problem WK : KK = −2 : 1. This ratio is known to be related to thepower n of the potential as φ ∝ xn . Important nature of the self-gravity is understood only with thisrelation without solving the hydrostatic balance equations. In the following, we describe the Virial
relation satisﬁed with isolated systems such as stars.
Hydrodynamic equation of motion using the Lagrangean time derivative [eq.(B.3)] isρdv= −∇p − ρ∇Φ.
dt

(2.104)For simplicity, consider a spherical symmetric conﬁguration. The equation of motion is expressed as
ρ

∂p
∂Φ
dv
=−
−ρ
.
dt
∂r
∂r

(2.105)Multiplying radius r to the equation and integrating by the volume dV = 4πr 2 dr over a volume fromr = 0 to r = R, we obtain the Virial relation as
1 d2 I− 2T = 3(γ − 1)U + W,
2 dt2
where




2

r 2 dM,ρr dV =

I=

(2.106)(2.107)

is an inertia of this body, and T and U are, respectively, the kinetic and thermal energies as

1 2
ρv dV =
2

T =



U=
Further,
W =−



ρ



1 2
v dM,
2

(2.108)

p
dV.
γ−1

GMr
dV = −r(2.109)

GMr
dM,
r


is a gravitational energy, where dM = 4πr 2 ρdr and Mr = 0r dM =Since
1 d2 r 2
dv=− u2 ,rdt
2 dt2(2.110)r

2
0 4πr ρdr

(2.111)

2.9. VIRIAL ANALYSIS

55

the lefthand-side of equation (2.105) is rewritten as
 R

ρr
0

dv
4πr 2 dr =
dt
=
=

 R

1 d2 r 2− u2 ρ4πr 2 dr,
2 dt2

0



1 d2 R 2
r dm −
2 dt2 01 d2 I
− 2T,
2 dt2

 R

v 2 dm0(2.112)using equations (2.107) and (2.108).
On the other hand, the ﬁrst term of the rhs of equation (2.105) becomes
−

 R
0



 R
R

dp
r 4πr 2 dr = − 4πr 3 p − 3
4πr 2 pdr ,
0
dr
0
 R

pdV,

= 3
0

= 3(γ − 1)U.

(2.113)

To derive this equation, we have assumed the pressure diminishes at a radius r = R and the surfacepressure term does not appear in the ﬁnal expression. This is valid for an isolated system such as astar.
The last term of the rhs of equation (2.105) is written as
−

 R
0



R GM
dΦ
r
rρ 4πr 2 dr = −ρ 2 4πr 2 dr,drr0(2.114)

where we used equation (2.13). This is rewritten as
−

 R



ρ
0

R
GMr
GMr
dM
4πr 2 dr =
−
2
r
r
0
= W.

(2.115)The energy +GMr /r per unit mass is necessary for a gas element to move from the radius r, insidewhich mass Mr is contained, to the inﬁnity. Adding the energy +GMr /rdM for all the gas, thepotential energy is obtained. In the case of a star composed of uniform density ρ0 ,W =−3 GM 2,5 R(2.116)

where M = 4πR3 ρ0 /3.
To obtain a condition for the mechanical equilibrium, we assume d2 I/dt2 = 0. Equation (2.106)
becomes
2T + 3(γ − 1)U + W = 0.
(2.117)
Assuming the system is static v = 0, the above equation reduces to3(γ − 1)U + W = 0.

(2.118)In the case of γ = 5/3 this reduces to 2U + W = 0. The total energy E = U + W is expressed asE= (4 − 3γ)U4−3γ
W
= − 3(γ−1)= −U (for γ = 5/3)
1
= W (for γ = 5/3).
2(2.119)

CHAPTER 2. PHYSICAL BACKGROUND

56

E

E

W

W

ΔE
ΔU

U

U
ΔW

Figure 2.10: Virial relation for stars composed of γ = 5/3 gas.
For the gas with γ > 4/3, equation (2.119) gives a negative total energy E < 0 and the system is in
a conﬁned state. However, if γ < 4/3, the gravity can not conﬁne the gas.
For γ = 5/3 > 4/3, equation (2.118) gives
2U + W = 0.

(2.120)This shows us a strange nature of the self-gravitating gas (see Fig. 2.10). That is, if the heatﬂux ﬂows outward, the total energy decreases ΔE < 0. The system must contract and equation(2.119) indicates ΔW = 2ΔE < 0 (the gravitational energy decreases: the absolute value of thegravitational energy increases). However, for the thermal energy, equation (2.120) indicates that
ΔU = −ΔE = −(1/2)ΔW > 0 for this system to be static. This shows that if the heat ﬂux ﬂows
out from the system the thermal energy increases in the self-gravitating system. This comes from thecontraction due to the gravity.

2.9.1

Virial Mass

Consider γ = 5/3 gas, the gravitational energy is equal to twice of the kinetic+thermal energy[eq.(2.120)] in the equilibrium state. This gives a deﬁnition of Virial mass. The kinetic energy is
obtained from the 3D velocity dispersion σ3D as
1
2
U = M σ3D2

(2.121)for globular clusters and elliptical galaxies, where M is the mass of the object. For molecular cloud
this is estimated from the width of the emission line. Equation (F.7) of Appendix F gives
2
=
σ3D

3 2
v
.
2 ln 2 x,HWHM(2.122)

2.10. RADIATIVE TRANSFER

57

where vx,HWHM is the half width of the half maximum of the emission line (if the radiation is opticallythin; see Fig.F.1). Since the gravitational energy is proportional to M 2 as
W =−

GM 2
,
R

(2.123)the Virial theorem gives

2
Rσ3D
G
if the surface term is ignored. This gives an estimation of mass of astronomical object.

M=

2.10

Radiative Transfer

2.10.1

Radiative Transfer Equation(2.124)

Speciﬁc intensity of radiation ﬁelds I(x, n, ν) is deﬁned by the radiation energy dE transferred byphotons with wavelength ν through a unit surface placed at x whose normal is directed to n per unit
time per unit wavelength, and per unit steradian. Average intensity of radiation J(x, ν) is deﬁned asJ(x, ν) =I(x, n, ν)

dΩ
,
4π

(2.125)where J is obtained by averaging I for the solid angle. This is related to the energy density ofradiation u(x, ν) as
4π
J(x, ν).
(2.126)u(x, ν) =c
If the radiation is absorbed in the displacement δs as δI(x, n, ν), δI must be proportional to δs
and I(x, n, ν) as
δI(x, n, ν) = −κI(x, n, ν)δs,(2.127)
where κ is a coeﬃcient and called volume absorption coeﬃcient. The dimension of κ is cm−1 . We
can rewrite the above to the diﬀerential equation as
dI(x, n, ν)= −κI(x, n, ν) = −κρ ρI(x, n, ν),(2.128)ds
where we used mass absorption coeﬃcient κρ which represents the absorption per mass. Equation(2.127) is reduced todIdI(x, n, ν)≡= −I(x, n, ν),(2.129)
κds
dτ
where
 s
κδs,
(2.130)τ=0

is called the optical depth. This means that τ is a measure for absorption as the intensity decreases
at a factor exp(−1) from τ = 0 to τ = 1.
If the ray runs δs crossing a volume whose volume emissivity equal to the intensity increases
δI(x, n, ν) = (x, ν)δs.

(2.131)The volume emissivity (x, ν) is the energy emitted by a unit volume at a position x per unit time
per unit solid angle per unit wavelength. From equations (2.127) and (2.131), the radiation transfer
is written as
dI(x, n, ν)= −κI(x, n, ν) + (x, ν).
(2.132)ds

CHAPTER 2. PHYSICAL BACKGROUND

58

Level

Energynu

Enuhν=Enu-En

nEn

hν=En-End
ndEnd

Figure 2.11:
Using the optical depth τ =s

0 κds, this givesdI(ν)= −I(ν) + S(ν),(2.133)dτ
where S(ν) is called the source function and is deﬁned as S(ν) = (ν)/κ(ν). Assuming the speciﬁc
intensity I = I0 at the point of τ = 0, that at the point τ on the same ray is given τ (x)I = I0 exp (−τ ) +

0





S(s(τ  )) exp τ  − τ (s) dτ  .

(2.134)In the case of constant source term S =const it reduces to
I = I0 exp (−τ ) + S (1 − exp (−τ )) .

(2.135)Equation (2.134) gives
I

S ...
τ
1,
τ S + I0 (1 − τ ) . . . τ  1.

(2.136)This indicates that if we see an optically thick cloud τ
1 the speciﬁc intensity reaches us represents
S, while if we see an transparent cloud τ  1, I represents that of background.
ProblemShow that equations (2.134) and (2.134) are solutions of equation (2.133).

2.10.2Einstein’s Coeﬃcients

In this section, we describe the absorption coeﬃcient κ and the emissivity for the line emission
and absorption. Here, consider a hypothetical atom which has only three levels named nd , n, nu .
and assume the energy of each level is, respectively, End , En , and Enu . When the atom changes
its state from nu to n, an photon with energy of hνnu →n = Enu − En is emitted, if this transitionis radiative. Besides the radiative process, atom changes its state exchanging its energy with otheratoms, molecules and electrons.
Spontaneous Emission
Atoms in an excited state nu changes its state from nu to n spontaneously. The number of transitions
is proportional to the number of atoms in the nu state Nnu as
Anu →n Nnu ,(2.137)

2.10. RADIATIVE TRANSFER

59

where Anu →n represents the probability of transition from nu to n in unit time. This coeﬃcient iscalled Einstein’s A coeﬃcient. The lifetime of nu level is related to this coeﬃcient asτlife(nu ) =1,
Anu →n + Anu →nd(2.138)for this three-level atom.
Absorption
Photon with energy Enu − En is absorbed by the atom in the state of n. At the same time the atom
is excited from the state n to nu . The number of this transition is proportional to the number ofatoms in the state n, Nn , the number of photons passing and a cross-section of this absorption. Thisis written as(2.139)Bn→nu Jνn→nu Nn ,where J represents the average speciﬁc intensity of radiation with corresponding wavelength. This
coeﬃcient is called Einstein’s B coeﬃcient for absorption.
Induced Emission
When atoms in the state nu are in the radiation ﬁeld with the average intensity of Jnu →n , downward
transition from nu to n is induced and atoms emit photons with energy hνnu →n . This is called theinduced emission or stimulated emission. This rate is expressed as
Bnu →n Jnu →n Nnu .

(2.140)Collisional Excitation and Deexcitation
Beside radiative processes, energy states of atoms are changed by collisions with other atoms andmolecules, which are called collisional excitation and deexcitation. Since such transitions are proportional to the density of atoms or molecules, the number of transitions is written as follows: for
excitation from n to nu
Cn→nu Nn nH ,(2.141)and deexcitation from nu to n
Cnu →n Nnu nH ,(2.142)where the coeﬃcients C are determined by laboratory experiments and quantum mechanics calculations.

2.10.3

Relation of Einstein’s Coeﬃcients to Absorption and Emissivity

Emissivity and absorption coeﬃcients shown in §2.10.1 are related to the Einstein’s coeﬃcients in§2.10.2. Net absorption is deﬁned as the absorption minus the induced emission. Energy absorbed inunit volume is expressed as follows: integrating equation (2.127) over the solid angle, the ﬁrst termof the rhs becomesκ(ν)I(ν)dνdΩ = hνnu →n (Bn→nu Nn − Bnu →n Nnu ) Jnu →n .

(2.143)

CHAPTER 2. PHYSICAL BACKGROUND

60
Assuming isotropic radiation ﬁeld, since

IdΩ = 4πJ,(2.144)the absorption coeﬃcient integrated over the frequency range of one transition is written asκ(ν)dν =1hνnu →n (Bn→nu Nn − Bnu →n Nnu ) .
4π(2.145)

As for the emission, the second term of equation (2.127) becomes (ν)dνdΩ = hνnu →n Anu →n Nnu .
Thus,



(ν)dν =

(2.146)hνnu →nAnu →n Nnu .
4π(2.147)

Relation between Einstein’s Coeﬃcients∗
Analysis using the Einstein’s coeﬃcients is applicable to the state in which the thermal equilibrium is
not achieved. However, if this is applied to the thermal equilibrium, we can obtain relations between
these coeﬃcients.
Assume the detailed balance is achieved between the spontaneous emission, the induced emission
and the absorption as
Nnu (Anu →n + Bnu →n Jnu →n ) = Nn Bn→nu Jnu →n .

(2.148)In the thermal equilibrium, the population of atoms is given by the Boltzmann distribution as




gn
En n
Nn u
= u exp − u
,
Nn
gn
kT
where gn represents the statistical weight for the state n.
Equation (2.148) gives(2.149)Nnu
Nnu Anu →n= Jnu →n Bn→nu −Bnu →n .
Nn
Nn(2.150)In the case of thermal equilibrium, the average speciﬁc intensity agrees with the Planck function
B(ν, T ) as
1
2hν 3
,
(2.151)
J(ν) = B(T, ν) = 2
c exp (hν/kT ) − 1
Using equations (2.149) and (2.151), equation (2.148) reduces toB
exp (hν/kT ) − ggnnu Bnu →n
2hν 3
gnu
n→nu
.
Anu →n = 2 Bn→nugncexp (hν/kT ) − 1(2.152)

Lefthand-side of the equation is not dependent on the temperature. For this equation to be valid forvarious temperature,gnu Bnu →n= 1,(2.153)gn Bn→nu
and
2hν 3
gnu
Anu →n = 2 Bn→nu .
(2.154)gnc2.10. RADIATIVE TRANSFER

61

Non-Local Thermal Equilibrium∗

2.10.4

Excitation Temperature

When the Boltzman distribution is achieved, the level population of 1 and 2 levels are written asN2g2=
exp(−E21 /kT ),N1
g1(2.155)where g1 and g2 represent the statistical weight of 1 and 2 levels, respectively. In the case that the
Boltzman distribution is not achieved, the temperature that satisfys eq(2.155) is called the excitation
temperature and is expressed as
E21/ logeTex = −
k

2.10.5





g1 N2
.
g2 N1

(2.156)

Critical Density

Consider a two-level atom (a hypothetical atom which has only two levels), in which the spontaneous
downward transitions and collisional excitations and deexcitations are in balance as(A21 + C21 n) N2 = C12 nN1 .

(2.157)

C12 g1 = C21 g2 exp [−(E2 − E1 )/kT ],(2.158)Using
the number fraction of upper level is written
g2
C12
exp [−(E2 − E1 )/kT ]
N2
g1
C21
,
=
=
ncr g2C12
A21
N1 + N2
+
1
+
exp
[−(E
−
E
)/kT
]
2
1
+
1+
n
g1
C21 n C21

(2.159)

where

A21
,
(2.160)C21
is called critical density. When nncr , the second term of the denominator is small and the level
population is given by the Boltzmann distribution. As long as n  ncr , the number of upper-level
population is much smaller than that expected for the Boltzmann distribution. Such a low-density
gas emits only weakly.
Since B-coeﬃcients, which has a meaning of the cross-section for the radiation, is proportional tothe electric dipole moment of the molecule, A-coeﬃcients are large for molecules with large electric
dipole moment (eq.[2.154]). In the case of rotational levels, A-coeﬃcients increase ∝ (J +1)3 and thus
the critical density increases for higher transition. In Table 2.1, the critical densities for rotationaltransitions of typical molecules are shown as well as A and C coeﬃcients. Comparing J = 1 − 0transitions of CO, CS, and HCO+ , CS and HCO+ trace higher-density gas than CO. And higher
transition J = 5 − 4 lines trace higher-density gas than lower transition J = 1 − 0 lines.
In the discussion above, we ignored the eﬀect of transition induced by absorption. The above
critical density is deﬁned for optical thin case. For a gas element with a ﬁnite optical depth, photons
are eﬀectively trapped in the gas element (photon trapping). If we use the probability for a photon
to escape from the gas element, βν , the critical density is reduced to
ncr ≡

ncr ∼

A21
βν .
C21

(2.161)CHAPTER 2. PHYSICAL BACKGROUND

62

Table 2.1: A and C coeﬃcients and critical densities for rotational transition of typical molecules(s−1 )A10
C10 (cm−3 s)
ncr 10 (cm−3 )A54 (s−1 )C54 (cm−3 s)
ncr 54 (cm−3 )CO
7.2 × 10−8
2.8 × 10−11
2.7 × 103
1.3 × 10−5
9.3 × 10−11
1.4 × 105CS
1.8 × 10−6
3.7 × 10−11
5 × 104
3.0 × 10−4
3.3 × 10−11
9 × 106

HCO+
4.3 × 10−5
2.6 × 10−10
1.6 × 105
7.2 × 10−3
4.2 × 10−10
1.7 × 107

The most abundant CO, 12 C16 O, is sometimes optically thick, while rare molecules 13 C16 O and12 C18 O, which have almost the same Einstein’s coeﬃcients as 12 C16 O, are optically thin. In such a
cloud, the critical density of 12 C16 O is smaller than that of 13 C16 O and 12 C18 O.

Chapter 3

Galactic Scale Star Formation
3.1

Schmidt Law

Schmidt (1959) speculated that the star formation rate is proportional to a power of the surface
density of the interstellar mediumΣ̇SF ∝ Σngas ,(3.1)
where the power n seems between 1 and 2 around the solar vicinity. If n = 2, the star formation
rate is thought to be determined by the collision rate of interstellar clouds. At that time Schmidt
showed us n  2. On the other hand, if the gas passing through the galactic arms forms stars, thestar formation rate seems proportional to the gas surface density and the arm-to-arm period. Thus
this predicts n = 1.

3.1.1

Global Star Formation

The star formation rate is estimated by the intensity of Hα emission (Kennicutt, Tamblyn, & Congdon
1994) as
L(Hα)
,
(3.2)
SFR(M yr−1 ) =1.26 × 1041 erg s−1
which is used for normal galaxies. While in the starburst galaxies which show much larger starformation rate than the normal galaxies, FIR luminosity seems a better indicator of star formation
rate
L(FIR)
L(FIR)=.
(3.3)
SFR(M yr−1 ) =43−1
2.2 × 10 erg s
5.8 × 109 L
Kennicutt (1998) summarized the relation between SFRs and the surface gas densities [Fig.3.1 (left)]for 61 normal spiral and 36 infrared-selected starburst galaxies. As seen in Figure 3.1, the star formation rate averaged over a galaxy (ΣSFR (M yr−1 kpc−2 )), which is called the global star formation
rate, is well correlated to the average gas surface density Σgas (M pc−2 ). He gave the power of the
global Schmidt law as n = 1.4 ± 0.15. That is,
ΣSFR  (1.5 ± 0.7) × 10−4Σgas
1M pc−21.4±0.15

M yr−1 k pc−2 .

(3.4)

The fact that the power is not far from 3/2 seems to be explained as follows: Star formation
rate should be proportional to the gas density (Σgas ) but it should also be inversely proportional to
the time scale of forming stars in respective gas clouds, which is essentially the free-fall time scale.
63

64

CHAPTER 3. GALACTIC SCALE STAR FORMATION

Figure 3.1: Taken from Figs.6 and 7 of Kennicutt (1998). Left: The x-axis means the total (HI+H2 )gas density and the y-axis does the global star formation rate. Right: The x-axis means the total(HI+H2 ) gas density divided by the orbital time-scale. The y-axis is the same.
Remember the fact that the free-fall time given in equation (2.26) is proportional to τﬀ ∝ 1/(Gρ)1/2 .
Therefore(3.5)
ρSFR ∝ ρgas × (Gρgas )1/2 ∝ ρ3/2
gas ,where ρgas and ρSFR are the volume densities of gas and star formation.
He found another correlation between a quantity of gas surface density divided by the orbital
period of galactic rotation and the star formation rate [Fig.3.1 (right)]. Although the actual slope isequal to 0.9 instead of 1, the correlation in Fig.3.1(right) is well expressed as
ΣSFR  0.017Σgas Ωgas = 0.21

Σgas,
τarm−to−arm(3.6)

where Ωgas represents the angular speed of galactic rotation. This means that 21 % of the gas mass isprocessed to form stars per orbit. These two correlations [eqs. (3.4) and (3.6)] implicitly ask another
1/2
relation of Ωgas ∝ Σgas .

3.1.2

Local Star Formation Rate

In Figure 3.2 (left), the correlation between star formation rate and gas density is plotted for speciﬁcgalaxies (NGC 4254 and NGC 2841). This shows us that Hα surface brightness (star formation

3.1. SCHMIDT LAW

65

Figure 3.2: Distributions of ΣSF and Σgas . (Left:) Relation between ΣSF and Σgas for an Sc galaxy
NGC 4254 and an Sb galaxy NGC 2841. (Right:) Relation between ΣSF and Σgas for various galaxies.
These are taken from Figs.7 and 8 of Kennicutt (1989).
rate) and the gas column density are well correlated each other. Figure 3.2 (left) also indicatesthat there seems a critical gas density below which star formation is not observed. The value of
this threshold column density is approximately  4M pc−2 for both galaxies in Figure 3.2 (left).
The same correlation is seen in other spiral galaxies [Fig.3.2(right)]. Fitting the correlation with apower-law, he obtained
1.3±0.3,(3.7)
ΣSFR ∝ Σgas
for the region active in star formation. Take notice that this power is very close to that of the
global Schmidt law [eq.(3.4)] The threshold surface gas density ranges from 1 M pc−2 to 10M pc−2(1020 − 1021 H cm−2 ). Therefore, theory of star formation must explain (1) the Schmidt law (clear
correlation between star formation rate and the gas surface density) above the threshold column
density and (2) the fact that there is no evidence for star formation in the gas deﬁcient region below
the threshold column density.
Star Formation Seen with Akari
Figure 3.3 represents the far-IR observation of M101 with Akari (Suzuki et al.2007). Using 4 bands of
Akari FIS and 3 bands of ISO, the SED of this galaxy is ﬁtted with two components of dust emissions+4
with temperatures Tw = 55+9
−25 K and Tc = 18−9 K. Spiral galaxies (not classiﬁed with either AGN or
starburst galaxy) have a far-IR temperature with 20K–50K as a whole from Revised Shapley-Ames
galaxy survey by IRAS De Jong et al 1984). As shown in Fig.3.4 (a), cold component well correlates

CHAPTER 3. GALACTIC SCALE STAR FORMATION

66

Figure 3.3: Akari far-IR observation of M101 (Suzuki et al.2007). (Left) SED of M101. Dotted and+4dashed lines represent best-ﬁt modiﬁed blackbody spectrum with Tw = 55+9
−25 K and Tc = 18−9 K.

(a)

(b)

Figure 3.4: Akari far-IR observation of M101 (Suzuki et al.2007). (a) Distribution of emissions fromcold dust with Tc (left) and warm dust with Tw (right). (b) The warm-to-cold intensity ratio isplotted.

3.2. DISTRIBUTIONS OF ATOMIC AND MOLECULAR HYDROGEN

67

with that of gas distribution, while warm one has a relation to star formation regions. Fig.3.4(b) plots
the warm-to-cold ratio of intensities, which seems to show the distribution of star formation eﬃciency≡ star formation rate/gas mass. Some peaks in this plot correspond to famous HII regions. Centerof the galaxy and the end of the galactic bar indicate peaks in this map. Since barred galaxies havelarger star formation rate than spiral galaxies without bars, the central peaks are regarded as a region
where gas is accumulated by the eﬀect of a bar potential. Equations (eq:Kennicutt-correlation-1) and(3.7) indicates that the star formation eﬃciency
0.3−0.4,
ΣSFR /Σgas ∝ Σgas

(3.8)

M101 galaxy do not seem to be explained by a simple law as the star formation eﬃciency is anincreasing function of gas density.

3.2

Distributions of Atomic and Molecular Hydrogen

We believe that stars form in molecular clouds, in which major component of Hydrogen is “molecule.”
MHI is larger than MH2 and HI gas distributes rather uniformly than H2 gas. Figure 3.5(a) shows
distributions integrated intensities of HI (red) and CO (blue) for Large Magellanic Cloud (Fukui etal. 2009). CO gas is found as GMCs with a mass M ∼ 105−6 M , while HI gas is distributed more
uniform than that of CO gas.
Figure 3.5(b) shows that line width of HI is much larger than that of CO. The comparison in thephase-space (position VS line-of-sight velocity) indicates us that the HI gas coexisting with H2 gas
was overestimated when the comparison is made for the column density σ. H2 gas associated with
GMCs may be a small fraction of HI gas.
Figure 3.6 displays the distributions of HI and CO gas. CO and HI data are divided into pixels
whose 3D size is 40pc× 40pc ×0.17km s−1 . The number of pixels where HI (> 7.2K) and CO(> 0.21K) are both detected are plotted by a histogram with dark color. While that of Tb (HI) > 7.2K
but Tb (CO) < 0.21K (light-shade histogram) shows a diﬀerent distribution. The ratio of pixels
associated with CO increases with Tb (HI) (Figure 3.6 bottom). This may indicates an accretion of
HI gas to form GMCs.
Another important ﬁnding of NANTEN survey is a classiﬁcation of GMCs from a stand point of
massive star formation. Kawamura et al (2009) indicate that “Type I shows no signature of massive
star formation, Type II is associated with relatively small HII region(s), and Type III with both HII
region(s) and young stellar cluster(s),” as shown in Figure 3.6(right). The lifetime of a GMC throughType I to Type III is estimated to be a few 107 yr. The average HI intensity for Types I, II, and III
are estimated as 34 ± 16K, 47 ± 17K, and 56 ± 19K. The HI intensity increases with the sequence ofstar formation.

68

CHAPTER 3. GALACTIC SCALE STAR FORMATION

Figure 3.5: (Left) Distributions integrated intensities of HI (red) and CO (blue) for LMC. The HI
integrated intensity is obtained from combination of aperture synthesis and single-dish observations.
(Kim et al. 03). (Right) Spectra of HI and CO for a direction are shown.

Figure 3.6: (Left, upper) Number of pixels are plotted against the brightness temperature of HI.
CO and HI data are divided into pixels whose 3D size is 40pc× 40pc×0.17km s−1 . (Left, bottom)
Number of pixels associated with CO emission increases with higher HI brightness temperature (Fukui
et al.2009). (Right) Classiﬁcation of GMCs. Type I GMCs are observed without any indicationof massive star formation; Type II’s are GMCs with HII regions. However, young clusters areembedded in the GMCs and not observed. Type III are associated with HII regions and young
clusters (Kawamura et al.2009)

3.3. GRAVITATIONAL INSTABILITY OF ROTATING THIN DISK

3.3

69

Gravitational Instability of Rotating Thin Disk

Here, we will derive the dispersion relation for the gravitational instability of a rotating thin disk. We
will see the spatial variation of Toomre’s Q parameter, which determines the stability of the rotating
disk, explains the nonlinearity of star formation rate, that is, there is a threshold density and no stars
are formed in the low density region.
Use the cylindrical coordinate (R, Z, φ) and the basic equations for thin disk in §2.6. In linear
analysis, we assume Σ(R, φ) = Σ0 (R) + δΣ(R, φ), u(R, φ) = 0 + δu(R, φ), v(R, φ) = v0 (R) + δv(R, φ),where u and v represent the radial and azimuthal components of the velocity. Linearized continuity
equation is∂δΣ1 ∂∂δΣ Σ0 ∂δv+(RΣ0 δu) + Ω
+
= 0,
(3.9)
∂t
R ∂R
∂φ
R ∂φ
where Ω = v0 /R.
Linearized equations of motion are


and



∂
∂
∂
+Ω
(δΦ + δh),δu − 2Ωδv = −
∂t
∂φ
∂R



(3.10)



∂
1 ∂
κ2
∂
+Ω
δu = −(δΦ + δh),
δv +
∂t
∂φ
2Ω
R ∂φ
where h is a speciﬁc enthalpy as dh = dp/Σ anddΩ2
κ = 4Ω + R
dR(3.11)

1/2

2(3.12)

is the epicyclic frequency.
We assume any solution of equations (3.9), (3.10) and(3.11) can be written as a sum of terms of
the formδu = ua exp[i(mφ − ωt)],(3.13)

δv = va exp[i(mφ − ωt)],(3.14)

δΣ = Σa exp[i(mφ − ωt)],(3.15)

δh = ha exp[i(mφ − ωt)],(3.16)

δΦ = Φa exp[i(mφ − ωt)].

(3.17)

Using the equation of state of p = KΣγ ,ha = c2s Σa /Σ0 .

(3.18)

Using equations (3.13)-(3.17), equations (3.9), (3.10), and (3.11) are rewritten asi(mΩ − ω)Σa +


1 ∂
Σ0 va
(RΣ0 ua ) + im
= 0,
R ∂R
R

ua [κ2 − (mΩ − ω)2 ] = −i (mΩ − ω)and



(3.19)d(Φa + ha )(Φa + ha ) + 2mΩ
,
dR
R


(Φa + ha )
κ2 d(Φa + ha ) + m(mΩ − ω),va [κ − (mΩ − ω) ] =
2Ω dR
R
2

2

(3.20)(3.21)

CHAPTER 3. GALACTIC SCALE STAR FORMATION

70

Figure 3.7: Tightly wound (left) vs loosely wound (right) spirals.

3.3.1

Tightly Wound Spirals

We assume the wave driven by the self-gravity has a form of tightly-wound spiral [Fig.3.7(left)].
When we move radially, the density δΣ varies rapidly. While, it changes its amplitude slowly in theazimuthal direction. In a mathematical expression, if we write the density perturbation δΣ asδΣ = A(R, t) exp[imφ + i f (R, t)],(3.22)

where the amplitude of spiral A(R, t) is a slowly varing function of R, a tightly wound spiral means
the shape function varies fast (the radial wavenumber k  df /dR is large enough). We consider thegravitational force from the vicinity of (R0 , φ0 ), since the δΣ oscillates and cancels even if we integrate
over large region. Thus,(3.23)
δΣ(R, φ, t)  Σa exp[ik(R0 , t)(R − R0 )],whereΣa = A(R0 , t) exp[imφ0 + f (R0 , t)].

(3.24)

Notice that the density perturbation [eq.(3.23)] is similar to that studied in §2.6. The potential should
be expressed in a similar form to equation (2.64) asδΦ  −

2πGΣaexp[ik(R0 , t)(R − R0 )],|k|(3.25)

2πGΣa
.
|k|

(3.26)

which simply meansΦa = −

If we set R = R0 , we obtain our ﬁnal result for the potential due to the surface density perturbation
δΦ(R, φ, t)  −2πG
A(R, t) exp[imφ + f (R, t)].
|k|

(3.27)

Diﬀerentiating this equation with R and ignoring the term dA(R, t)/dR compared to that of df (R, t)/dR =k, we obtainsign(k) dδΦ(R, φ, t),(3.28)
δΣ(R, φ, t) = i
2πG
dR

3.3. GRAVITATIONAL INSTABILITY OF ROTATING THIN DISK

71

Table 3.1: Epicyclic frequency vs rotation law.
Rotation
Rigid-body rotation Ω =const.
Flat rotation vφ =const.
Kepler rotation vφ ∝ r −1/2

κ
2Ω
√
2Ω
Ω

Neglecting the terms ∝ 1/R compared to the terms containing ∂/∂R, equations (3.19), (3.20),
and (3.21) are rewritten as(3.29)i(mΩ − ω)Σa + ikΣ0 ua = 0,ua [κ2 − (mΩ − ω)2 ] = (mΩ − ω)k(Φa + ha ),(3.30)andκ2k(Φa + ha ),(3.31)
2Ω
Using these equations [(3.29), (3.30), and (3.31)], Φa = −2πGΣa /|k|, and ha = c2s Σa /Σ0 , we obtain
the dispersion relation for the self-gravitating instability of the rotating gaseous thin diskva [κ2 − (mΩ − ω)2 ] = i(mΩ − ω)2 = k2 c2s − 2πGΣ0 |k| + κ2 .

(3.32)

Generally speaking, the epicyclic frequency depends on the rotation law but is in the range of Ω<∼κ<∼2Ω(see Table 3.1 for κ for typical rotation laws). It is shown that the system is stabilized due to theepicyclic frequency compared with a nonrotating thin disk [eq.(2.65)].

3.3.2

Toomre’s Q Value

Consider the case of m = 0 axisymmetric perturbations. Equation (3.32) becomesω 2 = k2 c2s − 2πGΣ0 |k| + κ2 = c2s k −

πGΣ0
c2s

2



+ κ2 −

πGΣ0
cs

2

.

(3.33)

If ω 2 > 0 the system is stable against the axisymmetric perturbation, while if ω 2 < 0 the system isunstable. Deﬁning
κcs,(3.34)Q=
πGΣ0
if Q > 1, ω 2 > 0 for all wavenumbers k. On the other hand, if Q < 1, ω 2 becomes negative forsome wavenumbers k1 < k < k2 . Therefore, the Toomre’s Q number gives us a criterion whether thesystem is unstable or not for the axisymmetric perturbation. [Recommendation for a reference bookof this section: Binney & Tremaine (1988).]
The condition is expressed as
κcs(Q < 1).
(3.35)
Σ0 > Σcr =πG
Kennicutt plotted Σ0 /Σcr against the normalized radius as R/RHII for various galaxies, where RHII
represents the maximum distance of HII regions from the center (Fig.3.8). Since Σ0 /Σcr = Q, Figure
3.8 shows that HII regions are observed mainly in the region with Q < 1 but those are seldom seenin the outer low-density Q > 1 region. This seems the gravitational instability plays an importantrole.

CHAPTER 3. GALACTIC SCALE STAR FORMATION

72

Figure 3.8: Σ0 /Σcr vs R/RHII . RHII represents the maximum distance of HII regions from the center.
The sound speed is assumed constant cs = 6km s−1 . Taken from Fig.11 of Kennicutt (1989).
The above discussion is for the gaseous disk. The Toomre’s Q value is also deﬁned for stellar
system asσR κ,(3.36)
Q=
3.36GΣ0where σR represents the radial velocity dispersion.
For non-axisymmetric waves, even if 1<∼ Q<∼2 the instability grows. To explain this, the swing
ampliﬁcation mechanism is proposed (Toomre 1981). If there is a leading spiral perturbation in thedisk with 1<∼ Q<∼2, the wave unwinds and ﬁnally becomes a trailing spiral pattern. At the same time,
the amplitude of the wave (perturbations) is ampliﬁed (see Fig.3.9).

3.4

Spiral Structure

Figure 3.10 shows the B- (left) and I-band images of M51. B-band light originates from the massive
early type stars. Although the image taken in B-band shows a number of spiral arms, that of I-band
shows clearly two arms. The I-band light seems to come from mainly less-massive long-lived stars,
while the B-band light is essentially coming from the massive short-lived stars which are formed in
the spiral arm. On the contrary, the less-massive stars are not necessarily born in the spiral arm.
This suggests that there are two kinds of spiral patterns: one is made by stars (mainly less-massive)
and the other is the gaseous spiral arm where massive stars are born and contribute to the B-band
image.
In this section, ﬁrst, we will brieﬂy describe the density wave theory which explains the former
spiral pattern in the stellar component. You will ﬁnd the amplitude of the spiral pattern in stellar
component is not so large. However, the response of gaseous components (HI and H2 gas) to thespiral density wave potential with small amplitude is much more nonlinear than that of the stellar
component and a high-contrast spiral pattern appears in the gaseous component .

3.5. DENSITY WAVE THEORY

73

Figure 3.9: Numerical simulation of the swing ampliﬁcation mechanism. The number attached each panel
shows the time sequence. This is obtained by the time-dependent linear analysis. First, perturbation withleading spiral pattern is added to the Mestel disk with Q = 1.5. The leading spiral gradually unwinds andbecome a trailing spiral. Loosely wound spiral pattern winds gradually and the last panel shows a tightly
wound leading spiral pattern. The ﬁnal amplitude is ∼ 100 times larger than that of the initial state.

3.5

Density Wave Theory

We have derived the dispersion relation of the gravitational instability in the rotating thin disk as(mΩ − ω)2 = k2 c2s − 2πGΣ0 |k| + κ2 ,(3.37)

where m represents the number of spiral arms. Although the stability of the stellar system is a little
diﬀerent, we assume this is valid for the stellar system after cs is replaced to the velocity dispersion.
Since




πGΣ0 2πGΣ0 2
2
2 2
2
2
2
+κ −
,
(3.38)(mΩ − ω) = k cs − 2πGΣ0 |k| + κ = cs k −c2scs
we obtain|k| = kT2
1±
1 − Q2 (1 − ν 2 ) ,Q2(3.39)

where ΩP = ω/m is a pattern speed, ν = m(ΩP − Ω)/κ is the normalized frequency, kT = κ2 /2πGΣ0
is the Toomre’s critical wavenumber for a cold (cs = 0) system. ν = ±1, which leads to |k| = 0,

CHAPTER 3. GALACTIC SCALE STAR FORMATION

74

Figure 3.10: B- (left), and I-band (right) images of M51. Taken from Fig. 3 of Elmegreen et al.
(1989).

k/kT

Q=1

4

Q=1.2

3.5

Q=1.5

3
2.5

Q=2

2
1.5
1
0.5
-2

-1

1

2

nu

Figure 3.11: (Left:) A plot of the dispersion relation eq.(2.65). x- and y-axes are ν and k/kT .
Respective lines are for Q = 1 (straight lines), Q = 1.2, Q = 1.5, and Q = 2. The trailing partk > 0 is only plotted. Points ν = −1, ν = 0, and ν = 1 correspond respectively to Inner Lindbrad
Resonance (ILR), Corotation Resonance (CR), and Outer Lindbrad Resonance (OLR). The relation
is symmetric against the x-axis and the curve of k < 0 represents the leading wave. (Right:) Leadingvs trailing spiral.

3.5. DENSITY WAVE THEORY

75

represents the Lindbrad resonance and is rewriten as
ΩP =

ω
κ
=Ω± .
m
2

(3.40)

Assuming m = 2, the resonance when ΩP = Ω + κ/2 is called outer Lindbrad resonance while that ofΩP = Ω − κ/2 is called inner Lindbrad resonance. ν = 0 means the co-rotation resonance ΩP = Ω.
Plotting the wavenumber k against the normalized frequency ν of equation (3.39) as Figure
3.11(left), it is shown that, in the case of Q = 1, the wavenumber exists for all ν. Since ν = −1, 0,
and +1 correspond to the points of ILR, CR, and OLR and these three resonance points appear inaccordance with the radial distance, the x-axis of Figure 3.11(left) seems to correspond to the radial
distance from the center. In the case of Q > 1, it is shown that a forbidden region appears aroundthe co-rotation resonance point. Waves cannot propagate into the region. Figure 3.11(left) showsthat the k/kT has two possible wavenumbers in the permitted region. The waves with larger k andsmaller k are called short waves and long waves, respectively.
Consider a wave expressed by Σ ∝ exp[imφ + ikr]. If k < 0, moving from a point (R0 , φ0 ) in the
direction Δφ > 0 and Δr > 0 the phase diﬀerence between the two points [m(φ0 + Δφ) + k(R0 +Δr)] − [mφ0 + kR0 ] can be equal to zero. That is, in the case of k < 0 the wave is leading. On
the other hand, if k > 0, moving in the direction Δφ < 0 and Δr > 0 the phase will be unchanged.
In this case, the wave pattern is trailing. Since the dispersion relation is symmetric for k > 0 andk < 0, there are two waves, trailing waves and leading waves. Therefore there are four waves: a short
trailing wave, a long trailing wave, a short leading wave, and a long leading wave.

3.5.1

Group Velocity

The wave transfers the energy with the group velocity. Whether the group velocity is positive or
negative is quite important considering the energy transfer. Using the dispersion relation, equation(3.38), The group velocity|k|c2s − πGΣ0dω
= sign(k).
(3.41)
vg (R) =dk
ω − mΩ
For a region R > RCR , ω − mΩ > 0 or ν > 0. On the other hand, for a region R < RCR , ω − mΩ < 0or ν < 0. Consider ﬁrst the trailing wave. In the region R > RCR , long-waves propagate inwardly
to the CR, since |k|c2s < πGΣ0 for long-waves and ν > 0. Short-waves propagate outwardly from theCR, since |k|c2s > πGΣ0 for short-waves and ν > 0. In the region R < RCR , long-waves propagate
outwardly to the CR since |k|c2s < πGΣ0 and ν < 0. Short-waves propagate inwardly from the CR,
since |k|c2s > πGΣ0 and ν < 0. As a result, it is concluded that the long-wave propagates toward the
CR and the short-wave does away from the CR. As for the leading wave, the short-wave propagates
toward the CR and the long-wave does away from the CR.
Assuming that the wave packet is made near the Lindblad resonance points, (1) the long-trailingwaves propagate toward the co-rotation resonance points; (2) they are reﬂected by the Q-barrier;(3) they change their character to short-waves and propagate away from the co-rotation resonancepoints; and ﬁnally (4) the waves are absorbed at the center or propagate away to the inﬁnity. (1) the
short-leading waves propagate toward the co-rotation resonance points; (2) they are reﬂected by the
Q-barrier; (3) they change their character to long-waves and propagate away from the co-rotation
resonance points; and ﬁnally (4) they reach the Lindblad resonance points and the energy may be
converted to the long-trailing waves there.
The wave obtains its energy at the resonance points. The density wave transfer the energy to
the co-rotation points. Therefore, the density wave theory predicts the galactic stellar disk has spiral

CHAPTER 3. GALACTIC SCALE STAR FORMATION

76

Figure 3.12: Spiral coordinate. Along a curve of constant η, r ∝ exp(θ tan i).
density pattern between the inner Lindblad resonance points and the outer Lindblad resonance point
if Q>
∼1.

3.6

Galactic ShockIn the preceding section, we have seen that a spiral pattern of density perturbation is made by thegravitational instability in the stellar system with Q>∼1. This forms a grand design spiral observedin the I-band images which represent the mass distribution which consists essentially in less-massive
stars. The amplitude of the pattern is smaller than that observed in B-band images. In this section
we will see how the distribution early type stars is explained. [Recommendation for a reference bookof this section: Spitzer (1978).]
Gas ﬂowing through a spiral gravitational potential, even if its amplitude is relatively small, acts
rather nonlinearly. Consider the gravitational potential in the sum of an axisymmetric term, Φ0 (R)and a spiral term Φ1 (R). The spiral gravitational ﬁeld is rotating with a pattern speed ΩP . We use
a reference frame rotating with ΩP . u = (u, v), where u and v are radial (R) and azimuthal (φ)
component of the ﬂow velocity. Consider the reference system rotating with the angular rotation
speed ΩP . We assume the variables with suﬃx 0 represent those without the spiral potential Φ1and variables with suﬃx 1 are used to express the diﬀerence between before and after Φ1 is added.
u = u0 + u1 = 0 + u1 , v = v0 + v1 , andv0 = R(Ω − ΩP ) = R

1 ∂Φ0
R ∂R

1/2



− ΩP ,(3.42)

where Ω means the circular rotation speed. The equation of motion in the r-direction is∂u v ∂u v12
c2 ∂Σ ∂Φ1
∂u
+u
+
−
=− s
−
+ 2Ωv1 ,
∂t
∂r
r ∂θ
r
Σ ∂r
∂r
and that in the φ direction is

∂v1 v1
∂v
+u
+
∂t
∂r
r



+

c2 ∂Σ 1 ∂Φ1
κ2
v ∂v1
=− s
−
−
u,
r ∂φ
Σr ∂φ
r ∂φ
2Ω

(3.43)(3.44)

3.6. GALACTIC SHOCK

77

and the continuity equation is∂ΣRu ∂Σv
+
= 0.
∂R
∂φ

(3.45)As seen in Figure 3.12, we introduce the spiral coordinate (η, ξ) in which ξ-axis is parallel to thespiral pattern which has a pitch angle i and η-axis is perpendicular to the ξ-axis.
rdη = dr cos i + rdφ sin i,(3.46)rdξ = −dr sin i + rdφ cos i.

(3.47)vη = u cos i + v sin i,(3.48)vξ = −u sin i + v cos i(3.49)

Assuming that i  1 (tightly wound spiral), equations (3.45) (3.43) and (3.44) becomec2 ∂Σ ∂Φ1
∂vη
∂vη
+ vη
=− s
−
+ 2Ω(vξ − vξ0 ),
∂t
∂η
Σ ∂η
∂η
∂vξ
∂vξ
κ2
+ vη
= − (vη − vη0 ),∂t∂η2Ω∂Σ ∂Σvη
+
= 0.
∂t
∂η

(3.50)(3.51)(3.52)

Similar to §2.8 we look for a steady state solution. Equation (3.52) becomes
Σ dvη
dΣ
=−
.
dη
vη dη

(3.53)

Using this, equation (3.50) reduces tovη2 − c2s

 dv

η

dη

= 2Ω(vξ − vξ0 ) −

dΦ1
.
dη

(3.54)

Equation (3.51) becomesdvξκ2= − (vη − vη0 ).
dη2Ω
In these equations we used following quantities:vη

(3.55)vη = vη0 + vη1 ,(3.56)vξ = vξ0 + vξ1 ,(3.57)vη0 = v0 sin i,(3.58)vξ0 = v0 cos i,(3.59)

v0 is given in equation (3.42).
Equations (3.53), (3.54), and (3.55) are solved under the periodic boundary condition: X(η =right end) = X(η = left end). Since equation (3.54) is similar to the equations in §2.8, you may think
the solution seems like Figure 2.6. However, it contains a term which expresses the eﬀect of Coliois
force 2Ω(vξ − vξ0 ), the ﬂow becomes much complicated.
Taking care of the point that a shock front exists for a range of parameters, the solution of
the above equations are shown in Figure 3.13. Numerical hydrodynamical calculations which solve
equations (3.50), (3.51), and (3.52) was done and steady state solutions are obtained (Woodward

78

CHAPTER 3. GALACTIC SCALE STAR FORMATION

1975). It is shown that F ≤ 0.7%, the velocity (vη : velocity perpendicular to the wave) does notshow any discontinuity. In contrast, for F > 2%, a shock wave appears. Steady state solution
is obtained from the ordinary diﬀerential equations (3.53) (3.54), and (3.55) by Shu, Milione, &Roberts (1973). Inside the CR, v0 > 0 (Gas has a faster rotation speed than the spiral pattern). As
long as vη > cs there is no shock. Increasing the amplitude of the spiral force F , an amplitude of thevariation in vη increases and ﬁnally vη becomes subsonic partially. When the ﬂow changes its naturefrom supersonic to subsonic, it is accompanied with a shock. (An inverse process, that is, changing
from subsonic to supersonic is not accompanied with a shock.)In the outer galaxy (still inside CR) since Ω − ΩP decreases with the the distance from the galactic
center, vη0 = R(Ω − ΩP ) sin i decreases. In this region, vη0 < cis and the ﬂow is subsonic if there isno spiral gravitational force, F = 0. In such a circumstance, increasing F ampliﬁes the variation invη and vη ﬁnally reaches the sound speed. Transonic ﬂow shows again a shock.
Summary of this section is:
1. There exists a spiral density pattern of the stellar component if the Toomre’s Q parameter is
1<
∼ Q<
∼2.
2. This is driven by the self-gravity of the rotating thin disk.
3. If the amplitude of the non-axisymmetric force is as large as ∼ 1% of the axisymmetric one,interstellar gas whose sound speed is as large as cis ∼ 8km s−1 forms the galactic shock. Thisis observed when the ﬂow is transonic.
4. The amplitude of the gas density ﬂuctuation is much larger than that in the stellar density.
5. If stars are formed preferentially in the postshock region of the spiral arm, we expect clear
spiral arms made by early-type stars, which are massive and short-lived, as seen in the B-band
images of spiral galaxies.

3.6. GALACTIC SHOCK

79

Figure 3.13: Galactic shocks. The velocity perpendicular to the spiral wave front is plotted against thephase of the spiral wave. Gas ﬂows from the left to the right. The spiral potential takes its minimumat θ = 90 deg. Parameters are taken as cs = 8.6km s−1 , i = 6. deg 7, R = 10kpc, RΩ = 250km s−1 ,Rκ = 313km s−1 , RΩP = 135km s−1 . The amplitude of spiral gravity force is taken F =0.4%, 0.7%,
2%, and 5% of the axisymmetric force ∂Φ0 /∂R. Taken from Fig.13.3 of Spizer (1978) [originally
Woodward (1975)].

80

CHAPTER 3. GALACTIC SCALE STAR FORMATION

Figure 3.14: Galactic shock. (Left:) vη vs vξ plot. Increasing F from 0.5% to 0.9%, the variationin vη increases. And at F = 0.97%, the minimum speed reaches the sound speed 8km s−1 . Further
increasing F , a subsonic region appears (vη < cis for F = 1.05% and 5%). When ﬂow changes its
nature from supersonic to subsonic, a shock appears. R = 10kpc, cs = 8km s−1 , and i = 6. deg 7 aretaken. Other parameters are the same as Fig.3.13. (Right:) Comparison of two density distributions:
there are no shocks only cusps for F < 0.97%, while shochs and high compressions in the postshock
region appear for models with F = 5%. Taken from Shu, Milione, & Roberts (1973).

Chapter 4

Local Star Formation Process
4.1

Hydrostatic Balance

Consider a hydrostatic balance of isothermal cloud. By the gas density, ρ, the isothermal soundspeed, cis , and the gravitational potential, Φ, the force balance is written as
c2 dρ dΦ
−
= 0,
− is
ρ dr
dr

(4.1)

and the gravity is calculated from a density distribution as
−

GMr
4πG
dΦ
=− 2 =− 2
dr
r
r

 r

ρr 2 dr,(4.2)

0

for a spherical symmetric cloud, where Mr represents the mass contained inside the radius r. The
expression for a cylindrical cloud is
2Gλr
4πG
dΦ
=−
=−
−
dr
r
r

 r

ρrdr,

(4.3)

0

where λr represents the mass per unit length within a cylinder of radius being r.
For the spherical symmetric case, the equation becomes the Lane-Emden equation with the polytropic index of ∞ (see Appendix D.1). This has no analytic solutions. However, the numericalintegration gives us a solution shown in Figure 4.1 (left). Only in a limiting case with the inﬁnite
central density, the solution is expressed as
ρ(r) =

c2is −2
r .
2πG

(4.4)

Increasing the central density, the solution reaches the above Singular Isothermal Sphere (SIS) solution.
On the other hand, a cylindrical cloud has an analytic solution (Ostriker 1964) as
ρ(r) = ρc 1 +

r2
8H 2

−2

,

(4.5)

where
H 2 = c2is /4πGρc .
81

(4.6)

CHAPTER 4. LOCAL STAR FORMATION PROCESS

82
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

1000

1e-05
0.01

0.1

1

10

100

1000

Figure 4.1: Radial density distribution. A spherical cloud (left) and a cylindrical cloud (right). In theright panel, solutions for polytropic gases with Γ = 1.1 (relatively compact) and Γ = 0.9 (relatively
extended) are plotted as well as the isothermal one.
Far from the cloud symmetric axis, the distribution of equation (4.5) gives
ρ(r) ∝ r −4 ,(4.7)

ρ(r) ∝ r −2

(4.8)

while the spherical symmetric cloud has

distribution.
Problem 1
Show that the SIS is a solution of equations (4.1) and (4.2).
Problem 2
Show that the density distribution of [equation (4.5)] is a solution of equations (4.1) and (4.3).

4.1.1

Bonnor-Ebert Mass

In the preceding section [Fig.4.1 (left)], we have seen the radial density distribution of a hydrostatic
conﬁguration of an isothermal gas. Consider a circumstance that such kind of cloud is immersed ina low-density medium with a pressure p0 . To establish a pressure equilibrium, the pressure at thesurface c2is ρ(R) must equal to p0 . This means that the density at the surface is constant ρ(R) = p0 /c2is .
Figure 4.2 (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 4.2 (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
.
(4.9)
Mmax  1.14
1/2
G3/2 p0

4.1. HYDROSTATIC BALANCE

83

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 .
Another important result from Figure 4.2 (right) is the stability of an isothermal cloud. Even fora cloud with Mcl < Mmax , any clouds on the part of ∂Mcl /∂ρc < 0 are unstable, whose clouds are
distributed on the branch with ρc > ρcr . This is understood as follows: For a hydrostatic cloud themass should be expressed with the external pressure and the central density [Fig.4.2 (right)] as
Mcl = Mcl (pext , ρc ).

(4.10)In this case, a relation between the partial derivatives such as


∂Mcl
∂pext





·
ρc

∂pext
∂ρc





·
Mcl

∂ρc
∂Mcl



= −1,

(4.11)

pext

is satisﬁed, unless each term is equal to zero. Figure 4.1 (left) shows that the cloud mass Mcl is adecreasing function of the external pressure pext = ρs c2is , if the central density is ﬁxed. Since this
means


∂Mcl
< 0,
(4.12)∂pext ρcequation (4.11) gives us



∂pext
∂ρc





·
Mcl

∂ρc
∂Mcl



> 0.

(4.13)

pext

For a cloud with ρc < ρcr = 16pext /c2is the mass is an increasing function of the central density as


∂Mcl
∂ρc



> 0.

(4.14)> 0,(4.15)

pext

Thus, equation (4.13) leads to the relation


∂ρc
∂pext


Mcl

for ρ < ρcr . This means that gas cloud contracts (the central density and pressure increase), whenthe external pressure increases. This is an ordinary reaction of a stable gas.
2
On the other hand, the cloud on the part of (∂Mcl /∂ρc )pext < 0 (for 16<∼ρc /(pext /cis )<∼2000)
behaves


∂ρc
< 0,
(4.16)∂pext Mcland this represents that an extra external pressure decreases the central density and thus the pressure.
Pressure encourages expansion. This reaction is unstable.

4.1.2

Equilibria of Cylindrical Cloud

In Figure 4.1 (right) we plotted the structure for a polytropic cloud. Inner structure is not dependentof Γ, it is clear the slope of the outer envelope is dependent on Γ.

CHAPTER 4. LOCAL STAR FORMATION PROCESS

84
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

1

10

1

x

10

100

1000

10000

rho

Figure 4.2: (Left) radial density distribution. Each√curve has diﬀerent ρc . The x-axis denotes theradial distance normalized by a scale-length as c2is / 4πGpext . It is shown that the radius increases
with decreasing ρc in this range (ρc = (10, 102 , 103 ) × pext /c2is ). (Right) The relation between mass
(Mcl ) and the central density (ρc ) is plotted, under the condition of constant external pressure. The
x-axis represents the central √density normalized by ρs ≡ pext /c2is . The y-axis represents the cloudmass normalized by 4πρs (cis / 4πGρs )3 .
1. In the case of the spherically symmetric, consider a polytrope (p ∝ ρΓ ) with Γ < 6/5 (at least
the envelope of Γ = 6/5 cloud extends to ∞.), in which gas extends to ∞. if ρ ∝ r −p , themass inside of r is proportional to Mr ∝ r 3−p . Thus, the gravity per unit volume at a radiusr, GM ρ/r 2 , is proportional to GM ρ/r 2 ∝ r 1−2p . On the other hand the pressure force is|(∂p/∂r)| = (∂p/∂ρ)|(∂ρ/∂r)| ∝ (r −p )Γ−1 r −p−1 ∝ r −pΓ−1 . These two powers become the same,
only if p = 2/(2 − Γ).
2. In the case of cylindrical cloud with Γ ≤ 1, the mass per unit length λ ∝ r 2−p . The gravity atr, Gλρ/r ∝ r 1−2p . Note that the power is the same as the spherical case. Since the power ofthe pressure force should be the same as the spherical case, the resultant p should be the samep = 2/(2 − Γ).
The case of Γ = 0.9, an envelope extending to a large radius indicates the power-law distribution
much shallower than that of the isothermal Γ = 1 one.

4.2

Virial Analysis

Hydrodynamic equation of motion using the Lagrangean time derivative [eq.(B.3)] is


ρ

dv
dt



= −∇p − ρ∇Φ.

(4.17)

Multiplying the position vector r and integrate over a volume of a cloud, we obtain the Virial relationas
1 d2 I= 2(T − T0 ) + W,
(4.18)
2 dt2
where

(4.19)
I = ρr 2 dV,

4.2. VIRIAL ANALYSIS

85

is an inertia of the cloud,

 

T =



3
1
3
pth + ρv 2 dV = p̄Vcl ,2
2
2(4.20)

is a term corresponding to the thermal pressure plus turbulent pressure,3
Pth r · ndS = P0 Vcl
2
S

T0 =

(4.21)

comes from a surface pressure, andW =−ρr · ∇ΦdV = −3 GM 25 R(4.22)

is a gravitational energy. To derive the last expression in each equation, we have assumed the cloudis spherical and uniform. Here we use a standard notation as the radius R, the volume Vcl = 4πR3 /3,
the average pressure P̄ , and the mass M .
To obtain a condition of mechanical equilibrium, we assume d2 I/dt2 = 0. Equation (4.18) becomes
4π p̄R3 − 4πp0 R3 −3 GM 2= 0.
5 R(4.23)

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

3M
.
4πR3

(4.24)

Using equation (4.24) to eliminate p̄ from equation (4.18), the external pressure is related to the massand the radius as
3c2 M
3GM 2
.
(4.25)
p0 = is 3 −
4πR
20πR4
Keeping M constant and increasing R from zero, p0 increases ﬁrst, but it takes a maximum, p0,max =3.15c8is /(G3 M 2 ), and ﬁnally declines. This indicates that the surface pressure must be smaller thanp0,max for a cloud to be in the equilibrium. In other words, keeping p0 and changing R, it is shownthat M has a maximum value to have a solution. The maximum mass is equal toMmax = 1.77

c4is
1/2

G3/2 p0.

(4.26)

The cloud massive than Mmax cannot be supported against the self-gravity. This corresponds to the
Bonnor-Ebert mass [eq.(4.9)], although the numerical factors are slightly diﬀerent.

4.2.1

Magnatohydrostatic Clouds

Consider here the eﬀect of the magnetic ﬁeld. In the magnetized medium, the Lorentz force
F =

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

works in the ionized medium. The ﬁrst term of equation (4.27), which is called the magnetic
pressure, has an eﬀect to support the cloud against the self-gravity.

CHAPTER 4. LOCAL STAR FORMATION PROCESS

86

The virial analysis is also applicable to the magnetohydrostatic clouds. The terms related to themagnetic ﬁelds are


M





B2
B2
dV + (r · B)B · ndS −r · ndS
=
8π
S
S 8π

B 2 − B02
Φ2B
1
Φ2

dV  2
− B ,
8π
6π
R
R0

(4.28)

where ΦB represents a magnetic ﬂux and it is assumed to be conserved if we change the radius, R,
that is ΦB = πB0 R02 = πBR2 . Equation (4.23) becomes
4π p̄R3 − 4πp0 R3 −1 Φ2
3 GM 2+ 2 B = 0,
5 R
6π R

(4.29)

Φ2where we ignored the term 6π1 2 RB0 . The last two terms are rewritten as3G
M 2 − MΦ2 ,5R(4.30)

where MΦ is deﬁned as 3GMΦ2 /5 = Φ2B /6π 2 .
This shows the eﬀects of the magnetic ﬁelds:
1. B-ﬁelds eﬀectively 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 a cloud contracts (decreasing its radius from R0 to R), the ratio of the gravitational 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 to predominate over the magnetic term.
If M < MΦ , a sum of last two terms in equation (4.29) is positive. Since the second term of rhs ofequation (4.25) is positive, there is one R which satisﬁes equation (4.29) 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 ﬁelds work to support the cloud or not. Morerealistic calculation [Mouschovias (1976a,1976b), Tomisaka et al (1988)] gives us a criterion
G1/2

1
G1/2 σdm= 0.17 
,
=
dΦB
B
2π

(4.31)

where σ and B means the column density and the magnetic ﬂux density. A cloud with a massM>ΦB
2πG1/2(4.32)

is sometimes called magnetically supercritical, while that withM<

ΦB
2πG1/2(4.33)

is magnetically subcritical.
More precisely speaking, the criterion showed in equations (4.32) and (4.33) should be applied fora cloud which has a much larger mass than the Bonnor-Ebert mass. That is, even without magnetic
ﬁelds, the cloud less-massive than the Bonnor-Ebert mass has a hydrostatic conﬁguration shown in4.2. VIRIAL ANALYSIS

87

Figure 4.3: Cloud mass is plotted against the central density. Respective curves correspond to clouds
with diﬀerent magnetic ﬂux. Plasma beta of the external medium (β0 ) speciﬁes the magnetic ﬂux.

Figure 4.2 (left). The cloud with central density of ρc = 10 has a stable density distribution. Magnetohydrostatic clouds with diﬀerent magnetic ﬂuxes are calculated by 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 ﬂux (Fig. 4.3). The maximum allowable mass (critial mass) supported
by some magnetic ﬂux increases with the magnetic ﬂux. To ﬁt the numerical results, Tomisaka et al(1988) obtained an expression for the critical mass when the cloud has a mass-to-ﬂux ratio dm/dΦB ,the isothermal sound speed cis , and the external pressure p0 as


2
1/2π
Mcr = 1.3 1 −G1/2 dm/dΦB |r=0

−3/2

c4s
1/2

p0 G3/2.

(4.34)

This shows that the critical mass is a decreasing function of the mass-to-ﬂux ratio or increasing
function of the magnetic ﬂux. And the critical mass becomes much larger than the Bonnor-Ebert
1/2
mass  c4s /(p0 G3/2 ) only when the mass-to-ﬂux 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.

CHAPTER 4. LOCAL STAR FORMATION PROCESS

88

4.3

Subcritical Cloud vs Supercritical Cloud

We have seen there is a critical mass above which the cloud has no (magneto)hydrostatic conﬁguration
but below which the cloud has at least an equilibrium state. This gives us an idea that there are
two kind of clouds/cloud cores: that with a mass larger than the critical mass which has to collapse
dynamically and that with a mass smaller than the critical mass which is in an equilibrium state.
10−3
In the density range of 104 cm−3 <∼n<∼10 cm , the interstellar gas is essentially isothermal. In thisregion a major cooling agent is dusts; that is, the dust is heated by the collision of molecules. The
excess energy liberated at the collision increases the dust temperature. Finally the thermal emission
from the dust cools down the dust again. By this process, the thermal energy of the gas is reduced.
Therefore, we consider the cloud/cloud core is isothermal and study the collapse of the isothermal
cloud.
Since the supercritical cloud has no hydrostatic conﬁguration, it must evolve in a dynamical way.
On the other hand, since the subcritical cloud is in a static state, it evolves in much longer time-scaleof the free-fall time. Such cloud evolves by the eﬀect called the ambipolar diﬀusion.

4.3.1

Observation of Magnetic Field Strength

To determine whether a cloud is magnetically subcritical or supercritical, we have to measure thestrength of the magnetic ﬁeld.
Zeeman Splitting Method
If the atoms and molecules are permeated by magnetic ﬁeld B, in the normal Zeeman splitting, thetransition splits into three lines
ν1 = ν0 − β|B|,(4.35)

ν2 = ν0 ,(4.36)

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

where β = eh̄/2me = 1.3996Hz μG−1 is the Bohr magneton constant. This comes from the factthat the upper and lower levels are splited 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 to measure the magnetic ﬁeld strength. Figure 4.4 represents the correlation between themagnetic ﬁeld strength and gas density compiled by Fiebig & Guesten (1989), which indicates that
B is approximately proportional to ρ1/2
Chandrasekhar-Fermi Method
Chandrasekhar & Fermi (1953) obtained the strength of interstellar magnetic ﬁeld. Applying the same
method to the molecular cloud, Crutcher et al. (2004) obtained the magnetic ﬁeld strength in theplane 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 velocity of transverse Alfvén wave, VA = B/(4πρ)1/2 , thewave is described as(4.38)
y = A cos [k(x − VA t)] ,where y is a displacement in the plane of the sky and perpendicular to the mean magnetic ﬁeld. This
gives


 2
∂y 2
∂y
=
,
(4.39)
VA2∂x∂t

4.4. AMBIPOLAR DIFFUSION

89

Figure 4.4:
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 .

(4.40)The lhs of the equation is given with a dispersion of polarization direction Δφ (rad)




∂y 2
= Δφ2 ,
(4.41)∂xin which the lhs is measured by the wave-pattern of the Alfvén wave. Therefore, the magnetic ﬁeld
strength is given as
ΔV(4.42)
B⊥ = (4πρ)1/2
Δφ

4.3.2

Observation of Mass-to-Flux Ratio

Figure 4.5 shows the mass-to-ﬂux ratio normalized with the critical value


λC ≡


where

M
Φ obs

M
Φ crit(4.43)

1
M=(4.44)
Φ crit 2πG1/2
is the critical mass-to-ﬂux ratio which is obtained numerically from equation (4.31). λC > 1 represents
a supercritical cloud while λC < 1 represents a subcritical cloud.
Figure 4.5 indicates cloud cores are more or less found near the critical mass-to-ﬂux ratio. They
1 (very supercritical) or λC  1 (very subcritical).
are not distributed either in the regions λC

4.4

Ambipolar Diﬀusion

4.4.1

Ionization Rate

In the dense clouds, the ionization fraction is low. Since the uv/optical radiations from stars can notreach the cloud center, potential ionization comes from the cosmic ray particles. In this case the rate

90

CHAPTER 4. LOCAL STAR FORMATION PROCESS

Figure 4.5: Mass-to-ﬂux ratio is plotted against the column density. Dots are obtained from Zeeman
splitting and stars are Chandrasekhar-Fermi method.

Figure 4.6: 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.

4.4. AMBIPOLAR DIFFUSION

91

of ionization per volume is given as ζnn , where ζ ∼ 10−17 s−1 . In Figure 4.6, the ionization fraction forvarious density is shown (Nakano, Nishi, & Umebayashi 2002). This clearly shows that the fraction of
−1/2
8
−3
for nH <ions decrease approximately in proportion to nH∼10 cm . This is understood as follows:+
Equilibrium balance between one kind of ion m and its neutral m0
m0 ↔ m+ + e− ,(4.45)

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 recom8−3
bination rate coeﬃcient. Nakano (1984) obtained in the range of 102 cm−3 <∼n<∼10 cmρi = Cρ1/2
n ,(4.46)

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

4.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 ﬁeld. Thus, the neutral molecules, a major
component of the gas, are coupled with the magnetic ﬁeld indirectly via ionized ions. Equation of
motion for the molecular gas isdun= −c2s ∇ρ − ρ∇ψ + ρrω 2 er + F ,(4.47)
ρdt
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 γ(v i − v n ).

(4.48)We have assumed ρ  ρn since the mass density of the charged component is much smaller than thatof the neutral one. On the other hand, that for ions is1
du=(∇ × B)×B − F ,(4.49)
ρidt
4πwhere we ignored the pressure and self-gravity forces compared with the Lorentz force. If the inertiaof the ions are ignored [lhs of eq.(4.49)=0], equations (4.47) and (4.49) give the equation of motion
similar to that of the one-ﬂuid as follows:1
du= −c2s ∇ρ − ρ∇ψ + ρrω 2 er +(∇ × B)×B.
(4.50)
ρdt
4πIgnoring the inertia of the ions, equation (4.49) indicates that the friction force should be balancedwith the Lorentz force as
1(∇ × B)×B.
(4.51)F =4π
Using equations (4.46) and (4.48), this gives the drift velocity of ions against the neutrals asv D ≡ vi − v n ,F,=ρi ρn γ1(∇ × B)×B,=
4πρi ρn γ
1(∇ × B)×B.
=
3/2
4πCγρn(4.52)

CHAPTER 4. LOCAL STAR FORMATION PROCESS

92

Figure 4.7: Relation central density ρc against mass-to-ﬂux ratio M/Φ. Plane-parallel and axisymmetric models are illustrated in the left and right panels, respectively.
Since the Lorentz force works outwardly, the drift velocity v D also points outwardly. Viewing fromthe ions and magnetic ﬁeld lines, v n − v i points inwardly. Thus, the mass in a ﬂux tube varies withtime; that of the cloud center increases while that of the cloud surface decreases. The ambipolar
diﬀusion is a process to change the mass distribution against the magnetic ﬂux tube.
It is useful to derive the timescale of ambipolar diﬀusion. The characteristic timescale of magnetic
ﬂux loss is given as τA = L/|v D |, where L represents the size of the cloud. This gives
τA 

4πCγρ3/2 L|(∇ × B) × B|

 3 × 106 yr(4.53)n4
10 H2 cm−3

3/2 

B
30μG

−2 

L
0.1pc

2

.

(4.54)

This ambipolar diﬀusion timescale is several times longer than the dynamical contraction timescale(eq.[2.26]) or free-fall timescale as


tﬀ

=

3π
32Gρ̄

1/2

(4.55)
n 0.5 × 10 yr
4
10 H2 cm−3
6

−1/2

.

(4.56)

Thus, the magnetic ﬂux of the cloud is not be lost during the dynamical contraction ( free-fall).
How is the structure aﬀected by the increase of mass/ﬂux ratio at the cloud center? Here, we will
derive the relation between the mass-to-ﬂux ratio and central density and answer the question. Numerical calculation is necessary for the axisymmetric (realistic) cloud. However, this can be obtained
for the case of a plane-parallel disk in which the magnetic ﬁeld is running parallel to the disk, For
the plane-parallel disk the hydrostatic balance is achieved as
∂
∂z

p+

B2
8π

= −ρ

∂ψ
= −2πGρσ,∂z(4.57)

4.4. AMBIPOLAR DIFFUSION

93



z
ρdz. Equation (4.57) is integrated to obtain the relation between the variables at thewhere σ = −z
center z = 0 as
B2
B2
πG 2
σ ,
(4.58)
pc + c − ps + s =
8π
8π
2 s

where we used dσ/dz = 2ρ. This is rewritten as
1
ρc
+
ρs 8π



Bc
ρc

2

ρs
c2s



ρc
ρs

2

1
−
8π



Bs
ρs

2

πGσs2
ρs
=
+ 1.
c2s
2ρs c2s(4.59)

1/2

In the second term of lhs, yc ≡ (ρc /Bc )(cs /ρs ) is a nondimensional mass ﬂux ratio at the center andthis is related to the center-to-surface density ratio ρc ≡ ρc /ρs . Using these nondimensional variables,
equation (4.59) is expressed as(ρ )2
1= ρ̄c ,(4.60)
ρc + c 2 −
8πyc
8πys2
where ρ̄c ≡ πGσs2 /(2ρs c2s ) + 1 represents the central density necessary for the disk to be supported
without a magnetic ﬁeld (y → ∞). Figure 4.7 (left) illustrates this relation. This shows that if thecentral mass-to-ﬂux ratio increases (moving upward), the central density increases. In low density,ρc  ρ̄c , equation (4.60) indicates that ρc ∝ yc ∝ ρc /Bc . This means that on the low density regime,
the density increases but the strength of magnetic ﬁeld does not increase. Equation (4.60) also
indicates that yc increase much as yc → ρ̄c ys when ρc → ρ̄c . This relation between mass-to-ﬂux ratio
and the central density is well ﬁtted to the numerical calculation of plane-parallel disk cloud driven
by the ambipolar diﬀusion by Mouschovias, Paleologou, & Fiedler (1985). Does this indicate that the
mass-to-ﬂux 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 ﬁrst studied by Nakano (1979)
adopting a method seeking for magnetohydrostatic conﬁguration (Mouschovias 1976a, 1976b). He
obtained that the magnetic ﬂux to mass ratio near the center of the cloud does not decrease much
below a critical density above which magnetohydrostatic conﬁgurations no more exist. Paleologou &Mouschovias (1983), Shu (1983), and Mouschovias, Paleologou, & Fiedler (1985) gave a completely
diﬀerent result. That is, the magnetic ﬂux to mass ratio near the center decreases much when the
6−31density becomes larger than n>∼10 cm . Critical discussion has been done between them . However,
we have now realized that the diﬀerence comes from the geometry assumed.
Relation of the mass-to-ﬂux ratio to the central density for the axisymmetric realistic (not planeparallel disk) cloud is plotted in Figure 4.7 (right) (Tomisaka, Ikeuchi, & Nakamura 1990). The
central mass-to-ﬂux ratio increases by the ambipolar diﬀusion. This increases the central density.
Lines with arrows indicate the evolutionary paths driven by the ambipolar diﬀusion. It is to be
noticed that the relation is completely diﬀerent from that of the plane-parallel disk. Increase of
the mass-ﬂux ratio is small although the central density increases much. This indicates that in thequasistatic evolution driven by the ambipolar diﬀusion the mass-ﬂux ratio does not increase much incontrast to the plane-parallel case.
In fully ionized plasma, the magnetic ﬁeld is coupled to the matter. However, in the low ionizedgas the magnetic ﬁeld is coupled only to the ions. Thus the induction equation of magnatic ﬁeld
1
One of the reasons why large increase in the mass/ﬂux ratio is favored is understood as follows: There havebeen a long-standing “magnetic ﬂux problem of stars” in which the magnetic ﬂux of, say, 1M main-sequence starΦ ∼ XXG cm2 is much smaller than that of the parent cloud Φ ∼ XXG cm2 . The magnetic ﬂux must be reduced in
the star formation process. If the mass-ﬂux ratio increases much in these density range, this might resolve the magnetic
ﬂux problem of stars.

CHAPTER 4. LOCAL STAR FORMATION PROCESS

94
should be written∂B+ ∇ × (vi × B) = 0.
∂t(4.61)

This leads to the ﬁnal expression as∂B+ ∇ × (vn × B) + ∇ ×∂t

1
3/2 [(∇ × B) ×B] × B
4πCγρn= 0,(4.62)

where we used the expression of the drift velocity [equation (4.52)]. In the isothermal cloud, there
is a maximum mass supported by the magnetic ﬁeld. Since in the subcritical cloud M < Mcr themagnetic ﬂux escapes from the center, the critical mass decreases with time. After the critical mass
becomes smaller than the actual mass (Mcl > Mcr ), there is no hydrostatic conﬁguration. The cloud
evolves into the supercritical cloud region and experiences dynamical collapse.

4.5

Dynamical Collapse

Figure 4.8: Evolution of isothermal clouds massive than the Bonnor-Ebert mass. Density (left) andvelocity (right) distributions are illustrated. Solid lines show the cores of the preprotostellar phase(prestellar core) and dashed lines show those of protostellar phase (protostellar core). The evolutionof the protostellar phase is studied by the sink-cell method, where we assume the gas that entered in
the sink-cells is removed from ﬁnite-diﬀerence grids and add to the point mass sitting at the centerof the sink-cells which corresponds to a protostar.
In 1969, Larson (1969) and Penston (1969) found a self-similar solution which is suited for the
dynamical contraction. Figure 4.8 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.
Although the ﬁgure 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 asfollows:

4.5. DYNAMICAL COLLAPSE

95

1. The cloud evolves in a self-similar way. That is, the spatial distribution of the density (left) at
t = 0.326 is well ﬁtted by that at t = 0.293 after shifting in the −x and the +y directions. As
for the infall velocity spatial distribution, only a shift in the −x direction is needed.
2. The density distribution in the envelope, which is ﬁtted 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 deﬁned 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 ﬁgure, at t − t0 = 0.003 (t = 0.326) the central density
reaches ρc  104 and at t − t0 = 0.036 (t = 0.293) the density is equal to ρc  102 . This shows
the maximum (central) density is approximately proportional to (t − t0 )2 , which is reasonablefrom the description of the free-fall time ∝ ρ−1/2 .
The basic equations of spherical symmetric isothermal ﬂow are1 ∂∂ρ
+ 2 (r 2 ρv) = 0,∂tr ∂r∂v c2is ∂ρ GM
∂v
+v
+
+ 2
= 0,
∂t
∂r
ρ ∂r
r

(4.63)(4.64) r

M (r, t) = M (0, t) +4πr 2 ρ(r  , t)dr  ,(4.65)

0where 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(ξ)
M (r, t) = is
G
r
,
ξ =
cis |t − t0 |ρ(r, t) =(4.66)(4.67)(4.68)(4.69)

should be found, where Ω and V are functions only on ξ. For example, equation (4.66) asks theshape of the density distribution is the same after resizing of equation (4.69) r → r/cs |t − t0 | andre-normalizing ρ → ρ · 4πG(t − t0 )2 as Ω(ξ). Since
d
∂
=
∂t
dξ



∂ξ
∂t



and



=
r

r
d
,
2
cis |t − t0 | dξ



d ∂ξ
d
1
∂
=
,
=
∂r
dξ ∂r t cis |t − t0 | dξ
the basic equations for the spherical symmetric model yieldm = (ξ − V )ξ 2 Ω,

 dV

(4.70)

(4.71)(4.72)
2(ξ − V ),ξ

(4.73)

2(ξ − V )2 − 1 dΩ= Ω − (ξ − V ) (ξ − V ),Ωdξ
ξ

(4.74)(ξ − V )2 − 1dξ= Ω(ξ − V ) −CHAPTER 4. LOCAL STAR FORMATION PROCESS

96

Equations (4.73) and (4.74) have a singular point at which (ξ−V )2 −1 = 0 or V = ξ±1. Since the pointof ξ =const moves with cis , the ﬂow 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 ﬂow
has to pass the sonic point smoothly, the rhs of equations (4.73) and (4.74) have to be equal to zeroat the singular point ξ∗ . This gives at the sonic pont ±ξ∗ > 0,Ω−22(ξ − V )= Ω(ξ − V ) − = 0,ξξ(4.75)

which leads toV∗ = ξ∗ ∓ 1,

(4.76)

2.
ξ∗(4.77)

Ω∗ = ±

These equations (4.72), (4.73) and (4.74) have an analytic solution
V (ξ) = 0, Ω =2, m = 2ξ, −∞ < ξ < ∞ξ2(4.78)

This is a solution which agrees with the Chandrasekhar’s SIS. Generally, solutions are obtained only
by numerical integration. |ξ → ∞| the solution have to converge to an asymptotic form of
A − 2 V∞ 4V∞ + (A − 2)(A − 6)+ 2 ++ O(ξ −4 ),ξξ6ξ 3
Ω∞ (A − 2)A−
+ O(ξ −6 ),
2
ξ
2ξ 4

V (ξ) = V∞ −(4.79)

Ω(ξ) =(4.80)

This shows that for suﬃciently large radius the gas ﬂows with a constant inﬂow 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 4.9 (left) and Figure 4.9 (right) witht < 0. This time evolution is expected from the self-similar solution. This shows thatρ ρc(in the central region)∝ r −2 (in the outer envelope)(4.81)v∝r(in the central region) 3.28cis (in the outer envelope)(4.82)

Reaching the outer boundary the numerical solution (Fig.4.8) diﬀers from the self-similar solution(Fig.4.9). For example, v is reduced to zero in the numerical simulations, while it reaches a ﬁnite 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.

4.5.1

Inside-out Collapse Solution

In 1977 Shu found another self-similar solution which is realized after a central protostar with inﬁnitesimal mass is formed in the singular isothermal sphere solution. The gas begins to accrete to theprotostar. Outside the region where the accretion occurs, the initial SIS is kept unchanged, since the4.5. DYNAMICAL COLLAPSE

97

16

3000

t=0y

14

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

log n [m-3]

2000
4

+105
10

-105
+106

8

+104 +105 +106
1000
t=0y

-106

-104
6
13

14

15
16
log r [m]

17

0
13

14

-105

-106

15
16
log r [m]17

Figure 4.9: A self-similar solution indicating a dynamical collapse of isothermal spherical cloud(Larson-Penston solution). Spatial distribution of the density and inﬂow velocity which are expectedfrom the self-similar solution are plotted. Dashed lines show the evolution prestellar core and solidlines show that of protostellar core. Taken from Hanawa (1999).

CHAPTER 4. LOCAL STAR FORMATION PROCESS

98

SIS is a hydrostatic solution. And the front of accretion expands radially outward in time. Since theinﬂow region expand outwardly, he called it the inside-out collapse solution. In Figure 4.10, theevolution is shown for density and inﬂow velocity. This solution gives the evolution of a protostellar
core formed in a subcritical cloud/cloud core.
Using equations (4.79) and (4.80) and assuming the inﬂow velocity should reduce at large radius,
we obtain
A − 2 (A − 2)(A − 6)−
+ ···,(4.83)V (ξ) = −ξ6ξ 3
Ω(ξ) =

A(A − 2)
A−
+ ···.
ξ22ξ 4(4.84)

Since ξ → ∞ means t → t0 (if r is ﬁnite), Ω(ξ) → A/ξ 2 means thatρ(r, t0 ) =Ac2is.
4πGr 2(4.85)

Comparing with the SIS, when A = 2 this gives the SIS and when A > 2 this gives a density
distribution in which the pressure is ineﬃcient and the cloud is contracting. The solution with A > 2
is obtained by a procedure as (1) at a suﬃciently large radius ξ1 , calculate V (ξ1 ) and Ω(ξ1 ). (2) fromthese values, integrate equations (4.73) and (4.74) inwardly. Figure 4.10 show the solution of thistype. The solution with A > 2 inﬂow 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 singular line2
V = −ξ + 1 at ξ = 1 (V = 0) . 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 inﬁnitesimally small amount of excess mass at thecenter of SIS, the accretion begins from the center while outside a radius the cloud is left static. The
1/2 and Ω ∝ ξ 3/2 .
inner part of the solution ξ <∼1, V and Ω are well expressed as V ∝ ξ
1/2
3/2
and Ω ∝ ξ
are explained as follows: Conservation of
The power-law distributions of V ∝ ξ
the total energy of an inﬂowing gas shell is expressed as
GMr0
v 2 GMr
−
=−
,
2
r
r0

(4.86)

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, this gives


v

2GMr
r

1/2

∝ r −1/2 ,(4.87)M∗ , thatis, Mr =where we assumed a major part of Mr comes from the mass of a protostar
r
r
r
2
2
M∗ + 0 ρ4πr dr  M∗ Since the average density inside the radius r ρ̄(< r) = 0 ρ4πr dr/ 0 4πr 2 dr =3ρ(r) for SIS distribution, the time necessary for a gas shell to reach the center is proportional to∝ ρ̄(< r)−1/2 ∝ ρ(r)−1/2 ∝ r 1 ∝ (t − t0 )1 , where we used a fact that the front of accretion expands
with a constant speed cs is. 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πr 2 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πr 2 ρv = const.
2

This means A = 2 +  and  > 0 and → 0(4.88)

4.5. DYNAMICAL COLLAPSE

99

Figure 4.10: 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 correspond to 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 ).
Since the mass swept by the sound wave per unit time is equal to Ṁsw = 4πr 2 ρSIS cs = 2c3s /G, theaccretion rate is proportional to c3s /G. Shu (1977) obtained the accretion ratec3
Ṁ = 0.975 s ,G(4.89)

for SIS A = 2+ . Using this equation, equations (4.87) and (4.88) indicate that the spatial density
distribution is expressed by a power-law as
ρ(r) ∝ r −3/2 ,(4.90)

which is valid for the region except for the vicinity of the front of the accretion.

4.5.2

Protostellar Evolution of Supercritical Clouds

What is a protostellar core formed in a supercritical cloud/cloud core? Is this diﬀerent from theinside-out solution of Shu (1977)? A solution which corresponds to the protostellar core is obtainedby Hunter (1977) and Whitworth and Summers (1985). This is a solution with t > t0 in equation(4.69). The asymptotic behaviors of the density and infall velocity reaching the center are diﬀerent
from that of the Larson-Penston self-similar solution for a prestellar collapse. That is,Ω→ ﬁnite(LP)→ inﬁnite (Inside − out)

(4.91)

V→ ﬁnite(LP)→ inﬁnite (Inside − out)

(4.92)

Using the boundary conditions suitable for the inside-out type solution, another self-similar solution
is obtained. In Figure 4.9, such kind of solution is also plotted for t > 0.

100

CHAPTER 4. LOCAL STAR FORMATION PROCESS

Table 4.1: Ages of stellar population observed for respective star forming regions.

4.6. STAR FORMATION TIME SCALE

101

Take notice that the solutions of t < t0 (prestellar) and t > t0 (protostellar) agree with eachother at t = t0 . Even if the boundary conditions at the center for the similarity variables, V andΩ, are completely diﬀerent, the diﬀerence between the two is small in the physical variable v andρ. Therefore, the evolution of a supercritical core is thought to be expressed by the Larson-Penston
self-similar solution extended to the protostellar core phase by Hunter (1977) and Whitworth andSummers (1985).
Assume that we observe a protostellar core and obtain their density and infall velocity spatial
distributions. Can we distinguish which solution is appropriate for the Shu’s inside-out solution or theextended Larson-Penston solution? This seems hard, because the structure of density and velocity
distributions are similar after the protostar is formed: the density and velocity show almost similar
power-law as ρ ∝ r −3/2 and v ∝ r −1/2 irrespective of the inside-out solution or the extended LarsonPenston solution. The region where the infall velocity is accelerated toward the center (accretiondominated region) is expanding after the protostar is formed. Therefore, to distinguish between thetwo solutions becomes harder and harder after the protostar is formed. The diﬀerence would be large
and we would have a deﬁnite answer which solution is appropriate to describe the cloud collapse, if we
can observe a very young protostellar core or a preprotostellar core which shows dynamical collapse.
However, since the time-scale of such a phase is much shorter than the evolved protostellar phaseor a younger preprotostellar core, the number of such kind of objects would be small (τﬀ ∝ ρ−1/2 ).
Therefore, we are looking for such objects just before or after the protostar formation.

4.6

Star Formation Time Scale

There are long-standing controversy on how fast stars are form in the molecular clouds. One idea isthat the interstellar cloud is in a hydrostatic balance in which magnetic ﬁeld plays an important role.
That is, the cloud is magnetically subcritical. In such a cloud, the magnetic ﬁeld escapes from thecenter in a timescale of τD ∼ 10 × tﬀ (see §4.4). After the mass-to-ﬂux ratio exceeds a critical value,
dynamical contraction begins. The timescale of cloud lifetime is equal to the time scale of ambipolar
diﬀusion, τD ∼ 3×107 yr(ρ/102 H2 cm−3 )−1/2 . Long lifetime of molecular clouds is necessary to explain
giant molecular clouds are formed by agglomeration process of small clouds. The necessary lifetime
is estimated 107 yr – 108 yr.
Another idea is that the star formation timescale is much shorter than this timescale. Table 4.1(Hartmann 2009) shows the ages of stellar population observed for respective star forming regions.
Stellar clusters and associations accompanied with molecular clouds have ages younger than 5 Myr.
And clusters and associations without accompanied molecular clouds are older than 5 Myr. Thismeans that a molecular cloud will disappear after 5 Myr since star formation begins in the cloud.
Since this timescale is shorter than τD , the clouds seem to continue to be formed and the turbulance
supports the cloud, the cloud core begins dynamical contraction after the turbulance decays away
within the sound crossing timescale ∼ Myr.

4.7

Accretion Rate

Using equation (2.26), the necessary time for a mass-shell at R to reach the center (free-fall time) isexpressed as
1/2R3π(4.93)
T (R) ≡
2GM (R)
2(for detail of this section see Ogino, Tomisaka, & Nakamura 1999).

102

CHAPTER 4. LOCAL STAR FORMATION PROCESS

Figure 4.11: Mass accretion rate against the typical density of the cloud.
Consider two shells whose initial radii are R and R + ΔR. The time diﬀerence for these two shells
to reach the center ΔT (R) can be written down using equation (4.93) as
R dM (R)
3
πR1/2
−
ΔR.
ΔT (R) = 3/2
1/2
2 2M (R) dR
2 (GM (R))

(4.94)

Mass in the shell between R and R + ΔR, ΔM ≡ M (R + ΔR) − M (R) = (dM/dR)ΔR, accretes
on the central object in ΔT (R). Thus, mass accretion rate for a pressure-free cloud is expressed as
ΔM/ΔT . This leads to the expression asdM (R)dR.
dM (R)
2M (R) dR
R

23/2 G1/2 M (R)3/2
dM
M (R)
(R) =
3/2
3
dT
π
R
− R
2

(4.95)

This gives time variation of the accretion rate. Consider two clouds with the same density distribution ∂ log ρ/∂r but diﬀerent absolute value. Since these two clouds have the same ∂ log M (R)/∂ log R,
the mass accretion rate depends only on M (R)/R, and is expressed asdM(R) ∝ M (R)3/2 .
dT(4.96)

This indicates that the accretion rate is proportional to ρ3/2 , while the time scale is to ρ−1/2 . This is
conﬁrmed by hydrodynamical simulations of spherical symmetric isothermal clouds (Ogino et al.1999).
When the initial density distribution is the SIS as ρ ∝ r −2 , the mass included inside R is proportional
to radius M (R) ∝ R. In this case, equation (4.95) gives a constant accretion rate in time. In Figure
4.11 we plot the mass accretion rate against the cloud density. α represents the cloud density relative
to that of a hydrostatic Bonnor-Ebert sphere. This shows clearly that the mass accretion rate is4.8. OUTFLOW

103

Figure 4.12: A model proposed to explain time variation in accretion rate by Henriksen, André, &Bontemps (1997). The density distribution ρ(r) at t = 0 (left) and expected accretion rate (right).
proportional to α3/2 for massive clouds α > 4. This is natural since the assumption of pressure-lessis valid only for a massive cloud in which the gravity force is predominant against the pressure force.
Similar discussion has been done by Henriksen, André, & Bontemps (1997) to explain a decline
in the accretion rate from Class 0 to Class I IR objects. They assumed initial density distribution ofρ⎧⎨ = ρ0(r ≤ rN ),⎩ = ρ0r
rN

−2/D1(r > rN ),(4.97)as shown in Figure 4.12. Since the free-fall-time of the gas contained in the inner core r ≤ rN is thesame, such gas reaches the center once. It makes a very large accretion rate at t = (3π/32GρN )1/2
3 δ(t − (3π/32Gρ )1/2 ). If D = 1, ρ ∝ r −2 for r >r . Since M ∝ R1 and t ∝ R,
as Ṁ = (4π/3)ρ0 rN
N
1
ﬀ
∼ N
equation (4.95) predicts Ṁ ∝ R0 ∝ t0 . A constant accretion rate is expected for this power-law andthe accretion rate is converged to a constant value after the stellar mass is much larger than than2
M (rN ). If D1 = 2, ρ ∝ r −1 for r >that was containd in rN , M∗
∼rN . Since for this power M ∝ R
and tﬀ ∝ R1/2 , equation(4.95) predicts Ṁ ∝ R3/2 ∝ t3 . They gave Ṁ ∝ t3D1 −3 for D1 > 2/3.
ProblemShow that Ṁ ∝ t3D1 −4/3D1 if ρ ∝ r −2/D1 and D1 < 2/3 [Henriksen, André, & Bontemps (1997)].

4.8

Outﬂow

In chapter 1, we mentioned ‘outﬂows’ ejected from the protostars and pre-main-sequence stars. L1551
IRS5 is a typical example which shows a number of outﬂows are ejected in the course of star formation.
One is molecular outﬂow, whose lobes extends L ∼ 0.5(D/160pc)pc in two opposite directions fromthe IRS-5 (Snell et al 1980). This is traced by CO J = 1 → 0 emission line. They estimated an
expansion velocity of the CO gas of VCO ∼ 15km s−1 , a mass of 0.3 M . The dynamical age is equal toτdyn = L/VCO ∼ 3 × 104 yr. Inside the CO outﬂow lobe, several Herbig-Haro objects are found, which

CHAPTER 4. LOCAL STAR FORMATION PROCESS

104

2

1.5

1

0.5

0
0

0.5

1

1.5

2

Figure 4.13: Magneto-centrifugal wind model. Figure shows the eﬀective potential for a gas parcel
with conserving the angular rotation speed ω0 of the disk (eq.4.99). 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.
are emission-line nebulosities. Proper motion studies have revealed the motion of such H-H objects.
Cudworth & Herbig (1979) reported HH28 and 29 have 0.15 − 0.20arcsec yr−1 , which corresponds
to 150 − 200km s−1 . In a neighbor of the IR source a stellar jet is found (Mundt & Fried 1983) by
optical emission lines. Its size is equal to 17arcsec = 4 × 1016 cm(D/160pc) and the jet indicatesrather well-collimated shape (the opening angle  10◦ ). The relationship between these two outﬂows(massive molecular outﬂow and less massive optical jet) are not clear yet.

4.8.1

Magneto-driven Model

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 conﬁguration.
Figure 4.13 explains the magneto-centrifugal wind mechanism (Fig. 1 of Blandford & Payne 1982).
Consider a disk under the gravity of a central star with mass M∗ . Assume that gas in the disk at theradius r0 rotates with a Kepler speed ofvφ (r0 ) =

GM∗
r0

1/2

.

(4.98)

Consider a gas parcel whose original position was (r0 , z = 0). The eﬀective potential for the particle
rotating with an angular rotation speed of ω0 = vφ (r0 )/r0 is written as


1
GM∗
1
r0
− r 2 ω02 = −+φ(r, z) = −
1/2
1/2
2
r0
2
(r 2 + z 2 )(r 2 + z 2 )
GM∗

3/2



r
r0

2 

,

(4.99)

where ω0 = (GM )1/2 /r0 . (Be careful that ordinary eﬀective potential is calculated for constantangular momentum. However, this is for constant angular rotation speed.) The gas will co-rotate

4.8. OUTFLOW

105

with the same angular rotation speed as the disk, if the magnetic ﬁeld is suﬃciently strong. Equation(4.99) corresponds the eﬀective potential for the gas co-rotating with the disk at r0 . Figure 4.13
plots the isopotential contour lines for the eﬀective potential [eq.(4.99)]. The contour line passing thepoint (r0 , 0) represents Φ = −3GM/2r0 . The angle between the isopotential line of Φ = −3GM/2r0
and the r-axis is equal to ±60◦ (see problem below). Further, it can be obtained that the eﬀective
potential increases upwardly while it decreases with reaching the central star or radially outwardly.
This leads to an important conclusion that
1. If the angle between the magnetic ﬁeld line and the r-axis θmag is larger than 60◦ (θmag > 60◦ ),gas parcel moving along the ﬁeld line must climb the eﬀective potential. It needs extra energy
to depart from the disk.
2. If θmag < 60◦ , gas parcel moving along the ﬁeld line slides down the eﬀective potential. It can
depart from the disk without extra energy.
Be careful that this is valid if the magnetic ﬁeld is suﬃciently strong. If the Kepler disk is threadedwith the magnetic ﬁeld lines with a shallow angle, and if the angular momentum is suﬃciently
transferred along the magnetic ﬁeld line, in this case, the gas is accelerated by the extra centrifugal
force and ﬁnally escapes from the potential well of the central star. This mechanism is called the“magneto-centrifugal wind” mechanism.
ProblemShow that the angle between the isopotential line which passes the point (r0 , 0) and the r-axis isequal to ±60◦ . Calculate the Taylor expansion of the eﬀective potential near the point (r0 , 0) andshow that to satisfy φ(r0 + Δr, Δz) = φ(r0 , 0), (Δz/Δr)2 must be equal to 3.
Angular Momentum Transfer
In the axisymmetric case, the poloidal (z, r) and toroidal (φ) components of the magnetic ﬁeld B andcurrent j are decoupled with each other. That is, the poloidal (Bp ) and toroidal (Bφ ) magnetic ﬁelds
are made by the toroidal (jφ ) and poloidal (jp ) electric currents, respectively. As for the Lorentz force
j × B, the poloidal component comes from jφ Bp or jp Bφ , while the toroidal component does fromj p ×B p . Even if there is no toroidal magnetic ﬁeld (thus no poloidal electric current), there exists thepoloidal component of the Lorentz force, which acts as a pressure to counter-balance the self-gravity(§4.2). On the other hand, the toroidal component of the Lorentz force appears only the case withthe poloidal electric current and thus toroidal component of magnetic ﬁeld Bφ . This means thatthe angular momentum is transferred by the magnetic ﬁeld only when Bφ exists. Equation (C.14)
explains how the angular momentum density ρvφ r is transferred. The left-hand side of equation(C.14) represents the advection of the angular momentum density, while the right-hand side


∂Bφ
r 1 ∂
(rBφ )Br +Bz =4π r ∂r
∂z



Bp
4π



· ∇(rBφ ),(4.100)represents the torque exerted on the gas parcel.
The induction equation of the magnetic ﬁeld [eq.(C.17)] shows that Bφ is generated from poloidal
magnetic ﬁeld by the eﬀect of rotational motion vφ . This indicates that the angular momentum istransferred as follows:
1. The toroidal component of magnetic ﬁeld Bφ is ampliﬁed by the dynamo process [equation(C.17)]. The magnetic ﬁeld lines run like spiral viewing from the top like Figure 4.14.

106

CHAPTER 4. LOCAL STAR FORMATION PROCESSSide View
Pole-on View

Figure 4.14: Pole-on (left) and side (right) views of the magnetic ﬁeld lines, which drive the outﬂow.
2. Negative torque works near region attached to the fast-rotating disk, while positive one appears
in the region connected to the radial magnetic ﬁeld. The angular momentum is transferredfrom inside to outside along the same ﬁeld line.
3. Angular momentum received by the gas near disk surface leads the excess centrifugal force which
accelerate the gas. The disk inside which is aﬀected by the negative torque loses its angularmomentum and insuﬃcient centrifugal force drives further accretion.
This angular momentum transfer on the same magnetic ﬁeld has a character of angular momentum
redistribution.
Centrifugal Radius
In a diﬀuse cloud, the centrifugal force does not play an important role. Consider a gas parcel whose
speciﬁc angular momentum is equal to j. When the gas element contracts till the distance r fromthe center, the centrifugal force per mass of
vφ2
j2
= 3,
r
r

(4.101)works. Contraction stops when the gravitational acceleration GM/r 2 is balanced by the centrifugal
force. Its radius (centrifugal radius) is expressed as
rc =

j2
.
GM

(4.102)

Assuming that the disk is near the hydrostatic balance in z-direction, total column density is related
to the volume density on the z = 0 plane asσ=

2c2s ρcπG

1/2.

(4.103)

We can show that a nondimensional ratio of speciﬁc angular momentum to the mass is equal to theratio of the free-fall time to the rotation period as follow:
cs ωc r 2
ωc
cs j
=
=
.
2
GM
Gσπr
(2πGρc )1/2

(4.104)

4.9. EVOLUTION TO STAR

107

Figure 4.15: (Left) Relation between density and temperature at the center of the cloud. (Right)
Cooling and heating rates, which are calculated for the cloud center are plotted against the density.
Many numerical simulations conﬁrm that the contracting disk in the runaway isothermal contraction
phase has some ‘universality.’ The ratio of the free-fall time to the rotation period is approximately
equal toωc
 0.3,(4.105)(2πGρc )1/2
regardless of the initial conditions (Matsumoto, Nakamura, & Hanawa 1997). Thus, cs j/GM  0.3.
Finally, we can see the centrifugal radius as
rc =

GM
c2s



cs j
GM

2

 0.3

GM
,
c2s

(4.106)

increases with time in proportional to the mass M , because gas element with large j contracts later.
As the centrifugal radius increases with time, the launching point of the outﬂow also expands with
time.

4.8.2

Entrainment Model

Models that a protostar is driving a jet and the jet transfers its momentum to the ambient material
and forms a molecular outﬂow are called as entrainment models. From hydrodynamical calculations,
De Young (1986) pointed out there are two mechanisms of entrainment. The ejected jet forms a bow
shock and matter hit by the bow shock is accelerated promptly (prompt entrainment). His numericalsimulations predicted that gas within the radius from the jet axis of (1 − 3)× jet radius is entrainedin this mechanism. After the bow shock has passed away, slower mechanism of entrainment works,
in which a turbulent layer is developed around the jet and the Reynolds stress in the layer transfers
the momentum (Stahler 1994). In this mechanism, the ratio of entrained ambient gas to the injectedmass by the jet is estimated as (1/3 − 3) (De Young 1986).

4.9

Evolution to Star

In the proceeding sections we have seen the evolution from gas cloud to stars. However we wererestricted to the isothermal gasses. To understand the change in the temperature, we have to consider the radiation which keeps the gas isothermal. Spherical symmetric radiative hydrodynamical
simulations have been studied to understand the evolutionary path from interstellar cloud to star.

108

CHAPTER 4. LOCAL STAR FORMATION PROCESS

Figure 4.16: Dynamical evolution of a spherical cloud studied with RHD simulation. Density (left)
and temperature (right) distributions are plotted against the radius.

4.9. EVOLUTION TO STAR

109

Figure 4.15 (left) illustrates the relation between ρ and T obtained by Masunaga & Inutsuka (2000).
Several characteristic power-laws are seen in the ﬁgure. Figure 4.16 plots the density and temperature
distribution.
Isothermal Phase−13−30
Below ρ<∼10 g cm , gas is essentially isothermal T ∝ ρ . This corresponds to points number 1-3
of Figure 4.16. Since internal energy of gas is transferred to the thermal energy of dusts by collisions,the main coolant in this regime is the dust thermal radiation. The cooling rate per mass isΛ = 4κp σT 4 ,(4.107)

where κp (10K) ∼ 0.01cm2 g−1 and σ represent Planck mean absorption coeﬃcient and the Stephan−14−3 is cosmic-ray heating:Boltzman constant. Main heating process for ρ<CR ∼ 2 ×∼10 g cm
−4
−1
−1
10 erg s g (Goldsmith & Langer 1978). Balance between these two majors asks that temperature
is constant as1/4 
1/2
κp
CR
)
)
.
(4.108)T ∼ 3K2 × 10−4
0.01cm2 g−1
Heating rate due to dynamical compression,Γ = −pd(1/ρ)
 c2s (4πGρ),
dt(4.109)increases according to the contraction and it balances with the above cooling at the density ρA ∼10−14 g cm−3 (Masunaga, Miyama, & Inutsuka 1998).
The cloud in this phase experiences the dynamical contraction as described in section 4.5. Structure of ρ(r) and vr (r) is well represented by the Larson-Penston self-similar solution.
First Core−13
−3 ρ<10−9 g cm−3 ,After ρ>∼ ∼∼ρA , gas is no more isothermal and becomes adiabatic. Between 10 g cm <,theoptical
depth
for
the
thermalradiationgas obeys the γ = 5/3 polytropes. Above the density ρ>ρA∼exceeds unity τ >1andradiative
coolingcan
notcompensate
thecompressional
heating.
Aslong as∼<the temperature is low as T ∼100K, neither rotation nor vibration is excited for H2 molecule. Even
H2 gas behaves like single-atom molecule. Thus γ  5/3.
−7.5 g cm−3 , the exponent becomes γ = 7/5, which characterizes that
Between 10−9 g cm−3 <∼ρ<∼10the gas consists of two-atom molecule H2 .
In this phase, relatively large gas pressure supports against the gravity and the cloud becomeshydrostatic (points number 4-6 of Figure 4.16). This is called as “ﬁrst core” made by the molecularhydrogen. The density structure of the ﬁrst core is well represented by a polytrope sphere with thespeciﬁc heat ratio of γ = 7/5 or the polytropic index n = 2 12 . From equation (D.11) in Appendix
D.1, such a polytrope has a mass-density relation as,
Mc1 ∝ ρ1/10
c(4.110)where Mc1 and ρc represent, respectively, the mass of the ﬁrst core and the central density. At the
beginning, the core mass is equal to Mc1  0.01M . As long as the mass increases a factor 3, thecentral density increases 5 orders of magnitude.

CHAPTER 4. LOCAL STAR FORMATION PROCESS

110

Figure 4.17: Explanation of nested grid method.
Second Collapse−2.5 g cm−3 , the exponent becomes γ = 1.1. In this density range of ρ>10−7.5 g cm−3
For 10−7.5 g cm−3 <
∼
∼ρ<
∼10
the temperature exceeds 103 K. Due to this temperature, the hydrogen molecules begin to dissociate
and the dissociation process absorbs the energy. By this, gas becomes much softer (γ → 1) than theﬁrst core. This introduces another dynamical contraction, the second collapse (points number 6-7 ofFigure 4.16).

Second Core
After the dissociation is completed, the exponent becomes γ = 5/3. This gas forms a second corecomposed of atomic hydrogen, which is called the second core (points number 7-8 of Figure 4.16).
Since the accretion rate [eq.(4.89)] is proportional to c3s or T 3/2 , the accretion rates onto the second
core Ṁ2nd is larger than that of the ﬁrst core Ṁ1st as3/2

T /103.5 K
Ṁ2nd 5 × 103.
Ṁ1st
(T /10K)3/2(4.111)Thus, the ﬁrst core disappears quickly and after that the gas begins to accrete onto the second
core, which will be a protostar. However, this is the case of non-rotating, spherical symmetric cloud
collapse. As previously seen, the angular momentum plays a crucial role and forms a disk. The
evolution must be diﬀerent completely for such a case.

4.10

Example of Numerical Simulation

We have described the cloud run-away collapse and succeeding accretion process In the latter phase,
there is a possibility that the outﬂow is driven by the eﬀect of magnetic Lorentz force. A restricted
number of numerical simulation can simulate such evolution throughout. The size of the molecularcores ∼ 0.1pc is more than 3 × 104 times as large as the typical size of the ﬁrst core ∼ 1AU.

4.10. EXAMPLE OF NUMERICAL SIMULATION
(a)

111
(b)

(c)

Figure 4.18: Ωc = 10−13 rad s−1 , Bc = 10μG, ρc = 104 H2 cm−3 , and cs = 190m s−1 (T /10K). Snapshots at the time of 0.6066τﬀ = 1.061Myr represented in diﬀerent 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 ﬁeld lines (dottedlines) and isodensity contours (solid lines) are presented as well as the velocity vectors by arrows.
Therefore numerical scheme to calculate such process must have a large dynamic-range. In the ﬁnitediﬀerence scheme, the spatial dynamic range is restricted by the number of the cells. Although at
least 3 × 105 × 3 × 105 grid points (in 2D) seems necessary to resolve a factor of 3 × 104 , to solve the
MHD equation with the Poisson equation of the self-gravity have not been done. This is done by
the Eularean nested-grid simulation. This method uses a ﬁnite number of grid systems with diﬀerent
spatial resolutions are prepared. 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 ﬁne-scale structure appeared near
the center. Figure 4.17 explains grid cells of the nested grid method. The size of the cells of the n-thlevel grid (Ln) is taken equal to a half of those of the Ln − 1.
A cylindrical cloud with coaxial magnetic ﬁeld in hydrostatic balance is studied. Rotation vector
ω and the magnetic ﬁeld 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 like shown in Figure
4.15, which consists of a number of power-laws for some density ranges.
p=⎧⎨ c2s ρ⎩ c2s ρA

(for ρ < ρA ),

ρ 7/5
ρA

(for ρ > ρA ),(4.112)where we take ρA = 1010 H2 cm−3 . Figure 4.18 is the model of Ωc = 10−13 rad s−1 , Bc = 10μG, ρc =

CHAPTER 4. LOCAL STAR FORMATION PROCESS

112

104 H2 cm−3 , cs = 190m s−1 (T /10K), and scale-length H = 0.341pc(cs /190m s−1 )(ρc /104 H2 cm−3 )−1/2 .
This is a snapshot of 0.6066τﬀ = 1.061Myr. The cylindrical cloud breaks into pieces with a lengthin z-direction equal to the most unstable wavelength of the gravitational instability (Matsumoto,
Nakamura, & Hanawa 1994) as
λmax 



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

0.72 (1 + α/2 + β)

1/2

− 0.6

cs
.
(4πGρc )1/2

(4.113)The magnetic ﬁelds prevent gas motions perpendicular to the ﬁeld 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 4.18 illustrates a snapshot at the time of
0.6066τﬀ = 1.061Myr. 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 ﬁner spatial resolution(upper right). Since speed, vz , of infalling gas is larger than the sound speed, a number of shock fronts
are formed near z  0.02H and z <∼0.01H. Further inner structure is represented in L10 grid (bottom)
which has 32 times ﬁner resolution than L5. Another discontinuity in density is now forming nearz ∼ 0.0002H. Since gas inside of this discontinuity has the density larger than the critical density
ρ > ρA , in this phase the ﬁrst core (p.109) begins to be formed. In the isothermal collapse phase, acharacteristic power-law for the Larson-Penston self-similar solution is seen in the radial distributionsof density as ρ(r, z = 0) ∝ r −2 , and Bz (r, z = 0) ∝ r −1 .
After t > 0.6067τﬀ the ﬁrst core is formed. Radial infall motion is accelerated toward the surfaceof the ﬁrst core and the rotational motion vφ also increases. This ampliﬁcation promotes the toroidal
magnetic ﬁeld, Bφ . Since the torque is proportional to (∇×Bφ eφ )×B p , the amplitude of the torqueincreases after the ﬁrst core formation.
Structure just after ∼500 yr has passed from Figure 4.18 is shown in Figure 4.18 In τ ∼ 500
yr from the epoch of Figure 4.18, ﬂow is completely changed and outﬂow is launched. Panel (b),which corresponds to a four-times close-up of panel (a), shows clearly that outﬂow is ejected from a
region near r  1.5 × 10−4 H and z  0.5 × 10−4 H. And the outﬂow vectors and the disk has angngle of approximately 45 deg. The magneto-centrifugal wind model favors a small angle between
the magnetic ﬁeld and the disk, θmag < 60deg. This model predicts that the angle between the ﬂow
and the disk is also smaller than 60 deg. Numerical results is not inconsistent with this prediction ofmagneto-centrifugal wind model.
To see which force is driving such outﬂow, amplitudes of three forces are compared: the Lorentz
force F m = (∇ × B)×B/4π, the centrifugal force F c = (ρvφ2 /r)er , and the thermal pressure gradientF t = −∇p. The parallel components to the poloidal magnetic ﬁeld of the three forces are as
Fm ·Bp|B p |(Ustyugova et al. 1999)∇ × B)×B B p·,
4π|B p |
1 Bp· ∇(rBφ )2= −8πr 2 B p |=(4.114)

Fc ·

Bp
Br
= ρvφ2 r ·,|B p ||B p |(4.115)

Fp ·

Bp
Bp
= −∇p ·.
|B p ||B p |(4.116)and4.11. EVOLUTION IN THE H-R DIAGRAM

113

Figure 4.19 illustrates the largest force at each grid point. Centrifugal force dominated region (regionC), which is shown by asterisk (*), overlaps the outﬂow region around point P1 (z = 10−4 H, 2.5 ×
10−4 H). Just radially exterior to this region C, there is another region near point P2 (z = 2 × 10−4 H,
3.5 × 10−4 H) ﬁlled with plus signs where the magnetic force is dominated (region M). Comparing
with the ﬂow vectors in the left, the strong outﬂow coincides with this region M and the above region
C. Therefore, the outﬂow is driven by the centrifugal force and the Lorentz force (pressure gradientby toroidal magnetic ﬁelds).
Diﬀerence between models with diﬀerent magnetic ﬁeld strength but the same rotation speed isillustrated in Figure 4.20. 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 the left column aresnapshots when the central density reaches ρA while panels (d), (e), (f) in the right column show
the structure after τ = 4.5 × 10−3 τﬀ = 8000yr passed. In the model with α = 1, the outﬂow gas
ﬂows through a region whose shape resembles a capital letter U. This is similar to that observed inmodel (α = 1 and Ω = 5) shown in Figure 4.19. Decreasing the magnetic ﬁeld strength α = 1 → 0.1
and → 0.1, the shape of disk changes from a ﬂat disk to a round ellipsoid. The model with neither
magnetic ﬁeld nor rotation shows a spherical contracting core. Outﬂow is also aﬀected by decrease inmagnetic ﬁeld strength. Globally folded magnetic ﬁeld lines are seen in the outﬂow region in panel(e). Smaller-scale structure folding the magnetic ﬁeld lines develops in the outﬂow region in panel(f). The strucure seen in the polidal magnetic ﬁeld is well correlated to the distribution of toroidal
magnetic ﬁeld strength. That is, a packet of strong Bφ pinch the gas and poloidal magnetic ﬁeld lines
radially inwardly. This means that the toroidal magnetic ﬁeld is ampliﬁed by the rotation motion andbecomes predominant over the poloidal one in the case with meak initial magnetic ﬁeld. Thus, thewhoop stress, which is expressed by the ﬁrst term of the Lorentz force of equation(C.13) and comesfrom the tension of toroidal magnetic ﬁeld, pinches the 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 ﬁeld. In this case, the magnetic
pressure gradient plays an important role. Gas of outﬂow is ejected perpendicularly from the disk
and forms ﬂow similar to the capital letter I. In conclusion, U-type and I-type outﬂows are generated
for strong magnetic model and weak magnetic model, respectively.

4.11

Evolution in the H-R diagram

4.11.1Main Accretion Phase

Before Deuterium Ignition
Gas accreting onto the protostar generates energy as follows
Lacc =GM∗ Ṁ∗R∗(4.117)

where M∗ , R∗ and Ṁ∗ represent the mass and radius of protostar and the mass accretion rate.
The central temperature increases with time. Finally thermonuclear fusion reaction of Deuterium
2 H(p, γ)3 He begins. Before Deuterium burning begins, the protostar is radiative, that is, the energy
is transported radiatively.
When a fresh gas of Δm accrets, thermal energy of ΔU = GM∗ Δm/R∗ increases. Virial theorem(eq.[2.120]) requires the potential energy must decrease (increase in the absolute volume) at
GM∗ Δm
,(4.118)
ΔW = −2ΔU  −R∗ /2

CHAPTER 4. LOCAL STAR FORMATION PROCESS

114

(a)

(b)

(c)
P3

P2

P1

Figure 4.19: Continuation of the evolution shown in Fig.4.18. Snapshots at the time of 0.6069τﬀ =
1.066Myr just ∼ 500 yr after from Fig.4.18. Panel (a) plots the structure of L10 which is similar to
Fig.4.18 (c). Spherical region inside ∼ 1.2 × 10−3 H has been swept by the outﬂow. Panel (b) is asnapshot of L12 grid which has 4-times ﬁner spatial resolution than L10 (a). In panel (c), we plottedmagnetic force dominated region with ‘+’, centrifugal force dominated with ‘*’, and thermal pressure
force dominated with blank ’ ’.

4.11. EVOLUTION IN THE H-R DIAGRAM

115

(a)

(d)

(b)

(e)

(c)

(f)

Figure 4.20: Diﬀerence in models with the same rotation speed Ω = 2 × 10−14 rad s−1 but diﬀerent
magnetic ﬁeld strength ((a) and (d): Bc = 10μG, (b) and (e): 3μG, and (c) and (f): 1μG).

CHAPTER 4. LOCAL STAR FORMATION PROCESS

116

Figure 4.21: Speciﬁc entropy distribution is plotted against the accumulated mass (Stahler, Shu, &Taam 1982). The bottom curve shows the distribution at t = 2.3 × 104 yr (M∗ = 0.23M ). The top
curve is for t = 3.4 × 104 yr. The dashed curve indicates the entropy distribution if no convection is
generated. While, the solid curves represent that of new distribution achieved with convection.
to achive a mechanical equilibrium. As a result, the fresh gas must contract to about half the radiusat which the gas ﬁrst joined the core (Stahler, Shu, & Taam 1980).
When the free-falling fresh gas accrets on the static star, an accretion shock forms. Since the gas
temperature is increased with passing the shock front, temperature of the postshock gas, Tg , is muchhigher than that of the radiation, Tr . Thus, the postshock gas cools very eﬀectively. The postshock
region with Tg > Tr is called radiative relaxation region. The outgoing luminosity at the accretion
shock is much larger than that inside of the radiative relaxation region.
In the case of low mass stars since the Kelvin-Helmholtz contraction time tK−H  GM∗2 /R∗ L(L represent the luminosity at the base of the radiative relaxation region) is much longer than theacctretion time scale tacc  M∗ /Ṁ∗ . This gives Lacc  GM∗ Ṁ∗ /R∗ is much larger than L, which isconsistent with the above statement that in the radiative relaxation region a large amount of accretiontacc , the speciﬁc entropy inside the relaxation region isluminosity is radiated away. Since tK−H
essentially frozen to that when the gas obtained passing through the relaxtaion region. Figure 4.21
taken from Stahler, Shu & Taam (1982) shows the distribution of the speciﬁc entropy against theaccumulated mass Mr . The bottom curve corresponds to the state before the nuclear burning begins
when no entropy generation occurs. The temperature increases with mass, bacause the star must becompressed to support an extra mass. After the temperature becomes high enough for Deuterium
burning reaction 2 H(p, γ)3 He as T ∼ 106 K, an extra energy is liberated by the nuclear fusion reaction.
This increases the speciﬁc entropy mainly in the oﬀcenter region. Figure 4.21 clearly shows that shell
Deuterium burning occurs at M ∼ 0.025M .
After Deuterium Ignition
One gram of interstellar matter generates energy by the reaction of 2 H(p, γ)3 He asδ≡1
mp μ H

nDnH

D→He 

(4.119)

4.11. EVOLUTION IN THE H-R DIAGRAM

117

Figure 4.22: Luminosities in protostars accreting at Ṁ∗ = 10−5 M yr−1 (left) and Ṁ∗ = 10−4 M yr−1(right) (Palla & Staler 1993). Accretion luminosity Lacc , Deuterium burning luminosity LD , maximum
radiatively transferable luminosity Lrad and luminosity at the surface Lsurf are plotted.
where we used an energy per reaction of 2 H(p, γ)3 He as D→He = 5.5MeV and a composition of
Deuterium of nD /nH = 2.5 × 10−5 . This gives an estimate of the Deuterium-produced luminosity as
LD = Ṁ δ = 15LṀ
10−5 M yr−1,

(4.120)for a star with the accretion rate of Ṁ  10−5 M yr−1 .
Entropy inversion ∂s/∂r < 0 is unstable for the convection (see §2.7 in p.47). Convection is
generated in the region of ∂s/∂M < 0. Convective motion transfers the energy and the entropy. The
speciﬁc entropy peak generated by the Deuterium burning is smoothed out by the convection. The
convection induced with the Deuterium burning brings a fresh gas containing Deuterium into the hot
Deuterium burning region. Thus, although the consent of Deuterium is small and Deuterium burning
could not continue long, Deuterium burning continues actually as long as the convection can mix theouter Deuterium rich fresh gas into the burning region due to the convection.
In Figure 4.22, luminosities in protostars accreting at Ṁ∗ = 10−5 M yr−1 (left) and Ṁ∗ =−4
10 M yr−1 (right) are plotted (Palla & Stahler 1993). Lrad represents the maximum luminosity
transferable radiatively for a star with mass M∗ and radius R∗ as

Lrad ∼ 1L

M∗
1M

11/2 

R∗
R

−1/2

,

(4.121)

where the Krammers opacity κ ∝ ρT −3.5 is applied for stars with M <∼5M . Actual luminosity at thesurface of the star4,
(4.122)
Lsurf = 4πR∗2 Teﬀ

118

CHAPTER 4. LOCAL STAR FORMATION PROCESS

is larger than Lrad for low mass stars (M∗ < M1 = 2.4M for Ṁ∗ = 10−5 M yr−1 and M1 = 3.4M
for Ṁ∗ = 10−4 M yr−1 ). Excess energy is transferred by the convection, Lsurf = Lrad + Lcon . Thisﬁgure also shows that except for extremely low-mass object the accretion luminosity is the majorenergy source over Deuterium burning Lacc > LD . Since massive objects with M∗ > M3 satisﬁes
Lrad  Lsurf > Lacc , interior of such stars are fully radiative (energy is transferred by the radiation).
M3 = 3.9M for Ṁ∗ = 10−5 M yr−1 and M3 = 5.9M for Ṁ∗ = 10−5 M yr−1 .
Figure 4.23 (left) illustrates how the Deuterium burning in protostars evolves. Accumulating
infalling gas and increasing M∗ , the structure of protostars is changed. (a) in low-mass stars, convection feeds the Deuterium of the accreting matter to the burning center. (b) as increasing themass and thus the temperature, the opacity of the stellar interior decreases. Finally a radiative zone
appears and it prevents from mixing Deuterium. (c) as a consequence, central Deuterium is depleted
and energy generation rate is reduced. Thus the central region becomes radiative. (d) For higher
mass stars, Deuterium ignites as a shell source just outside the radiative central region. Due to thisDeuterium shell burning, which is similar to the Hydrogen shell burning in the red giants, radius of
the protostar expands. Figure 4.23 (right) plots the mass-radius relation of a protostar with accretion
rate of Ṁ∗ = 10−5 M yr−1 . The left open circle represents the beginning of central Deuterium burning (state a). On the second circle, radiative barrier appears (state b) and Deuterium shell burning
begins. Beyond the mass, the radius swells dramatically. Further increase in the mass increases thegravity and the radius shrinks. At the third open circle, Hydrogen begins to burn at the center ofthe protostar.
In the HR diagram, the mass-radius relation for the protostars shown in Figure 4.23(right) corresponds to the dotted line (Fig.4.24). Observed premain-sequence stars: T Tauri stars (M∗ <∼2M )<<and Herbig Ae/Be stars (2M ∼M∗ ∼5M ) distribute below the dotted line. Thus, this line is often
called the birth line for stars. In Figure 4.24, evolutionary tracks without mass accretion Ṁ = 0 are
plotted from the birth line. Modulation of the birth line is related to the swell and shrink of radius
shown in Figure 4.23. That is, increase of luminosity around M  1M corresponds to the ﬁrst swellafter Deuterium burning begins (between the ﬁrst and second open circles). Another increase aroundM  2.5M is related to the Deuterium shell burning (after the second open circle).

4.11.2Premain-sequence Evolution

Premain-sequence evolutionary path is understood as follows: Consider a situation that a star withM∗ is left when accretion stops. Figure 4.22 indicates that a less-massive star with M∗ < M1 isfully convective. Thus, this evolves along the convective Hayashi track (Hayashi 1961) down to themain-sequence. In contrast, a massive star with M∗ > M3 (Fig.4.22) is fully radiative and evolves
along the radiative Henyey track (Henyey, LeLevier, & Levee 1955) to the main-sequence. These
diﬀerences are clearly indicated in Figure 4.24. We can see veritical Hayashi track is seen essentiallyfor stars with M <∼2.5M evolve to the upper-left direction along the∼1.5M . Massive stars with M >Henyey track.

4.11. EVOLUTION IN THE H-R DIAGRAM

119

Figure 4.23: Schematic view of Deuterium burning in protostars (left). Evolution of the radius of a
protostar is plotted against its mass (right). Both ﬁgures are taken from Palla & Stahler (1990).

120

CHAPTER 4. LOCAL STAR FORMATION PROCESS

Figure 4.24: Premain-sequence evolutionary path in the H-R diagram. (Palla & Staler 1993). Evolution of stars with diﬀerent masses from 0.6M to 6M are calculated from the birth line (dottedline).

Bibliography[1] Alves, J., Lada, C.J., & Lada, E.A. 2001, Internal Structure of a Cold Dark Molecular Cloud
Inferred from Extinction of background Starlight, Nature, 409, 159[2] Alves, J., Lada, C.J., Lada, E.A., Kenyon, S.J., & Phelps, R. 1998, Dust Extinction and Molecular
Cloud Structure: L977, ApJ, 506, 292[3] André, P. 1994, Observations of Protostars and Protostellar Stages in “The Cold Universe,” eds.
T. Montmerle, C.J.Lada, I.F.Mirabel, & J.Tran Thanh Van, Editions Frontières, p.179[4] André, P.m Ward-Thompson, D. & Barsony, M. 2000, From Prestellar Cores to Protostars: The
Initial Conditions of Star Formation, Protostars and Planets IV, 59[5] Binney, J., & Tremaine, S. 1988, Galactic Dynamics, Princeton Univ. Pr. Chapter 6.
[6] Bonnor W. B. 1956, Boyle’s Law and Gravitational Instability, MNRAS, 116, 351[7] Blandford, R. D., & Payne, D. G. 1982, Hydromagnetic Flows from Accretion Discs and theProduction of Radio Jets, MNRAS, 199, 883[8] Chandler, C. J., & Sargent, A. I. 1993, The Small-Scale Structure and Kinematics of B335, ApJ,
414, L29[9] Chandrasekhar, S. 1939, An Introduction to the Study of Stellar Structure, (Univ. of Chicago Pr.;
Chicago)[10] Chandrasekhar, S. & Fermi, E. 1953, Magnetic Fields in Spiral Arms, ApJ, 118, 113[11] Cox, A.N. (editor) 2000, Allen’s Astrophyisical Quantities 4-th Edition, Springer Verlag[12] Cox, D.P. & Smith, B.W. 1974, Large-Scale Eﬀects of Supernova Remnants on the Galaxy:
Generation and maintenance of a Hot Network of Tunnels, ApJ, 189, L105
[13] 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.
[14] Crutcher, R. M., Nutter, D. J., Ward-Thompson, D. & Kirk, J. M. 2004, SCUBA Polarization
Measurements of the Magnetic Field Strengths in the L183, L1544, and L43 Prestellar Cores, ApJ,
600, 279[15] Cudworth, K. M., & Herbig, G. 1979, In Two Large-Proper-Motion Herbig-Haro Objects, AJ,
84, 548
121

122

BIBLIOGRAPHY[16] De Young, D. S. 1986, Mass Entrainment in Astrophysical Jets, ApJ, 307, 62[17] Draine, B.T. 2011, Physics of the Interstellar and Intergalactic Medium, Princeton U. Pr.
[18] Ebert R. 1955, er die Verdichtung von H I-Gebieten. Mit 5 Textabbildungen, Z. Astrophys., 37,
217[19] Elmegreen, B.G.; Seiden, P.E. & Elmegreen, D.M. 1989, Spiral Arm Amplitude Variations andPattern Speeds in the Grand Design Galaxies M51, M81, and M100, ApJ, 343, 602[20] Fiebig, D. & Guesten, R. 1989, Strong magnetic ﬁelds in interstellar H2O maser clumps, AAp,
214, 333[21] Fiege, J.D., & Pudritz, R.E. 2001, Helical Fields and Filamentary Molecular Clouds - I, MNRAS,
311, 85[22] Field, G.B., Goldsmith, D.W., & Habing, H.J. 1969, Cosmic-Ray Heating of the Interstellar Gas,
ApJ Letterm 155, L 149[23] Fridlund, C.V.M., & Liseau, R. 1998, Two Jets from the Protostellar System L1551 IRS5, ApJL,
499, 75[24] Girart, J. M., Beltrán, M. T., Zhang, Q., Rao, R. & Estalella, R., 2009, Magnetic Fields in theFormation of Massive Stars, Science, 324, 1408
[25] Girart, J. M., Rao, R. & Marrone, D. P. 2006, Magnetic Fields in the Formation of Sun-LikeStars, Science, 313, 812[26] Goldsmith, P. F., & Langer, W. D. 1978, Molecular Cooling and Thermal Balance of Dense
Interstellar Clouds, ApJ, 222, 881[27] Gordon, M.A., & Burton, W.B. 1976, Carbon Monoxide in the Galaxy. I - The Radial Distribution of CO, H2, and Nucleons, ApJ, 208, 346[28] Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, Disk Frequencies and Lifetimes in Young Clusters,
ApJ, 553, L153
[29] Hanawa, T. 1999, Star Formation, in “Active Universe” in Japanese (Shokabo, Tokyo) ed. by K.
Shibata, J. Fukue, R. Matsumoto, & S. Mineshige[30] Hartmann, L. 1998 (1st edition), 2009 (2nd edition) Accretion Processes in Star Formation,
Cambridge U.Pr.
[31] Hayashi, C. 1961, Stellar Evolution in Early Phases of Gravitational Contraction, PASJ, 13, 450[32] Hayashi, C. et al. 1978, Utyubutsurigaku(in Japanese), Iwanami Shoten[33] Heithausen, A., Bensch, F., Stutzki, J., Falgarone, E., & Panis, J.F. 1998, The IRAM Key
Project: Small-Scale Structure of Pre-Star Forming Regions Combined Mass Spectra and Scaling
Laws, AAp, 331, L65[34] Henriksen, R., Andre, P., & Bontemps, S. 1997, Time-Dependent Accretion and Ejection Implied
by Pre-Stellar Density Proﬁles, AAp, 323, 549

BIBLIOGRAPHY

123[35] Henyey, L. G., Lelevier, R. Levee, R. D. 1955, The Early Phases of Stellar Evolution, PASP, 67,154[36] Hirano, N., Kameya, O., Nakayama, M., & Takakubo, K. 1988, Bipolar outﬂow in B335, ApJL,
327, L69[37] Honma, M., Sofue, Y., & Arimoto, N. 1995, The Molecular Front in Galaxies. II. Galactic-scale
Gas Phase Transition of HI and H2 , AAp, 304, 1[38] Hunter, C. 1977, The Collapse of Unstable Isothermal Spheres, ApJ, 218, 834[39] Ikeda, N. & Kitamura, Y. 2009, A C18O Study of the Origin of the Power-Law Nature in theInitial Mass Function, ApJ, 705, L95[40] Ikeda, N., Sunada, K., & Kitamura, Y. 2007, A Survey of Dense Cores in the Orion A Cloud,
ApJ, 665, 1194[41] Itoh, Y. et al. 2000, A Pair of Twisted Jets of Ionized Iron from L 1551 IRS 5, PASJ, 52, 81[42] Johnstone, D., Wilson, C. D.; Moriarty-Schieven, G., Joncas, G., Smith, G., Gregersen, E., &Fich, M. 2000, Large-Area Mapping at 850 Microns. II. Analysis of the Clump Distribution in theOphiuchi Molecular Cloud, ApJ, 545, 327[43] Johnstone, D., Fich, M., Mitchell, G. F.& Moriarty-Schieven, G. 2001, Large Area Mapping at
850 Microns. III. Analysis of the Clump Distribution in the Orion B Molecular Cloud, ApJ, 559,
307[44] Kandori, R. et al. 2009, Distortion of magnetic ﬁelds in the pre-stellar core Barnard 68, Cosmic
Magnetic Fields: From Planets, to Stars and Galaxies Proceedings IAU Symposium No. 259, 2008
K.G. Strassmeier, A.G. Kosovichev & J.E. Beckman, eds.
[45] Kawamura, A. et al. 2009, The Second Survey of the Molecular Clouds in the Large Magellanic
Cloud by NANTEN. II. Star Formation, ApJS, 184, 1[46] Kennicutt, R.C. 1998, The Global Schmidt Law in Star-forming Galaxies, ApJ, 498, 541[47] Kennicutt, R.C., Tamblyn, P., & Congdon, C.W. 1994, Past and Future Star Formation in Disk
Galaxies, ApJ, 435, 22[48] Kim, S.; Staveley-Smith, L., Dopita, M. A., Sault, R. J., Freeman, K. C., Lee, Y., Chu, Y.-H.
2003, A Neutral Hydrogen Survey of the Large Magellanic Cloud: Aperture Synthesis and Multibeam
Data Combined, ApJS, 148, 473[49] Kramer, C., Stutzki, J., Rohrig, R., & Corneliussen, U. 1998, Clump Mass Spectra of Molecular
Clouds, AAp, 329, 249[50] Kroupa, P., Tout, C. A. & Gilmore, G. 1993, The distribution of low-mass stars in the Galactic
disc, MNRAS, 262, 545[51] Lada, E. A., Bally, J., & Stark, A. A. 1991, An unbiased survey for dense cores in the LYNDS
1630 molecular cloud, ApJ, 368, L432

124

BIBLIOGRAPHY[52] Lada, C.J. 1999, The Formation of Low Mass Stars, in “The Origin of Stars and Planetary
Systems” ed. C.J. Lada & N.D. Kylaﬁs, (Kluwer) p.143[53] Larson, R. B. 1981, Turbulence and Star Formation in Molecular Clouds, MNRAS, 194, 809[54] Launhardt, R., Pavlyuchenkov, Ya., Gueth, F., Chen, X., Dutrey, A., Guilloteau, S., Henning,Th., Pitu, V., Schreyer, K. & Semenov, D. 2009, Rotating molecular outﬂows: the young T Tauri
star in CB 26, AAp, 494, 147L
[55] Launhardt, R. & Sargent, A. I. 2001, A Young Protoplanetary Disk in the Bok Globule CB 26?,
ApJ, 562, L173
[56] Loren, R. B. 1989, The Cobwebs of Ophiuchus. I - Strands of 13 CO - The Mass Distribution,
ApJ, 338, 902[57] 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[58] Masunaga, H., Miyama, S. M., & Inutsuka, S. 1998, A Radiation Hydrodynamic Model for
Protostellar Collapse. I. The First Collapse, ApJ, 495, 346[59] Matsumoto, T., Hanawa, T., & Nakamura, F. 1997, Gravitational Contraction of Rotating
Clouds: Formation of Self-similarly Collapsing Disks, ApJ, 478, 569[60] Matsumoto, T., Nakamura, F., Hanawa, T. 1994, Fragmentation of a Magnetized Filamentary
Molecular Cloud Rotating around its Axis, PASJ, 46, 243[61] McKee, C.F. & Ostriker, J.P. 1977, A Teory of the Interstellar Medium – Three Components
Regulated by Supernova Explosions in an Inhomogeneous Substrate, ApJ, 218, 148[62] Meyer, M. R., Adams, F. C., Hillenbrand, L. A., Carpenter, J. M. & Larson R. B. 2000, The
Stellar Initial Mass Function: Constraunts from Young Clusters, and Theoretical Perspectives,Protostars and Planets IV, 121[63] Miller, G. E., Scalo, J. M. 1979, The initial mass function and stellar birthrate in the solar
neighborhood, ApJS, 41, 513[64] Mizuno, A., Onishi, T., Hayashi, M., Ohashi, N., Sunada, K., Hasegawa, T., & Fukui, Y. 1994,
Molecular Cloud Condensation as a Tracer of Low Mass Star Formation, Nature, 368, 719[65] Mizuno, A. et al. 1995, Overall Distribution of Dense Molecular Gas and Star Formation in the
Taurus Cloud Complex, ApJL, 445, L161
[66] Mouschovias, T. Ch. 1976a, Nonhomologous contraction and equilibria of self-gravitating, magnetic interstellar clouds embedded in an intercloud medium: Star formation. I Formulation of the
problem and method of solution, ApJ, 206, 753[67] Mouschovias, T. C. 1976b, Nonhomologous Contraction and Equilibria of Self-Gravitating, Magnetic Interstellar Clouds Embedded in an Intercloud Medium: Star Formation. II - Results, ApJ,
207, 141

BIBLIOGRAPHY

125[68] Mouschovias, T. Ch., Paleologou, E. V., & Fiedler, R. A. 1985, The Magnetic Flux Problem andAmbipolar Diﬀusion during Star Formation - One-Dimensional Collapse. II - Results, ApJ, 291,
772
[69] Motte, F., & André, P. 2001, The Circumstellar Environment of Low-Mass Protostars: A Millimeter Continuum Mapping Survey, AAp, 365, 440[70] Motte, F., André, P., Ward-Thompson, D.,& Bontemps, S. 2001, A SCUBA Survey of the NGC
2068/2071 Protoclusters, AAp, 372, L41[71] Myers, P.C. 1978, A Compilation of Interstellar Gas Properties, ApJ, 225, 380[72] Mundt, R., & Fried, J. W. 1983, Jets from Young Stars, ApJ, 274, 83[73] Nakano, T. 1979, Quasistatic Contraction of Magnetic Protostars due to Magnetic Flux Leakage- Part One - Formulation and an Example, PASJ, 31, 713[74] Nakano, T. Nishi, R., Umebayashi, T. 2002, Mechanism of Magnetic Flux Loss in Molecular
Clouds, ApJ, 572, 199[75] Ogino, S., Tomisaka, K., & Nakamura, F. 1999, Gravitational Collapse of Spherical Interstellar
Clouds PASJ, 51, 637.
[76] Ohashi, N., Hayashi, M., Ho, P.T.P., Momose, & M., Hirano, N. 1996, Possible Infall in the Gas
Disk around L1551 IRS 5, ApJ, 466, 957[77] Ohashi, N. Lee, S. W., Wilner, D. J. & Hayashi, M. 1999, CS Imaging of the Starless Core
L1544: an Envelope with Infall and Rotation, ApJL, 518, L410
[78] Onishi, T., Mizuno, A., Kawamura, A., Ogawa, H., Fukui, Y. 1996, A C18 O Survey of Dense
Cloud Cores in Taurus: Core Properties, ApJ, 468, 815[79] Onishi, T., Mizuno, A., Kawamura, A., Tachihara, K., & Fukui, Y. 2002, A Complete Search for
Dense Cloud Cores in Taurus, ApJ, 575, 950[80] 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[81] Palla, F. & Stahler, S. W. 1990, The Birthline for Intermediate-Mass Stars, ApJ, 360,47[82] Palla, F. & Stahler, S. W. 1993, The Pre-Main-Sequence Evolution of Intermediate-Mass Stars,
ApJ, 418, 414[83] Reipurth, B., Yu, K.C., Rodrı́guez, L.F., Heathcote, S., & Bally, J. 1999, Multiplicity of the HH
111 Jet Source: it Hubble Space Telescope NICMOS Images and VLA Maps, AAp, 352, L83[84] Saito, M., Sunada, K., Kawabe, R., Kitamura, Y., & Hirano, N. 1999, he Initial Conditions for
Formation of Low-Mass Stars: Kinematics and Density Structure of the Protostellar Envelope in
B335, ApJ, 518, 334[85] Salpeter, E. E. 1955, The Luminosity Function and Stellar Evolution, ApJ, 121, 161

126

BIBLIOGRAPHY[86] Sanders, D.B., Scoville, N.A., & Solomon, P.M. 1985, Giant Molecular Clouds in the Galaxy. II.
Characteristics of Discrete Features, ApJ, 289, 373[87] Scalo, J. M. 1986, The stellar initial mass function Fundamentals of Cosmic Physics, 11, 1[88] Schmidt, M. 1959, The Rate of Star Formation, ApJ, 129, 243[89] Shu, F. H. 1978, Self-similar Collapse of Isothermal Spheres and Star Formation, ApJ, 214, 488[90] Shu, F. H. 1983, Ambipolar Diﬀusion in Self-Gravitating Isothermal Layers, ApJ, 273, 202[91] Shu, F.H., Milione, V., & Roberts, W.W.,Jr. 1973, Nonlinear Gaseous Density Waves and Galactic Shocks, ApJ, 183, 819[92] Snell, R. L., Loren, R. B., & Plambeck, R. L. 1980, Observations of CO in L1551 - Evidence forStellar Wind Driven Shocks, ApJ, 239, 17[93] Spitzer, L. Jr. 1978, Physical Processes in the Interstellar Medium, Wiley, Cahper 13.
[94] Stahler, S. W. 1983, The Birthline for Low-mass Stars, ApJ, 274, 822[95] Stahler, S. W. 1994, The Kinematics of Molecular Outﬂows, ApJ, 422, 616[96] Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, The Evolution of Protostars. II - The HydrostaticCore, ApJ, 242, 226[97] Strong, A.W. et al. 1988, ICO to N (H2 ) Conversion, AAp, 207,1[98] Stutzki, J. & Guesten, R. 1990, High Spatial Resolution Isotopic CO and CS Observations of
M17 SW - The Clumpy Structure of the Molecular Cloud Core, ApJ, 356, 513[99] Sugitani, K. et al. 2011, Near-infrared-imaging Polarimetry Toward Serpens South: Revealing
the Importance of the Magnetic Field, ApJ, 734, 63[100] Tafalla, M., Mardones, D., Myers, P. C., Caselli, P., Bachiller, R. & Benson, P. J. 1998, L1544:
A Starless Dense Core with Extended Inward Motions, ApJ, 504, 900[101] Tamura, M., Hough, J.H., Hayashi, S.S., 1995, Millimeter Polarimetry of Young Stellar Objects:
Low-Mass Protostars and T Tauri Stars, ApJ, 448, 346[102] Testi, L. & Sargent, A. I. 1998, Star Formation in Clusters: A Survey of Compact MillimeterWave Sources in the Serpens Core, ApJ, 508, L91[103] Tomisaka, K. 2002, Collapse of Rotating Magnetized Molecular Cloud Cores and Mass Outﬂows,
ApJ, 575, 306[104] Tomisaka, K., Ikeuchi, S., & Nakamura, T. 1988, The Equilibria and Evolutions of Magnetized,
Rotating, Isothermal Clouds. I - Basic Equations and Numerical Methods, ApJ, 326, 208[105] Tomisaka, K., Ikeuchi, S., & Nakamura, T. 1988, Equilibria and Evolutions of Magnetized,
Rotating, Isothermal Clouds. II - The Extreme Case: Nonrotating Cloud, ApJ, 335, 293[106] Tomisaka, K., Ikeuchi, S., & Nakamura, T. The Equilibria and Evolutions of Magnetized, Rotating, Isothermal Clouds. IV - Quasi-Static Evolution, ApJ, 362, 202

BIBLIOGRAPHY

127[107] Toomre, A. 1981, What Ampliﬁes the Spirals, in “The Structure and Evolution of Normal
Galaxies” Cambridge U.Pr., 1981, pp.111-136[108] Tothill, N. F. H., White, G. J., Matthews, H. E., McCutcheon, W. H., McCaughrean, M. J.,& Kenworthy, M. A. 2002, The Structure and Evolution of the Lagoon Nebula. I. Submillimeter
Continuum and CO Line Mapping, ApJ, 580, 285[109] 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[110] Ward-Thompson, D., Kirk, J.M., Crutcher, R.M., Greaves, J.S., Holland, W.S., André, P. 2000,
First Observations of the Magnetic Field Geometry in Prestellar Cores, ApJL, 537, L135
[111] Weintraub, D.A., Goodman, A.A., & Akeson, R.L. 1999, Polarized Light from Star-Forming
Regions, in Protostars and Planets IV, ed. V.Mannings, A.P.Boss, & S.S. Russell, U. Arizona Pr.
pp.247-271[112] Whitworth, A., & Summers, D. 1985, Self-Similar Condensation of Spherically Symmetric SelfGravitating Isothermal Gas Clouds, MNRAS, 214, 1[113] Wolf, S., Launhardt, R. & Henning, T., 2003, Magnetic Field Evolution in Bok Globules, ApJ,
592, 233[114] Yun, J. L., & Clemens, D. P. 1994, Outﬂows from Young Stellar Objects in BOK Globules:
Maps, ApJS, 92, 145[115] Zhou, S. 1995, Line Formation in Collapsing Cloud Cores with Rotation and Applications to
B335 and IRAS 16293-2422, ApJ, 442, 685[116] Zhou, S., Evans, N.J.,II, Koempe, C., & Walmsley, C.M. 1993, Evidence for Protostellar Collapse in B335, ApJ, 404, 232

128

BIBLIOGRAPHY

Appendix A

Basic Astrophysical Quantities
M
L
R
mp
me
AU
pc
amu
G
k
h
c
σ

the Solar mass
the Solar luminosity
the Solar radius
proton masselectron mass
the astronomical unit
parsecatomic mass unit
gravitational constantBoltzmann constant
Plank constantlight speed
Stefan-Bolzmann constant

1.989 × 1033 g
3.83 × 1033 erg s−1
6.96 × 1010 cm
1.67262 × 10−24 g
9.10938 × 10−28 g
1.49598 × 1013 cm
3.0857 × 1018 cm
1.66054 × 10−24 g
6.67 × 10−8 dyn cm2 g−2
1.38065 × 10−16 erg K−1
6.62607 × 10−27 erg s
2.99792458 × 1010 cm s−1
5.67040 × 10−5 erg s−1 cm−3 K−4Table A.1: Astrophysical quantities, units, and constants.

129

130

APPENDIX A. BASIC ASTROPHYSICAL QUANTITIESAppendix B

Basic Equation of Fluid Dynamics
B.1

What is ﬂuid?

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

B.2

Equation of Motion

Consider a ﬂuid element with mass m as in B.1. From the Newton’s second law of mechanics, theacceleration of an element is written as
Fdv= 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 ﬂuid dynamics, using the mass density ρ and the force working on the unit volume f v equation(B.1) is rewritten asdv= f v.
(B.2)
ρ
dt
Which kind of force works in a ﬂuid? Gas pressure force does work in any ﬂuids. Beside this, if
there is the gravity, ρg should be included in f v . If the electric currents is running in the ﬂuid andthe magnetic ﬁelds exist, the Lorentz force j × B should be added.
To write down the expression of the gas pressure force, consider a ﬂuid element between x andx + Δx as shown in Figure B.2. Pressure force exerting on the surface S at x is p(x)S, while that

m
V
ρ

m
V’
ρ’

t

t+Δt

Figure B.1: Movement of ﬂuid element
131

APPENDIX B. BASIC EQUATION OF FLUID DYNAMICS

132

S
p(x)

p(x+Δx)x

x+Δx

Figure B.2: Pressure force exerted on an element with volume SΔxon the opposite side is −p(x + Δx)S. The net pressure force working on the volume of SΔx is equalto (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 Euler 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 diﬀerent.
However, considering the velocity ﬁeld in the space, the time derivative of the velocity should be
calculated staying at a ﬁxed point x.
These two time derivative are diﬀerent 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 diﬀerence dF
dt Δt, using the Lagrangian time derivative of F , represents the the diﬀerenceof F (t + Δt) from F (t) focusing on a speciﬁc ﬂuid element, whose positions are diﬀerent 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 diﬀerence is expressed asdF
Δt = F (x0 + v 0 Δt, t0 + Δt) − F (x0 , t0 ),
dt



 
∂F 
∂F 
· v0 Δt +Δt,
∂x t x0 ,t0∂t x x0 ,t0(B.4)

where we used the Taylor expansion of F . The diﬀerence 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, theB.4. CONTINUITY EQUATION

133

Lagrangian 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:


ρ

B.4



∂v
+ v · ∇v = −∇p + ρg.
∂t

(B.7)

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 ﬂuid element whose volume is equal
to ΔV . The mass contained in the volume is constant. Thus
dρ
dΔV
dρΔV
=
ΔV +
ρ = 0.
dt
dt
dt

(B.8)

The variation of the volume dΔV
dt is rewritten as
dΔV
=
dt


∂ΔV

v · dS =

divvdV,(B.9)

ΔV

where ∂ΔV represents the surface of the ﬂuid 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 Euler derivative are equations (B.7) and (B.11).

B.4.1

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 cordinate (x, y, z) = (xi , x2 , x3 ), thisbecomes
⎛

⎞

⎛

⎞⎛

⎞

⎛

⎞⎛

⎞

⎛

⎞

ρv
∂/∂x
ρvx vx ρvx vy ρvx vz
∂/∂x
p 0 0
gx
∂ ⎜ x ⎟ ⎜
⎟⎜
⎟ ⎜
⎟⎜
⎟
⎜
⎟
⎝ ρvy ⎠ + ⎝ ∂/∂y ⎠ ⎝ ρvy vx ρvy vy ρvy vz ⎠ + ⎝ ∂/∂y ⎠ ⎝ 0 p 0 ⎠ = ρ ⎝ gy ⎠ .
∂tρvz
ρvz vx ρvz vy ρvz vz
gz
∂/∂z
∂/∂z
0 0 p
(B.13)

APPENDIX B. BASIC EQUATION OF FLUID DYNAMICS

134
Using the momentum stress tensor

πi j = ρvi vj + δi j p,

(B.14)∂ρvi ∂πi j= ρgi ,+
∂t∂xj

(B.15)

the above equation is written as

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 ﬂuid and thus the momentumof the ﬂuid element is not conserved.

B.5

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.5.1

Polytropic Relation

If the pressure of the ﬂuid, 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 suﬃcient to describe
the dynamics of the ﬂuid. Occasionally, the presure 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 ﬂuid, in which the entropy is kept constant. In this case thepressure is proportional to ργ , as(B.20)
p = c2s ργ ,where γ is the speciﬁc heat ratio = cp /cv and cs = (γkT /μmp )1/2 represents the adiabatic soundspeed. In these cases the polytropic replations of equations (B.19) and (B.20) plays a role as thethird equation of basic equations of hydrodynamics.

B.6. SHOCK WAVE

B.5.2

135

Energy Equation from the First Law of Theromodynamics

In more general cases, the last equation comes from the ﬁrst 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 or subtracted
from the system. Using the total energy per unit volume
U1= ρv 2 +2V(B.22)

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

(B.23)

Since the total energy per unit volume is expressed using the basic physical quantities as11p,= ρv 2 +2γ−1(B.24)

equation (B.23) is the ﬁnal basic equation for hydrodynamics. Equations (B.7), (B.11), and (B.23)
are basic equations hydrodynamics using the Eulerian time derivative.
Problem: 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 enerygy 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).

B.6

Shock Wave

B.6.1

Rankine-Hugoniot Relation

Passing through a shock front moving with a speed Vs , the physical variables ρ, p, and u change
abruptly. Since the basic equations of hydrodynamics is unchanged after chosing a system moving
Vs , the continuity equation∂ρ ∂ρu+=0(B.26)∂t∂x
gives an equation for a steady state as∂ρu= 0.
(B.27)∂x
Considering a region containing the shock front at xs which extends from x = xs − Δx to x = xs + Δx
and integrating the above equation, we get xs +Δx
xs −Δxdρudx = [ρu]xxss +Δx
−Δx .
dx

(B.28)

Thus, we obtain the jump condition conserning the mass coservation as
ρ1 u1 = ρ2 u2 ,(B.29)

136

APPENDIX B. BASIC EQUATION OF FLUID DYNAMICS

where the quatities with suﬃx 1 are for preshock and those with suﬃx 2 are for postshock.
Equation of motion for steady state
∂p
∂u
=− ,
∂x
∂x

(B.30)

p1 + ρ1 u21 = p2 + ρ2 u22 ,(B.31)

ρu
gives
where we used equation(B.29).
Isotharmal shock

In the case of the gas is isothermal p = c2is ρ, equation (B.31) becomes

ρ1 c2is + u21 = ρ2 c2is + u22 .

(B.32)

Eliminating ρ from equations (B.29) and (B.32), we obtainu1 u2 − c2s (u1 − u2 ) = 0,

(B.33)

u1 u2 = c2s .

(B.34)u1
u2
ρ2== 21 .
ρ1
u2
cs(B.35)

which means
From equation (B.29),

This indicates the postshock velocity u1
density becomes large.

cs the ratio of the postshock density to the preshockAppendix C

Basic Equations ofMagnetohydrodynamics
C.1

Magnetohydrodynamics

Here, we derive the basic equations of magnetohydrodynamics. Diﬀerences from hydrodynamics are
the Lorentz force in the equation of motion and the induction equation of magnetic ﬁeld. The Lorentz
force is written as
1(∇ × B) × B,
(C.1)
j×B =
4π
where we used the Ampere’s law ∇ × B = (4π/c)j and the displacement current is ignored. Thus,
the equatin of motion becomes as




1
∂v
+ (∇ · v)v = −∇p − ρ∇ψ +(∇ × B) × B,
ρ
∂t
4π

(C.2)

in the Eulerian form.
The Faraday’s law ∇ × E = −(1/c)(∂B/∂t) and the Ohm’s law j = σ(E + v × B/c) (σ is theelectric conductivity) lead the induction equation of the magnetic ﬁeld as
c2
∇×B
∂B
= ∇ × (v × B) −
∇×
.
∂t
4π
σ

(C.3)

The factor c2 /4πσ is the coeﬃcient of magnetic diﬀusibility. When we assume the electric conductivity
σ → ∞, diﬀusive MHD equation, equation(C.3), reduces to the ideal MHD equation as∂B= ∇ × (v × B).
∂t

C.1.1

(C.4)

Flux Freezing

An important nature of the ideal MHD is derived from equatio (C.4), that is, the matter and themagnetic ﬁeld is coupled with each other. Consider a gas parcel threaded with the magnetic ﬁeld B.
The magnetic ﬂux anchored to the gas parcel ΦB ≡ B · dS is changes
d
dΦB
=
dt
dt



B · dS =∂B· dS +∂t

B · v×ds,(C.5)

where s = ∂S. The ﬁrst term of the lhs of the equation comes from a time derivative of the magnetic
ﬂux density, while the last term represents the change of the integral region due to gas motion (see137

APPENDIX C. BASIC EQUATIONS OF MAGNETOHYDRODYNAMICS

138

Figure C.1: The change of ﬂux anchored to a surface dS. One part of the shange comes from thetemporal change in magnetic ﬂux density. The other comes from the change of the boundary of thesurface due to gas motion. Time variation of the area is given by vdt × ds.
Fig.C.1).
(dS = v × ds) Using the Stokes’ theorem ∇ × A · dS =to B × v · ds = ∇(B × v) · dS. Finally, we obtain
dΦB
=
dt

 



A · ds, the last term is reduced∂B− ∇(v × B) · dS = 0,
∂t

(C.6)

where we used equation (C.4). This means that the magnetic ﬂux anchored to a gas parcel does notchange in the ideal MHD regime.

C.1.2

Basic Equations of Ideal MHD

Basic equations are as follows: The mass continuity as∂ρ+ ∇ · (ρv) = 0,∂t
the equation of motion as


ρ

(C.7)



1
∂v
+ (∇ · v)v = −∇p − ρ∇ψ +(∇ × B) × B,
∂t
4π

(C.8)

the equation of thermal energy as∂+ div( + p)v = ρv · g,∂tor some barotropic relation p = P (ρ) and the induction equation as∂B= ∇ × (v × B).
∂t

C.1.3

(C.9)

(C.10)

Axisymmetric Case

The basic equations to be solved are the magnetohydrodynamical equations and the Poisson equation
for the gravitational potential. In cylindrical coordinates (z, r, φ) with ∂/∂φ = 0, the equations are
expressed as follows:
∂
1 ∂
∂ρ
+
(ρvz ) +(rρvr ) = 0,
(C.11)
∂t
∂z
r ∂r
∂ρvz
∂t

∂
1 ∂
(ρvz vz ) +(rρvz vr ) =∂zr ∂r 


∂ψ
1
∂Bz
∂Bφ
∂ρ
∂Br
−ρ
+
Bφ −
−
−
Br ,
− c2s
∂z
∂z
4π
∂z
∂z
∂r
+

(C.12)

C.1. MAGNETOHYDRODYNAMICS
∂ρvr
∂t

139

∂
1 ∂
(ρvr vz ) +(rρvr vr ) =∂zr ∂r 


1 ∂
∂ψ
1
∂Bz
∂ρ
∂Br
−
Bz ,
−ρ
+
(rBφ )Bφ +
−
− c2s
∂r
∂r
4π
r ∂r
∂z
∂r
+

(C.13)



∂Bφ
∂ρrvφ
1 ∂
r 1 ∂
∂
+
(ρrvφ vz ) +(rρrvφ vr ) =(rBφ )Br +
Bz
∂t
∂z
r ∂r
4π r ∂r
∂z
∂Bz
1 ∂
=
[r(vz Br − vr Bz )],
∂t
r ∂r
∂
∂Br
= − (vz Br − vr Bz ),
∂t
∂z
∂Bφ
∂
∂
=
(vφ Bz − vz Bφ ) −(vr Bφ − vφ Br ),
∂t
∂z
∂r


∂ψ
∂2ψ 1 ∂+
r
= 4πGρ,
∂z 2
r ∂r
∂r

(C.14)(C.15)(C.16)(C.17)(C.18)

where the variables have their ordinary meanings. Equation (C.11) is the continuity equation; equations (C.12), (C.13) and (C.14) are the equations of motion. The induction equations for the poloidal
magnetic ﬁelds are equations (C.15) and (C.16) and for the toroidal magnetic ﬁeld is equation (C.17).
The last equation (C.18) is the Poisson equation.

140

APPENDIX C. BASIC EQUATIONS OF MAGNETOHYDRODYNAMICS

Appendix D

Hydrostatic Equilibrium
D.1

Polytrope

If we choose the polytropic equation of state1 ,p = Kρ1+1/n ,(D.1)

the hydrostatic balance is expressed by1 ∂r 2 ∂rr 2 dpρ dr= −4πGρ,(D.2)

where we used equations (4.1) and (4.2). A hydrostatic gaeous star composed with a polytropic gas
is called polytrope. Normalizing the density, pressure and radius asρ ≡ ρc θ,
p ≡

(D.3)

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

1/2

ξ

(D.5)

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



ξ2

dθ
dξ



= −θ n ,(D.6)

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

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

141

(D.7)
(D.8)

APPENDIX D. HYDROSTATIC EQUILIBRIUM

142
1
0.8
0.6
0.4
0.2

1

2

4

3

5

7

6

Figure D.1: Solutions of eq.(D.6) are plotted for indices of n = 3 or Γ = 4/3 (the most extendedone), n = 5/2 or Γ = 7/5 (the middle), and n = 3/2 or Γ = 5/3(the most compact one). For n ≥ 5
gas extends inﬁnitely and the solution has no zero-point.
Mass of the polytrope is written down as
 r

M

=

4πρr 2 dr,0= 4πH 3 ξ

θ n ξ 2 dξ,0 ξd 2 dθ
ξ
dξ,
dξ
0 dξ

dθ 
= −4πH 3 ξ 2dξ ξ1= 4πH 3



= −4π



(n + 1)K 3/2 (3−n)/2n 2 dθ 
ρc
ξ
,
4πG
dξ ξ1(D.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 (D.9) reduces toK
M = −4π
πG

3/2



ξ
dξ 
2 dθ 

.

(D.10)

ξ1

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

7K
M = −4π
8πG

3/2




dθ 
ρ1/10
ξ2
,
c
dξ ξ1(D.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 as5K
M = −4π
8πG
where ξ1 = 3.65375 and −ξ12 dθ/dξ|ξ1 = 2.71406.

3/2




2 dθ 
ρ1/2
c ξ


dξ ξ1

,

(D.12)

D.2. MAGNETOHYDROSTATIC CONFIGURATION

D.2

143

Magnetohydrostatic Conﬁguration

In section 4.2, we obtained the maximum mass which is supported against the self-gravity using thevirial analysis. In this section, we will survey result of more realistic calculation.
Formalism was obtained by Mouschovias (1976a,6), which was extended by Tomisaka, Ikeuchi, &Nakamura (1988) to include the eﬀect of rotation. Magnetohydrostatic equlibrium is achived on abalance between the Lorentz force, gravity, thermal pressure force, and the centrifugal force as
1(∇ × B) × B − ρ∇ψ − ∇p − ρω × (ω × r) = 0.
4π

(D.13)In the axisymmetric case, the poloidal magnetic ﬁelds is obtained by the magnetic ﬂux function, Φ,or the φ-component of the vector potential as
1 ∂Φ
,
r ∂r
1 ∂Φ
,
= −
r ∂z
bφ
.
=
r

Bz =
Br
Bφ

(D.14)

Equation (D.13) leads tobφ ∂bφ
∂Φ
1
−
Δ1 Φ
4πr 2
∂z
4πr 2 ∂zbφ ∂bφ
∂Φ
1
−
Δ1 Φ
−
2
4πr
∂r
4πr 2 ∂r∂bφ ∂Φ ∂bφ ∂Φ
+
−
∂r ∂z
∂z ∂r
−

∂ψ ∂p
+
,
∂z
∂z
∂ψ ∂p
= ρ
+
− ρω 2 r,
∂r
∂r
= ρ

(D.15)(D.16)= 0,with(D.17)

∂2Φ
∂ 1 ∂Φ
+r.
(D.18)
Δ1 Φ ≡
2∂z∂r r ∂r
Equation (D.17) indicates bφ is a function of Φ as bφ (Φ), which is constant along one magnetic ﬁeld
line. Ferraro’s isorotation law demands, that is, to satisfy the stead-state induction equation ω isconstant along a magnetic ﬁeld. This means ω is also constant along one magnetic ﬁeld line, ω(Φ).
From this, the density distribution in one ﬂux tube is written






1
q
ρ = 2 exp − ψ − r 2 ω 2 /c2s .
cs2(D.19)

This means q is also constant along one magnetic ﬁeld line, q(Φ) . Since the forces are expressed by
the defrivative of function q






1
(∇q) exp − ψ − r 2 ω 2 /c2s = ∇p + ρ∇ψ − ρrω 2 − ρr 2 ω∇ω,2

(D.20)

where equation (D.19) is used, equation (D.15) and (D.16) are rewritten as∂Φ ∂(b2φ /2)
+
Δ1 Φ
∂z
∂z
∂Φ ∂(b2φ /2)
+
Δ1 Φ
∂r
∂r











1
∂ω
∂q
= −4πrexp − ψ − r 2 ω 2 /c2s + ρr 2 ω
,
∂z
2
∂z

 


1
∂ω
∂q
= −4πr 2
exp − ψ − r 2 ω 2 /c2s + ρr 2 ω
.
∂r
2
∂r
2

(D.21)(D.22)

APPENDIX D. HYDROSTATIC EQUILIBRIUM

144

Finally, using the fact that b, q, and ω are functions of Φ, these two equations are reduced to

 


d(b2φ /2)
1 2 2
dq
2
2
2 dω
− 4πr
exp − ψ − r ω /cs + ρr ω.
Δ1 Φ = −
dΦ
dΦ
2
dΦ

(D.23)

Another equation to be coupled is the Poisson equation as






1
q
Δψ = 4πG 2 exp − ψ − r 2 ω 2 /c2s .
cs2(D.24)The source terms of equations (D.23) and (D.24) are given by determining the mass δm(Φ) andthe angular momentum δL(Φ) contained in a ﬂux tube Φ–Φ + δΦ. Mass and angular momentum
distribution of zs (Φ)

 r(z,Φ+δΦ)dz

δm(Φ) = 20

δL(Φ) = 2

 zs (Φ)
0r(z,Φ)
 r(z,Φ+δΦ)dz

dr2πrρ(r, z),(D.25)

dr2πrρ(r, z)r 2 ω(Φ)(D.26)r(z,Φ)

is chosen artitrary in nature, where zs (Φ) is the the height of the cloud surface where the magnetic
potential is equal to Φ. For example, m(Φ) and L(Φ) are chosen as a uniformly rorating uniform r(Φ+δΦ)
2πrdr = δΦπ∂r 2 /∂Φdensity spherical cloud threaded by uniform magnetic ﬁeld. Since r(Φ)
 zs (Φ)

 1







1
dz2π ∂r /∂Φ 2 exp − ψ − r 2 ω 2 /c2s ,q(Φ) = (dm/dΦ)/
c
2
0
s
 

 zs (Φ)
 1
1 2 2
2
2
dz2π ∂r /∂Φ 2 q(Φ)r exp − ψ − r ω /c2s ,ω(Φ) = (dL/dΦ)/
cs
2
0
2

(D.27)(D.28)The source terms of PDEs [eqs (D.23) and (D.24)] are given from equations (D.27) and (D.28). While
the functons q(Φ) and ω(Φ) are determined from the solution of these PDEs after m(Φ) and L(Φ)
are chosen. This can be solved by a self-consistent ﬁeld method.

Appendix E

Basic Equations for Radiative
Hydrodynamics
E.1

Radiative Hydrodynamics

The basic equation for gas which is aﬀected by the radiation is as follows:dρ+ ρ∇ · v = 0,
dt
χF 0
dv
= ρ∇φ − ∇ψ +
F 0,
ρ
dt
c
d
ρ + p∇ · v = hCR ρ + cχE0 E0 − 4πχP 0 B,
dt
∇2 ψ = 4πGρ,p,=(γ − 1)ρ

(E.1)
(E.2)
(E.3)
(E.4)
(E.5)

where χ, hCR , and B are the absorption coeﬃcient, the cosmic-ray heating rate perunit mass, andteh Planck function B = σT 4 /π. In equation (E.2), the term χFc 0 F 0 represents the acceleration ofgas due to the photon pressure. In equation (E.3), hCR ρ, cχE0 E0 and 4πχP 0 B represent, respectively
the heating due to the CR particles, heating due to the absorption of radiation and cooling die to theemission. The frequence-integrated radiation energy density and radiation ﬂux are deﬁned as
 ∞

F0 =

dΩI(ν, n),

(E.6)

dν0

dΩI(ν, n)n,

(E.7)

0
 ∞
0



dν0

E0 =



where I(ν, n) denotes the speciﬁc intensity of radiation at frequency ν along the direction vector
n. These equation could be solved if the the radiation transfer is solved. The frequency-averaged
absorption coeﬃcients are deﬁned as follows:χF 0 =dνχ(ν)F 0 (ν) /F 0 ,χE 0 =
χP 0 =(E.8)

dνχ(ν)E0 (ν) /E0 ,


(E.9)



dνχ(ν)B(ν) /B,
145(E.10)

146

APPENDIX E. BASIC EQUATIONS FOR RADIATIVE HYDRODYNAMICS
While the equation for the radiation transfer is basically as follows:
1 ∂I(ν, r, n)+ n · ∇I(ν, r, n) = −χ(ν, r) [I(ν, r, n) − S(ν, r)] ,c∂twhere S(ν, r) represent the source term.

(E.11)

Appendix F

Random Velocity
Considering gas in Maxwellian velocity distribution, the distribution function for the velocity is asfollows:v2(F.1)
f (vx ) = A exp − x .
σThis gives the one-dimensional random velocity as

∞ 2
vx2
v
exp
−
x
0
σ2 dvx,∞vx2< vx2 > =

0
2
σ

=

2

exp − σ2 dvx

.

(F.2)

If we observe emissions from such a gas, the emission line is broaden due to the Doppler shift. Using
equation (F.2), the HWHM (half width of half maximum: the line width measured from the thecenter of the emission line to the point of the half intensity; see Fig.F.1) of the emission line is




2 

1
,
2

(F.3)

vx,HWHM = (2 ln 2 < vx2 >)1/2 ,(F.4)

vx,FWHM = (23 ln 2 < vx2 >)1/2 .

(F.5)

vx,HWHMexp −2 < vx2 >=which leads toandThus, if we assume isotropic distribution, three-dimensional random velocity of gas
2>=< vx2 > + < vy2 > + < vz2 >= 3 < vx2 >
< v3D

(F.6)

is obtained with the line width as




1/2
3
vx,FWHM ,23 ln 2
1/2
3vx,HWHM .
=2 ln 22>1/2 =
< v3D

147

(F.7)

APPENDIX F. RANDOM VELOCITY

148

Center

FWHM

Half Maxmum

HWHM

V

Figure F.1: HWHM (half width of half maximum: the line width measured from the the center ofthe emission line to the point of the half intensity) and FWHM (full width of half maximum: the line
width measured between the points of the half intensity).

Appendix G

X-factor
Since the H2 molecule is the homonuclear diatomic molecule, the H2 molecule has no electric dipole.
Thus, the electric dipole radiation is not expected for this molecule. While the second abundant
molecule CO is the heteronuclear diatomic molecules and has electric dipole. Lower levels of rotationaltransition of CO are relatively easily excited even in the cold interstellar medium of T ∼ 10K.
Therefore, the rotational transition of CO molecule is the ﬁrst possibility to observe cold moleculargas. The X-factor represents the ratio of the column density of H2 molecules N (H2 ) to the integrated
intensity of 12 CO line ICO as
NH2 (cm−2 )
.
(G.1)
X≡
ICO (K km s−1 )There are several empirical but physical estimations of the X-factor: Using γ-rays, X-factor is estimated as X  2.3 × 1020 (Strong et al. 1988). This is based on the idea that the emissivity of
γ-ray emission is proportional to the cosmic-ray intensity times the number of target nuclei. The
γ-ray intensity is proportional to the column density if the cosmic-ray intensity is uniform. From the
distribution of γ-ray intensity, Strong et al. (1988) estimated the total column density NH2 (cm−2 )using some model. This gives the value of X-factor.
Estimations using the Virial mass of the interstellar cloud and the visual extinction AV .

149

150

APPENDIX G. X-FACTOR

Appendix H

Abundance
H.1

Solar system abundance

Table H.1: Solar system abundance taken from meteorites and the Solar photosphere. Taken from
Table 3.2 of Allen’s Tables (Cox 2000).

151APPENDIX H. ABUNDANCE

152

H.2

Isotope ratio (Isotopic abundance)

Table H.2: Isotopic abundances for H, C, N, and O for various places. Taken from Table 21.10 ofAllen’s Tables (Cox 2000).

H.3

Molecular Weight

Table H.3: Group abundance ratios. Taken from Table 3.1 of Allen’s Tables (Cox 2000).