Thesis

A Radiation Hydrodynamical
Model
for Protostar Formation

Hirohiko Masunaga
Department of Astronomy
University of Tokyo
March 1999

Acknowledgment
I greatly appreciate valuable comments and continuous encouragement of Prof. Shoken
M. Miyama, the supervisor, and of Prof. Kazunari Shibata, the co-supervisor. I am very
grateful to Dr. Shu-ichiro Inutsuka for fruitful discussions and comments on scienticmatters, and for helpful advice on technical problems in numerical calculations. I would
like to express my thanks to Prof. Yasuo Fukui, Dr. Toshikazu Onishi, and Dr. Munetake
Momose for their helpful comments and advice. I also thank Prof. Harold W. Yorke for hiskindly sending me tabulated data of the dust opacities. I am grateful to Prof. Takenori
Nakano and Prof. Masahiko Hayashi for their help and encouragement as well as fordiscussions and comments. My research is supported by the Research Fellowships of theJapan Society for the Promotion of Science for Young Scientists.
Finally, I thank ism(Pentium II), yso(Solaris), luna(HP 9000/735), and bd(HP9000/782) for their excellent work for my jobs with little (but not without any) complaint.

1

Abstract

Investigation of star formation processes has been a signicant subject both for theoretical astrophysics and for observational astronomy. Nevertheless, many problems areleft unresolved still to date especially in an early phase in the evolution where a protostar is under dynamical growth. The lack of our knowledge for protostar formation is a
\missing link" in the evolutionary picture of star formation.
In order to construct a theoretical model for protostar formation which is capable
to account for recent observations, we have developed a numerical code for radiationhydrodynamic calculations which includes all the necessary physical processes, such as
dissociation of molecules, ionization of atoms, non-ideal equation of state, adequate opacities for gas and dust, and so forth. The numerical code is designed to yield the evolutionof spectral energy distributions as well as the dynamical evolution. We established in thiswork a theoretical model for the whole evolution of protostar formation in a consistentscheme which is accountable for observations. No such studies have been done to date.
First, we perform numerical calculations with simplied input physics for clarifying
the evolution before the central protostar is formed. Evolutions before hydrogen molecules
begin to dissociate (so-called \the rst collapse") are examined for dierent masses andinitial temperatures of the parent cloud cores and for dierent opacities. Numerical results
for a typical case (Tinit = 10K and P (10K ) 0:01cm2  g;1) show that the radiusand mass of the rst core are 5AU and 0:05M, respectively. These values are
independent of the mass of the parent cloud core. The analytical expressions for theradius, mass, and accretion luminosity of the rst core are also obtained. The SED
contains only cold components of a few10K throughout the rst collapse phase because
the opaque envelope veils the rst core from observers. We suggest that the molecular
cloud cores with luminosities higher than 0:1L should contain young protostars deepin the center even if they show no evidence for the existence of central stars in near-infraredand optical observations.
Following the study on the rst collapse, numerical calculations are carried out toreveal physical processes in the formation of a 1M protostar. The whole evolution are
pursued from the beginning of the rst collapse to the end of the main accretion phase.
We conrm that the typical features in the evolution are in good agreement with previous
studies.
We consider two dierent initial conditions for the density distribution: homogeneousand hydrostatic cloud cores. For the initially homogeneous model, the accretion luminosity
rapidly rises to the maximum value of 25L just after the formation of a protostar,and declines gradually as the mass accretion rate decreases. In contrast, the luminosity2

increases monotonically with time for the initially hydrostatic model. This discrepancy
arises because the mass accretion rate varies depending on the inward acceleration in theinitial condition, which aects the luminosity curve.
We conrm that the SED evolves from a 10K greybody spectrum to hotter spectra
typical for class I and II sources. The SED for the class 0 sources corresponds to the ageof 2  104yr, which is smaller by an order of magnitude than the typical age of class I
objects. Considering possible non-spherically-symmetric eects, we suggest that observed
class 0 sources should be the compound of the \genuine" class 0 that is as young as 104
yr and more evolved protostars on edge-on view. The contamination of older protostars
are not negligible because they are intrinsically abundant than genuine class 0 objects.
Since observations indicate that the class 0 sources are typically more luminous thanclass I sources, we exclude the initially hydrostatic model where the luminosity increases
monotonically with time. The initially homogeneous model, in contrast, is found to show
the tendencies consistent with observations.
We have also developed a numerical code for non-LTE line transfer problems andapplied it to protostar formation. For dynamical models, results from the radiation hydrodynamical calculations are used.
We conrm that the computational results show double-peaked proles in a stronger
blue peak for optically thick molecular lines. Optically thin lines show single-peakedproles with a slight blue asymmetry. These qualitative features are consistent with paststudies.
On the contrary to the remarks by a previous work we do not nd overestimation ofline widths. Furthermore, the infall motion produces wings extending to v = 2km/s inline spectra. Wings could not be produced in previous studies. These results imply that
simplied infall models such as the isothermal self-similar solutions adopted by previous
authors are not suitable to the detailed modeling of line spectra.
On the basis of the results and other theoretical and observational evidence, we illustrated an evolutionary picture of protostar formation. In terms of the evolutionary
time and the inclination to an observer, we nd that protostellar objects are successfully
categorized.

3

4

Contents
1 INTRODUCTION

9

2 NUMERICAL SCHEME

23

3 INPUT PHYSICS

39

1.1 Theoretical Studies of Protostellar Collapse . . . . . . . . . . . . . . . . . . 10
1.1.1 Historical Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.1.2 Technical Problems in Radiation Hydrodynamic Calculations . . . . 12
1.2 Observations of Young Protostars . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.1 General Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.2 Spectral Energy Distributions . . . . . . . . . . . . . . . . . . . . . 16
1.2.3 Molecular Line Spectra . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.3 The Aim of This Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1 The Basic Equations of Radiation Hydrodynamics . . . . . . . . . . . . . . 24
2.2 The Computational Method in Spherical Symmetry . . . . . . . . . . . . . 28
2.2.1 The Transfer Equation Solver . . . . . . . . . . . . . . . . . . . . . 28
2.2.2 The Moment Equations Solver . . . . . . . . . . . . . . . . . . . . . 30
2.3 The Energy Mean and Flux Mean of Opacities . . . . . . . . . . . . . . . . 34
2.4 Coordinates for Computation . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4.1 Lagrangian or Eulerian . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4.2 Grid Conguration . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1 Chemical Abundances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.1.1 General Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.1.2 Equations for Dissociation and Ionization Degrees . . . . . . . . . . 40
3.2 Specic Heat for Molecular Hydrogen . . . . . . . . . . . . . . . . . . . . . 45
3.3 Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.3.1 EOS for Ideal Gas . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.3.2 EOS for Non-ideal Gas . . . . . . . . . . . . . . . . . . . . . . . . . 50
5

CONTENTS
3.4 Opacities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4.1 Dust Opacities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4.2 Gas Opacities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4 THE FIRST COLLAPSE

57

5 THE WHOLE EVOLUTION OF PROTOSTAR FORMATION89

6 MOLECULAR LINE SPECTRA

113

4.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.1.1 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.1.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.1.3 Opacity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.4 Equation of State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2.1 Overview for a Typical Case . . . . . . . . . . . . . . . . . . . . . . 61
4.2.2 Characteristic Densities . . . . . . . . . . . . . . . . . . . . . . . . 68
4.3 Properties of the First Core . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.1 Denitions of the Core Radius . . . . . . . . . . . . . . . . . . . . . 71
4.3.2 Dependence on Cloud Mass . . . . . . . . . . . . . . . . . . . . . . 71
4.3.3 Dependence on Initial Temperature . . . . . . . . . . . . . . . . . . 73
4.3.4 Dependence on Opacity . . . . . . . . . . . . . . . . . . . . . . . . 76
4.4 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.4.1 Analytical Derivation of the Core Radius . . . . . . . . . . . . . . . 80
4.4.2 Estimate of the Accretion Luminosity . . . . . . . . . . . . . . . . . 82
4.4.3 Case for  = 7=5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.4.4 Comparison with Observations . . . . . . . . . . . . . . . . . . . . . 83
4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.1 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 90
5.1.1 Initial Density Structure . . . . . . . . . . . . . . . . . . . . . . . . 90
5.1.2 Initial Temperature Structure . . . . . . . . . . . . . . . . . . . . . 91
5.2 Overview of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3 Luminosity Curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.4 Spectral Energy Distributions . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.1 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6.1.1 Beyond LVG and Microturbulence Approximations . . . . . . . . . 114
6

CONTENTS
6.1.2 Basic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.1.3 Computational Procedure . . . . . . . . . . . . . . . . . . . . . . . 116
6.2 Calculations and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.2.1 Molecular Tracers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

7 AN EVOLUTIONARY PICTURE OF PROTOSTAR FORMATION 129
8 SUMMARY AND CONCLUSIONS

133

A Test Calculations for Radiation

137

B Thermal Relaxation Mode in Spherical Symmetry

141

C Conditions for Violation of Isothermality

145

C.1 Derivations for the Critical Densities . . . . . . . . . . . . . . . . . . . . . 147
C.2 The Minimum Jeans Mass for Fragmentation . . . . . . . . . . . . . . . . . 155
C.2.1 Comments on the Previous Studies . . . . . . . . . . . . . . . . . . 155
C.2.2 A Revised Denition of the Minimum Jeans Mass . . . . . . . . . . 156
C.3 SUMMARY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

7

CONTENTS

8

Chapter 1
INTRODUCTION
Investigation of star formation processes has been signicant subjects both for theoretical astrophysics and for observational astronomy. Nevertheless, many problems are left
unresolved still to date especially in an early phase in the evolution where a protostar isunder dynamical growth. The lack of our knowledge for protostar formation is a \missinglink" in the evolutionary picture of star formation.
Objects at a very early stage of the evolution are dicult to observe because of thefollowing two reasons: (1) a growing protostar is veiled by the infalling envelope, and (2)the timescale is relatively short for such young objects and then the number of the objects
expected to be observed is small. Recently, observations have made remarkable progressin the studies of protostar formation. Past theoretical studies on protostar formation,
however, were not designed to explain the recent observations by means of observable
quantities such as spectral energy distributions and molecular line spectra. Now a theoretical model for protostar formation is required which includes spectral synthesis toaccount for latest observations. This motivates me to do this work.
In this chapter, we present historical review on theoretical (x1.1) and observational(x1.2) backgrounds behind this work in order to clarify the problems to be resolved. In
x1.3 the aim of this work is summarized.

9CHAPTER 1. INTRODUCTION

1.1 Theoretical Studies of Protostellar Collapse
1.1.1 Historical Review

Protostellar evolution has been investigated actively during the past two decades. A
standard evolutionary scenario for the formation of isolated low-mass stars can be brieysummarized as follows
1. The scenario begins with the gravitational collapse of a molecular cloud core. Initially, the compressional heating rate by the collapse is much smaller than the coolingrate by the thermal radiation of dust grains. The cloud therefore collapses isothermally.
2. The isothermal condition is broken when the compressional heating rate overwhelms
the radiative cooling rate. The collapsing cloud thus forms an adiabatic core at its
center, which is the so-called \rst core". The evolution at this phase is called the\rst collapse".
3. When the central temperature reaches about 2000K, hydrogen molecules begin todissociate into atoms. Because released gravitational energy is consumed by the
dissociation, gas pressure cannot increase rapidly enough to support the rst coreagainst its self-gravity, and then the \second collapse" begins. After the dissociationis completed, the \second core" is born.
4. The surrounding envelope continues to infall through the \main accretion phase".
The gas cloud with signicant angular momentum forms a protoplanetary disk
around the central star. The infalling gas is thermalized by the shock at the surfaceof the star or the disk. The thermal radiation from the shocked region is absorbed
and reemitted by the envelope. Therefore, an observer can only see infrared emission
from the outer envelope until material in the envelope becomes suciently transparent to the radiation from the central objects. Bipolar outows are supposed to
be associated with this phase of evolution.
5. When the accreting envelope becomes transparent to optical light, the star andprotoplanetary disk appear as a T Tauri star. The dynamical evolution of a protostar
is followed by the quasi-static contraction of the central star and the hypotheticalviscous evolution of the protoplanetary disk.
Theoretical modeling for the above scenario are traced back to Larson(1969, hereafter
L69). He supposed, as a standard case, a slightly Jeans-unstable cloud of M = 1M
10

1.1. THEORETICAL STUDIES OF PROTOSTELLAR COLLAPSE
and T = 10K, and pursued gravitational collapse of the cloud by numerical calculation.
Dynamical evolution in the rst collapse phase was found to be well approximated by
the isothermal similarity solution, which was obtained independently by L69 and by Penston(1969) and is therefore called the \Larson-Penston (or LP) solution". In the mainaccretion phase, L69 solved the central protostar and the infalling envelope separately,
and connected them with each other by a discontinuous jump condition at the accretion
shock front. He showed that a protostar of 1M experiences the maximum values of
10R and 30L for the radius and luminosity, which nally settle at R 2:0R andL 1:3L. For an energy conservation relation at the shock front, he supposed that thethermal energy is radiated away with the eective temperature equal to the post-shock
temperature, which left a seed of controversy in the subsequent decade.
Hayashi(1970) and Narita, Nakano, & Hayashi(1970) pointed out that the isothermal
jump condition by L69 could overestimate the energy loss by radiation because it neglects
the self-absorption of radiation in the shocked layer. A protostar might be possibly larger
and more luminous with a more elaborate treatment for the shock front. Indeed, Narita,
Nakano, & Hayashi found that a protostar of 1M has a radius of 2  103R and a peakluminosity of 6  103L, which are greater by a few orders of magnitude than Larson's
results. This huge discrepancy, however, may be attributed to the dierence in the initialconditions, rather than to the dierence in treatments for the shock front. Narita, Nakano,& Hayashi initiated their calculation with a highly Jeans-unstable cloud and thus a large
amount of gravitational energy is available to form a very luminous protostar, while amarginally unstable cloud was supposed as the initial condition in Larson's case, where
the luminosity was found not to exceed a modest value.
In order to settle the controversy, a number of authors (Appenzeller & Tscharnuter
1975 Westbrook & Tarter 1975) attempted to resolve formidable diculties which arisein handling the accretion shock in numerical calculations. Tscharnuter & Winkler(1979,
hereafter TW) mentioned that the particular form of the von Neumann & Richtmyer articial viscosity was responsible for many of the diculties, and instead suggested a tensor
articial viscosity, which generalizes the von Neumann & Richtmyer viscosity to be applicable to hydrodynamical ows in general geometry. In use of the tensor viscosity, Winkler& Newman(1980a,b) (collectively WN hereafter) succeeded in overcoming the numericaldiculties and they obtained the results without any simplications on the treatment forthe shock front such as the isothermal jump condition. WN found their results in the
protostellar properties in good agreement with L69. Stahler, Shu, & Taam(1980, hereafter SST) also reached a consistent conclusion, although the policies in computational
procedure were totally dierent among Larson, WN, and SST (see x1.1.2). The dispute
11

CHAPTER 1. INTRODUCTION
seemed to cease eventually.
Tscharnuter(1987), however, objected to the authentic scenario for protostar formation. He found that a formed protostar is subject to a non-linear oscillation over twenty
orders of magnitudes(!) in central density during the whole evolutionary timescale, andthe steady accretion phase never appears. The only case that a protostar can remain
dynamically stable is when the initial cloud core is more massive by several orders ofmagnitude than the Jeans mass. Balluch(1988) conrmed the result by Tscharnuterand analyzed in detail the physical mechanism which destabilizes protostars (Balluch
1991a,b,c), but the problem remains unresolved to date. There are two possible sourcesof instabilities: (1) updated tables of EOS(equation of state) and opacities and/or improved numerical schemes enabled to discover physical instability of protostars, which had
not been found in the former calculations, and (2) calculations performed by Tscharnuterand Balluch suer from a kind of numerical instabilities which were possibly due to acomputational technique newly introduced in their numerical codes. We will try to settle
this controversial problem in this work.
In the next subsection we review technical problems that lie in numerical calculations
for protostellar collapse, especially from the viewpoint of radiation hydrodynamics.

1.1.2 Technical Problems in Radiation Hydrodynamic Calculations
Theoretical investigations of the above scenario involve accurate solutions of twodimensional or three-dimensional radiation magneto-hydrodynamical equations, which
inevitably require numerical calculations. At present, however, the computational facilities do not have sucient power for this task, and therefore a number of simplications are supposed in such investigations. A possible way is multi-dimensional (magneto)hydrodynamic approach, supposing the isothermal or polytropic equation of state insteadof taking radiation transport into account. This approach is suited to investigate angularmomentum transport through gravitational or magnetic torque in collapsing clouds, but
is inapplicable when main interests of the author lie in the thermal evolution such asformation of protostars and temperature structures in the protostellar envelopes of gas
and dust. In the latter case, radiation hydrodynamic calculations are inevitable. In thiswork, we concentrate our interests on the radiation hydrodynamic approach.
In the past studies which dealt with the radiative energy transport, the diusion (orEddington) approximations and/or the grey approximation are often introduced for simplicity. For instance, both were assumed by L69, Appenzeller & Tscharnuter(1974). In
recent works, Boss & Myhill(1992) introduced the grey 3D Eddington approximation(Boss121.1. THEORETICAL STUDIES OF PROTOSTELLAR COLLAPSE
1984) for incorporating the radiative energy transport into their radiation hydrodynamical code. The diusion or Eddington approximation, however, is justied only in optically thick regions such as stellar interiors. Considering this point, Narita, Nakano, &Hayashi(1970) adopted the diusion scheme only for the inner opaque region and solved
the angle-dependent radiative transfer in the tenuous and extended envelope. Besides theclassical diusion approximation, the ux-limited diusion (FLD) approximation is introduced by some recent works (Yorke, Bodenheimer, & Laughlin 1993,1995). FLD schemesaves the radiative diusion equation from yielding unphysical solutions in optically thinregions, but still it is not ensured to give the correct solutions.
The exact treatment of the radiation elds in radiation hydrodynamics (hereafterRHD) requires some drastic simplications in order to reduce the computational eort.
Indeed, the past works which dealt with the radiation elds can be classied into two kindsof categories, i.e., (i) two-dimensional radiative transfer without hydrodynamics or (ii)
one-dimensional RHD. The former (i) includes Efstathiou & Rowan-Robinson(1990,1991)
and Chick, Pollack, & Cassen (1996). They calculated the temperature distributions
which are consistent with the radiation elds in axisymmetry. In their calculations thedensity distributions were xed as snapshots of the dynamical evolution, instead of solving
the hydrodynamical equations. Moreover, the density distributions assumed by them were
solutions of isothermal collapse and therefore the temperature and density distributions
are not self-consistent in their results.
The latter (ii) gives the self-consistent solutions with respect to all the radiation hydrodynamical quantities. A couple of comprehensive works which belong to this category
were published in the same year, SST and WN. SST assumed steady-state inow in theenvelope and divided the problem into several parts, which were solved individually andthe solutions in each part were connected at the interfaces. Although SST successfully
overcame the diculties that appear in numerical calculations for protostar formation(x1.1.1), their treatment for the radiation eld, especially in the radiative relaxation layer
ahead the accretion shock front, still included a number of simplications. In contrast,
WN pursued the dynamical evolution by solving the radiation hydrodynamical equations
straightforwardly using the variable Eddington factor method, which is employed also inthis work. Both of them yielded consistent results with each other, but cases other than a
1M cloud with a single initial condition were left to be examined for future investigators.
Furthermore, both of them introduced the grey approximation, in which the energy-meanand ux-mean opacities are replaced by the Rosseland mean opacity instead of solving thefrequency-dependent radiative transfer equation. The grey approximation is valid where
the temperature varies spatially in a scale much larger than a mean free path of thermal
13CHAPTER 1. INTRODUCTION
photons. This is, however, not the case at the latter stage of the protostellar evolution(i.e., in the main accretion phase). Indeed, Preibisch, Sonnhalter, & Yorke(1995) showedthat the grey approximation can cause considerable errors in estimating the temperature
distribution and the radiative acceleration where a optically thin envelope is heated upby the central hot star.
Frequency-dependent calculations for radiation, without the grey approximation, were
performed by Yorke(1979,1980) for massive star formation. He solved the hydrodynamical
equations simultaneously to see the structural evolution of a circumstellar envelope, but
did not resolve the structure of the central star, putting a central radiative source there
instead.
In this work we adopt spherical symmetry, which inevitably requires that the eectsof rotation and magnetic elds should be neglected. Both the eects, of course, actually
have important roles in protostellar evolution. Nevertheless, the spherically symmetriccalculation still has signicance because of the following two reasons. One is the dicultyin the exact treatment of radiation eld in general geometry. Another reason is thatmoving coordinates such as the Lagrangian mesh are available only for one-dimensional
calculations. A molecular cloud core of  10;19 g cm;3 and R 0:1pc contracts to a starof  1g  cm;3 and R 1R, where the density increases by 20 orders of magnitude andthe radius decreases by 6 orders. Any multi-dimensional numerical calculations to date,that include adapted mesh renement technique and smoothed particle hydrodynamics,
have not succeeded to cover such a wide dynamic range in a self-consistent way.
As a kind of moving coordinates, one may prefer an adaptive grid scheme, where
grid points are moved freely to achieve a most desirable grid conguration to resolve thephysical structure. WN carried out a calculation with an adaptive grid based on thescheme developed by TW. While their work successfully pursued the whole evolution of
protostellar collapse, TW reported that a number of empirical techniques for stabilization
were essentially important for avoiding numerical instabilities. More recently, Tscharnuter(1987) and Balluch(1988,1991a,b,c) adopted an adaptive grid scheme developed byDor & Drury(1987). This sophisticated scheme is quite helpful to resolve detailed structures in shock fronts, but its stiness against numerical instabilities is still uncertain.
As Balluch(1991a) pointed out, the possibility cannot be denied that the freely moving
coordinate may be responsible for the dynamical instabilities found by Tscharnuter and
Balluch (x1.1.1). In order to avoid unknown numerical instabilities, we use no adaptive
mesh schemes and instead employ either Eulerian or Lagrangian coordinates (x2.4).

14

1.2. OBSERVATIONS OF YOUNG PROTOSTARS

1.2 Observations of Young Protostars
1.2.1 General Remarks

In early 1980's Infrared Astronomical Satellite (IRAS) discovered a number of sources
which radiate mostly in the far- to near-infrared and are invisible in optical light. Heavy
reddening and large infrared excess implies that the central star and circumstellar disk
are highly obscured by the opaque envelope of gas and dust, where visible photons fromthe central objects are reprocessed into colder thermal emissions in the infrared. Such ayoung star under growth through the mass accretion from the infalling envelope is called a\protostar". This picture is supported by the idea that the objects are so young that themass infall onto the central star has not ceased yet (Adams, Lada, & Shu 1987 Kenyon,
Calvet, & Hartmann 1993). The evolutionary sequence from embedded protostars to
evolved T Tauri stars is well categorized in terms of spectral energy distribution (Lada
1987):
classication spectral index
class I0<a<3class II
;2 < a  0class III
;3 < a < ;2typical correspondents
embedded protostars
classical T Tauri stars
weak-line T Tauri stars

where the spectral index is dened by a  d log(F )=d log() based on the near- andmid-infrared ux F .
Class I sources, however, were found to already have an evolved protostar in the center
while only a minor fraction of the total mass remains in the infalling envelope. Younger
protostars which are actively accreting the mass from the infalling gas were not known
until the discovery of VLA 1623 (Andre et al. 1993). VLA 1623 is the driving source ofa remarkable CO bipolar outow, which strongly indicates the existence of a protostar
in the center, but a protostar itself is undetected in the near- or far-infrared and thenVLA 1623 shows only cold components of T =15-20K in the spectral energy distribution.
Andre et al. proposed a new category for protostellar objects younger than class I as\class 0", which is dened as submillimeter sources with Lbol=L1:3 < 2  104, where L1:3
is the 1.3mm luminosity. Since the discovery of VLA 1623, more than 20 objects havebeen identied as \class 0" (e.g., Gregersen et al. 1997).
In addition to class 0 objects, observations of \star-less" cloud cores are also importantto reveal early stages in protostellar evolution. Mizuno et al.(1994) showed that the cloud
cores without young stars are less dense and more extended than those with stars, andderived the timescale of core collapse to be a few 105 yr. Onishi et al.(1998) presented
results of a complete survey of C18 O cores in Taurus molecular cloud and discussed on
15

CHAPTER 1. INTRODUCTION
general properties of star-forming cloud cores. Observational studies are now ready forfurther understanding of early protostellar evolution from star-less cores to class I objects
through class 0 stage. In contrast, theoretical support is far less sucient due to lack ofreliable models which can be compared directly with observations.
Recent studies using radio interferometers have succeeded in direct imaging of thephysical structures in the infalling envelopes around protostars (e.g., Saito et al. 1996
Ohashi et al. 1997 Momose et al. 1998). In addition, spectral energy distributions(SEDs)
and molecular line spectra are ecient tools for detailed analysis of observational data on
the basis of theoretical models. We summarize characteristic features in SEDs and line
spectra observed from protostellar objects in the rest of this section.

1.2.2 Spectral Energy Distributions
For YSOs in early phases of evolution, direct observation of the stellar surface is impossible
due to large extinction of the circumstellar envelope and the evolutionary tracks on HR
diagram loose its signicance in contrast to T Tauri stars and more evolved stars. WN
and SST drew the evolutionary tracks for embedded protostars on HR diagram using
the eective temperature of the dust photosphere (i.e., the temperature where the opticaldepth measured from an observer by the dust opacity equals unity or 2/3) for the abscissa.
However, this denition of Te cannot be interpreted directly to observable quantities, as
pointed out by WN themselves, and then another well-dened measure is required.
In order to examine the physical state of the whole system of the central objects andsurrounding envelope, spectral energy distributions (SEDs) are helpful because SEDs,
in principle, sensitively reect physical structures of the objects such as temperatureand density distributions. Indeed, the well-known classication of YSOs by Lada(1987)
mentioned above | i.e., class I,II, and III| is based on the shapes of SEDs in the nearand mid-infrared.
class 0 sources show cold SEDs which can be tted with greybody spectra of
T a few  10K (Andre et al. 1993 Ward-Thompson et al. 1995) although they
have evidence of protostellar activity such as molecular outows (Barsony 1994). These
cold SEDs seem to indicate that class 0 objects are substantially younger than class I
protostars. However, the absence of infrared counterparts disables direct extension of
Lada's classication to class 0. As a new quantitative indicator of evolutionary status of
YSOs, Myers & Ladd(1993) proposed the bolometric temperature, Tbol , dened byR1

0 B (Tbol )d = 
0 B (Tbol )d

R1

16

(1.1)

1.2. OBSERVATIONS OF YOUNG PROTOSTARS
where

R1

 d
  R0 1F
F d

(1.2)

h = 1:25  10;11 Hz;1K
Tbol  4 (4)
(5) k

(1.3)

0



on the basis of an observed SED, F . Combining equations (1.1) and (1.2), one obtains

where (n) is the Riemann zeta function of argument n:
1 Z 1 tn;1
(n) 
dt:
;(n) 0 et ; 1(1.4)

Myers & Ladd(1993) showed that Tbol successfully distinguishes evolutionary positions of
YSOs from protostars to weak-line T Tauri stars: Tbol = 285  54K for 23 protostars,
721  68K for 12 protostars with faint optical counterparts, 2128  105K for 57 classical T
Tauri stars, and 3552  139K for 37 weak-line T Tauri stars. Gregersen et al.(1997) foundthat the mean value of Tbol is 44K for 16 class 0 sources. The fact that class 0 sources
show the coldest bolometric temperature supports that class 0 is in an earlier phase thanclass I in protostellar evolution.
Nevertheless, the evolutionary position of class 0 is still ambiguous. An edge-on viewof a class I object would be observed as if it were class 0 because of the larger extinction
than a pole-on view. Indeed, almost all class 0 sources are associated with outows whoseaxes are nearly perpendicular to the line of sight (Hirano et al. 1997), which implies thatclass 0 sources are observed typically on edge-on views. Andre et al.(1993) evaluated104yr for the age of VLA 1623 but the estimate includes large ambiguity. The actual ageof class 0 sources should be determined by more detailed theoretical analysis.
While the activities occurring deeply in the opaque envelopes do not appear in SEDs,
molecular line spectra reect them more sensitively. The next subsection is devoted to
the review of observations of line spectra in YSOs.

1.2.3 Molecular Line Spectra
Observations have revealed that protostellar cloud cores often show characteristic prolesin molecular line spectra. In particular, a double-peaked prole with a stronger blue peakin an optically thick line is often observed, which is considered as strong evidence for theinfall motion. Figure 1.1 illustrates a schematic picture to explain how an asymmetricdouble-peaked spectrum is formed. Iso-velocity contour for the velocity projected to theline of sight is shown, where solid curves represent the red-shifted component and dashedcurves are the blue-shifted when the cloud is collapsing (or infalling). For brevity we
17

CHAPTER 1. INTRODUCTION

Blue

Red

Fig. 1.1.| A schematic picture for describing the reason why a blue asymmetry appears
in line spectra of a collapsing (or infalling) cloud. Contour curves represent the velocity
eld projected to the line of sight, where both the red-shifted (solid) and blue-shifted(dashed) components are shown. See text for more detail.

18

1.2. OBSERVATIONS OF YOUNG PROTOSTARS
assume spherically symmetric collapse and neglect the rotation of the cloud and bipolaroutows. Paying attention to a certain pair of contour curves of positive and negative
velocities with the same absolute value, one nds that there are a couple of intersections
by the line of sight for each contour curve of the red-shifted and blue-shifted velocities.
The nearer intersection to an observer locates a further side from the cloud center for the
red-shifted component, while it is closer to the center for the blue-shifted. The excitation
temperature decreases outward in general, and hence the blue component in a spectrum
should be more intensive than the red one if the observed molecular line is optically thickenough. Spectral proles of an opaque line should have a deep absorption feature at theline center because the bulk of gas that remains cold and almost at rest near the edge ofthe cloud absorbs the radiation from the active region in the cloud center.
Although the original idea is traced back to Snell & Loren(1977) and Leung &Brown(1977), it was not until Zhou(1992) that line features as evidence of infall motions began to be investigated systematically. Detailed modeling of line proles for B335
was done and a double-peaked spectrum with a stronger blue peak was conrmed to be agood qualitative signature of infall (Choi et al. 1995). Recently, some statistical investigations were done for protostellar candidates including class 0 and class I sources (Gregersen
et al. 1997 Mardones et al. 1997). Mardones et al. (1997) revealed that line proles inclass 0 sources are statistically more skewed toward the blue-shifted velocity than in class
I, which implies that the line asymmetry is a good indicator of the evolutionary status atearly stages in star formation. On the other hand, current models for infall motions do
not account for a red asymmetry and a wing component, and hence spectral proles ofindividual sources are not always modeled successfully. Those features can be attributed
to outows, as is usually done, but no theoretical models are available to t observations
including contamination of outows.
As infall models to reproduce the observed line proles, current works adopt theisothermal similarity solutions such as the LP solution (see x1.1.1, denoted by LP model
hereafter) and/or the Expansion-Wave collapse solution derived by Shu (1977, Shu model
hereafter). Zhou (1992) calculated emergent line spectra for these two models, and hereached the conclusion that LP model overestimates the velocity width of CS lines incomparison with observations 1 . This seems, however, to arise from an unrealistic condition in \naive" LP model, where the infall velocity distribution extends to a nite value
1 The denotations of \L-P ow" and \Shu model" by Zhou might cause a little confusion. \L-P ow"
by Zhou includes, in addition to the original LP solution derived by L69 and Penston(1969) for t < 0,
the continuous counterpart of the LP solution after t exceeds zero. This latter solution was obtainedby Hunter(1977). On the other hand, \Shu model" contains a sequence of the hydrostatic equilibrium
solutions for t < 0 and the Expansion-Wave collapse solution for t > 0. Here t = 0 denotes the instancewhen the central density diverges to innity, or when the central star is supposed to be born.

19

CHAPTER 1. INTRODUCTION(= 3:26cs) in the limit of innity. On the other hand, numerical experiments by radiationhydrodynamic (RHD) calculations (L69 Masunaga, Miyama, & Inutsuka 1998, hereafter
MMI) have found that the initial phase of protostellar evolution is excellently approximated by LP model under reasonable boundary conditions supposing a xed boundary ora constant external pressure at a nite radius. This \realistic" LP model may provide adierent spectral prole from that obtained by Zhou(1992). The Expansion-Wave collapse
solution is believed to provide a more reasonable model, but it might be still too idealized
for detailed spectral synthesis. The applicabilities of the isothermal similarity solutions
remain unexamined, and thus another study on spectral line proles is necessary using
an elaborate dynamical model. Detailed modeling of line spectra is one of the purposesof this work, as is mentioned in the succeeding section.

20

1.3. THE AIM OF THIS WORK

1.3 The Aim of This Work
Main goals of the thesis are summarized as follows.
1. The nal goal of this work is to give an evolutionary picture of protostar formation
which have both theoretical reliability and accountability for latest observations.
For this purpose, radiation hydrodynamic numerical calculations are performed.
We have developed the numerical scheme so that none of diusion, Eddington, andgray approximations are required, by solving the frequency- and angle-dependent
radiative transfer equation.
2. Input physical processes such as the non-ideal equation of state, the opacity tables
for atomic, molecular, and dust opacities have been updated since the previous
studies done before 1980. Incorporating the newest data, we carry out numericalcalculations and discuss on the consistency with past studies.
3. Computational results of RHD numerical calculations provide constraints on theoretical models for protostar formation, through comparing the results with observational data. The spectral energy distributions and the luminosity curve are helpful
tools to the comparison of models with observations, which gives, for instance, aconstraint on the initial conditions of protostellar collapse.
4. On the basis of the results of RHD calculations, we perform line spectral synthesis
by non-LTE line transfer calculations. While all previous studies adopted simplied dynamical models such as the isothermal self-similar solutions, more elaboratemodeling is enabled in our calculations. The computational results provide realistic
constraints on the possibility that observed line spectra are explained by the infall
motions.
In order to discuss on the above subjects, we rst introduce our numerical scheme
for RHD calculations in Chapter 2. The physical processes that are incorporated into
the numerical code are described in Chapter 3. The subsequent chapter (Chap. 4) isdevoted to discussions on the properties of an early stage of protostellar collapse, i.e.,the rst collapse, under some simplications on the input physics, and we will show thatthe numerical results are supported by analytic estimations. Chapter 5 is a main partin the thesis, where the computational results are presented for the whole evolution of
protostar formation and discussions are made about possible constraints on theoreticalmodels that reasonably account for observations. In Chapter 6, numerical calculations ofline spectral synthesis are carried out using the dynamical model presented in Chapter 5.
21

CHAPTER 1. INTRODUCTION
We investigate there typical features in spectral proles which are produced by the infall
model, and applicability of the isothermal self-similar solutions to line spectral synthesis
is critically examined.
Considering also possible non-spherically-symmetric eects, we illustrate an evolutionary picture of protostar formation based on our results in Chapter 7. The ndings in thewhole paper are summarized in Chapter 8.
In appendix A test calculations for the numerical codes are given, where the thermal relaxation mode for radiation is employed. The analytical solution of the thermal
relaxation mode for the spherically symmetric eigen function is derived in Appendix B.
Appendix C is devoted to the derivation of the critical density at which the isothermal
evolution is terminated. Some important consequences that arise from the considerationon the critical density are discussed as well.

22

Chapter 2
NUMERICAL SCHEME
In this chapter, our computational scheme for radiation hydrodynamic (RHD) calculations
is described in detail. In ordinary hydrodynamics the mean free path of particles issuciently small compared to the typical length scales with which physical quantities
vary. The Boltzmann equation therefore falls onto more feasible equations such as theEuler or Navier-Stokes equations. Analogously, the radiative transfer equation falls onto
the radiative diusion equation as long as the mean free path of photons is much smaller
than the typical scale of the system, in other words, if the medium is optically thick. This,
however, condition is not satised in general, which makes RHD hard to deal with.
In this chapter we introduce the formulation of the RHD equations and the computational scheme to solve them with best accuracy and least computational cost.

23

CHAPTER 2. NUMERICAL SCHEME

2.1 The Basic Equations of Radiation Hydrodynamics
The basic equations of radiation hydrodynamics (RHD) are introduced in this section.
Hydrodynamical equations which account for the interactions with radiation areD + r v = 0
(2.1)
Dt
 DDtv = ;r! ; rp + cF 0 F 0
(2.2)
 De
(2.3)
Dt + pr v =  CR ; "gd + (c E0E0 ; 4 P 0B )
r2! = 4 G
(2.4)
p  p( e)
(2.5)
where ,v, e, and T are the uid mass density, the velocity, the internal energy density
per unit mass, and the temperature, respectively. , CR , "gd , and  are the absorptioncoecient, the cosmic ray heating rate per unit mass, the rate of energy exchange perunit mass between gas and dust, and the ratio of specic heats, respectively. B  T 4=
denotes the Planck function. E0 and F 0 , where subscript 0 denotes comoving variables,
are the frequency-integrated radiation energy density and radiation ux, which are dened
as
Z 1I
E0 = d0 d 0 I0(0  n0 )
0

Z 1

I

d0 d 0 I0(0  n0 )n0
where I0(0  n0) denotes the specic intensity of radiation at frequency 0 along the direction vector n0 in the comoving frame. The frequency-averaged absorption coecientsare dened by
R1= 0 0(0 )F 0 (0)d0 
F0

=

0

F0

E0 =

F0

R1

0 (0 )E0 (0 )d0E0R10 0 (0 )B (0 )d00:BEquations (2.1),(2.2),(2.3),(2.4), and (2.5) correspond to the equation of continuity,
the equation of motion, the rst law of thermodynamics for material, Poisson's equation,and the equation of state, respectively. They are the ordinary hydrodynamical equations
P0 =

24

2.1. THE BASIC EQUATIONS OF RADIATION HYDRODYNAMICS
except for the last terms on the RHS of equations (2.2) and (2.3) which describe interactions between material and radiation. To account for the radiation eld exactly, thetransfer equation of radiation (2.6) must be solved.
1 @I ( r n) + n rI ( r n) = ; ( r)(I ( r n) ; S ( r))
(2.6)
c @t
where S ( r) is the source function. Equation (2.6) is written in the lab frame.
Because the specic intensity of radiation depends on frequency and on direction angle
as well as on location, it requires considerable computational eort to solve equation(2.6) directly. Moreover, equation (2.6) must be transformed into the comoving frame tocouple with the hydrodynamical equations, which would modify the equation into a quite
complicated form. In contrast, it reduces the computational eort to an acceptable level tosolve the moment equations of radiation in the comoving frame (eqs. #2.7] and #2.8] below)
instead of the transfer equation itself. The moment equations of radiation in the comoving
frame were originally derived by Castor(1972) and Buchler(1979) and their studies weresucceeded by Mihalas & Mihalas in a comprehensive work(1984, hereafter MM), and by
Stone, Mihalas, & Norman(1992, hereafter SMN) as a two dimensional numerical RHD
code (the nal work of the \ZEUS-2D" trilogy). The formulation we adopt is taken
after MM and SMN, who dropped the terms which are O(v=c) in all physical regimes.
Therefore the adopted equations can account for all the dynamical eects of radiation ina non-relativistic uid (see SMN).

 D ( F 0 ) + r P = ; F 0F 0 0
c2 Dt c
D ( E0 ) + rv : P + r F = ;(c E ; 4
 Dt
0
0
E0 0


(2.7)
P 0 B )(2.8)

where P 0 is the radiation stress tensor:
P0

=

Z 1

0

I

d0 d 0I0(0  n0)n0 n0:Equations (2.7) and (2.8) correspond to the momentum equation and the energy equation of radiation, respectively.
In our scheme, equation (2.8) has been replaced with the following one.

D (e + E0 ) + rv : P + r F + pr v =  
 Dt
0
0
CR


(2.9)

which is given by taking the sum of equations (2.3) and (2.8). Equation (2.9) ensures thatthe solutions will be more accurate than equation (2.8) when the uid is near radiative25

CHAPTER 2. NUMERICAL SCHEME
equilibrium because the interaction terms between material and radiation have cancelled(SMN).
Equations (2.7) and (2.8) are not closed in themselves since they have been integrated
over direction angle. Therefore, a closure relation must be given which compensates thelack of knowledge for the angular distributions of the specic intensity. As a closure
relation, we introduce the tensor variable Eddington factor f , which is dened byP= f E:(2.10)The radiation stress tensor P is eliminated using the tensor variable Eddington factor,
which is a measure of the degree of anisotropy of the radiation eld, to close the systemof the RHD equations.
When we calculate f , the radiation eld at each time step is assumed as a static snapshot under the instantaneous density and temperature distributions. In other words, wesolve equation (2.6), neglecting the term including the time derivative. This assumption
is justied only when the time scale of radiative diusion or propagation is suciently
short compared to the evolutionary time scale of interest. Fortunately, such a conditionis satised in many astrophysical problems including protostellar collapse. Even if thecondition is not satised, the static transfer equation is still applicable for most cases
in astrophysics because the solutions from the transfer equation contribute only to the
Eddington factors, which should vary less than the individual values of E and F (MM).
We also note that the terms of O(v=c), which would appear when transforming equation(2.6) into the comoving frame, have been dropped because the relativistic eects such as
aberration or Doppler shift hardly aect the Eddington factors for the problems of ourpresent interest. Thus, the transfer equation is reduced to a fairly simpler form and thenno longer bothers us for the unacceptable computational eort.
Once I ( r n) is obtained, the monochromatic tensor variable Eddington factor is
calculated using the relation:
P ( r )f ( r ) =E ( r) :The frequency-integrated Eddington factor is dened by
R1( )E ( )d
:E
The computational method, which is called the variable Eddington factor method,
was established as a technique for astrophysical RHD rst by TW. There is no essential
dierence in the computational scheme between TW and us, except that they supposed thegrey approximation while we solve the frequency-dependent radiative transfer equation.
f

= 0

f

26

2.1. THE BASIC EQUATIONS OF RADIATION HYDRODYNAMICS
This method can be divided into two parts. In the rst part, the transfer equation (2.6)
with @=@t  0 is solved under the density and temperature distributions assumed tobe known. Thus this part is called the transfer equation solver. The tensor variable
Eddington factor at each spatial point is obtained from the solutions I ( r n). In thesecond part, the full set of the RHD equations (2.1)-(2.4), (2.7), and (2.9) is solved using
the given Eddington factors. This part is called the moment equations solver. The newly
obtained density and temperature distributions provide improved Eddington factors forthe rst part, which is employed again. The two parts are iterated in this manner, forconsistency. In x2.2, we will show the details of each part for the spherical symmetricalcase. Test calculations are presented in Appendix A.

27

CHAPTER 2. NUMERICAL SCHEME

2.2 The Computational Method in Spherical Symmetry
2.2.1 The Transfer Equation Solver

We solve equation (2.6) numerically, dropping the time derivative term, in the transfer
equation solver to determine the variable Eddington factors. The transfer equation to besolved is reduced to@j = ;h 
(2.11)@@h = ;j + S (2.12)


@
where
(2.13)j  12 (I (+) + I (;))h  21 (I (+) ; I (;))
(2.14)
Z  ext(2.15)
 ds
X
X(i)sca J ):(2.16)S  1ext ( (i)abs B (T (i) ) +




(i)

(i)

(i)abs
The extinction opacity, ext
 , is the sum of the absorption and scattering opacities, and (i)sca respectively:
ext = X( (i)abs + (i)sca )



(i)

where the superscript (i) represents each of gas and three components of dust grains (seex2.2.2). An independent variable  is noted as a subscript hereafter.
Figure 2.1 illustrates the paths along which equations (2.11) and (2.12) should be
integrated. The specic intensity on each point at which an integral path intersects aradial grid is integrated over  on each ring to calculate the rst three moments forradiation,
Z 1
J = j d
(2.17)
0

H =

Z 1

h d

0
Z 1
K = 2j d
0(2.18)(2.19)
where  is the direction cosine, i.e.,  = cos , where  is the angle with respect to theoutward normal(see Figure 1). The rst two moments for radiation, J and H , correspond28

2.2. THE COMPUTATIONAL METHOD IN SPHERICAL SYMMETRY

p

I+
r
I-

ο

θ

x

Fig. 2.1.| The grid employed for the discretized transfer equation. The transfer equation
is integrated along the horizontal line, according to the orthodox manner in solving radiative transfer numerically in spherical symmetry. The cell-centered points (lled circles)
and cell-faced points (open squares) are employed for the discretized transfer equations
for j and h , respectively. Direction cosine, , is dened as  = cos , where  is theangle between the lines o-x and o-r.

29

CHAPTER 2. NUMERICAL SCHEME
to the monochromatic radiation energy density and radiation ux, respectively. The third,P or (4 =c)K , is the radiation pressure, which is the rr-component of the radiation stresstensor in spherical coordinates. Both the - and -components equal P ; (3P ; E )=2(see MM).
The scalar variable Eddington factor isf (r) = KJ ((rr)) :When the radiation eld is nearly isotropic, i.e., I () = I0 + I1, K is found to equal(1=3)J using equations (2.17) and (2.19). On the other hand, the streaming limit condition
I () = I( ; 1) yields K = J . The Eddington factor thus varies between 1/3 and 1,
except for some special cases where f goes down below 1/3 (cf., x5.2), according to thedegree of anisotropy of the radiation eld.
We adopt the monochromatic transfer equation solver in our scheme. It has two
advantages over the grey transfer equation solver. First, the energy mean and ux mean
opacities E and F for the moment equations are obtained (see x2.3), whereas the Planckmean or Rosseland mean opacities are substituted in the grey approximation. Anotheradvantage is that one can obtain SEDs directly from the monochromatic radiation uxes.
The obtained SEDs will help us to compare theoretical models with observations.

2.2.2 The Moment Equations Solver
In spherical symmetrical Lagrangian (comoving) coordinates, equations (2.1)-(2.4), (2.7),
and (2.9) lead to(MM)dr = v
(2.20)dtr = 43 @@M(2.21)(r3 )Dv = ; GMr ; 4 r2 @p + gas
F F ; vQ 
(2.22)
Dt
r2
@Mr c

D (e + E ) + (fE + p) D ( 1 ) ; (3f ; 1) v ]E = + ; " ; @ (4 r2F )  (2.23)
CR
Q
gd
Dt

Dt 
r
@Mr
1 DF + 4 r2 @ (fqE ) = ; F F (2.24)
c2 Dt
q @Mr
c
De + p D ( 1 ) = (c gas E ; 4 gasB ) + + ; " (2.25)
CRQ
gd
EP
Dt Dt 30

2.2. THE COMPUTATIONAL METHOD IN SPHERICAL SYMMETRY
Here vQ and Q are taken from the tensor articial viscosity in the formulation by TW,2@ (r2 v) )]
vQ = 9r @r@ 3 #r3l2  @ (@rr 3v) ( @v;@r @r3Q = ;4:5l

22
2 @ (r v ) ( @v ; @ (r v ) )2 
@r3 @r @r3

where l, a coecient with dimension of length, is of the order of the local width of gridto smear out the shock front articially so that the discretized hydrodynamic equations
can handle shock waves. Both vQ and Q are switched into zero where r  v > 0, that is,
the articial viscosity works only when the gas is compressed.
The sphericality factor q is dened by(Auer,1971)
3f ; 1 dr0:
Rin fr0
These equations are closed when the scalar Eddington factor and the equation of state(2.5) are given. The internal energy density e and the gas pressure p are replaced by T
and  using the equation of state. The gas opacity, gas , are treated as absorption opacity.
Heating or cooling by collisions between gas and dust is taken into account throughthe energy exchange rate, "gd:X
"gd = "(gdn) 
ln q =

Z r(n)where the superscript (n) denotes each of the dust components. We consider three components of dust grains, i.e., amorphous carbon, silicate, and dirty ice-coated silicate, on thebasis of the dust model by Preibisch et al.(1993) (see also x3.4.1). The energy exchange
rate for each dust component is expressed as (Hayashi & Nakano 1965 Leung 1975)"(gdn) = 3 ( 8k )1=2 nH2 n(grn) ( a(n)2 )T 1=2 k(T ; Tgr(n))
2 mH2


(2.26)

where n(grn) and Tgr(n) are the number density and temperature of dust grains of component(n). The dust temperature is determined by (Hayashi & Nakano 1965),"(gdn) + (En)abs cE ; 4 (Pn)abs Tgr(n)4 = 0:(2.27)The nuclear energy productions and the convective energy transport are not includedin the present version of the numerical code.
The nine variables r  v T E F and Tgr(n) (n = 1 2 3) are solved as functions ofmass. For spatial dierencing of the equations,  T Tgr(n) and E are dened at the centerof each grid cell, and r v and F are on the interfaces which bound the cell. We thus31CHAPTER 2. NUMERICAL SCHEME)
obtain vectors of the variables #ri ,i+ 12 ,vi ,Ti+ 21 ,Tgr(ni+
1 , Ei+ 1 ,Fi ], which are to be solved
22
in a fully implicit scheme. The system is linearized around the trial estimates #ri,i+ 12 ,vi),Ti+ 21 ,Tgr(ni+
1 , Ei+ 1 ,Fi ] to obtain a (I  I ) block tridiagonal system of (9  9) matrices,22)where I is the number of grid cells. The solution vectors #ri ,i+ 21 ,vi ,Ti+ 12 ,Tgr(ni+1,2Ei+ 12 ,Fi] are solved by Gaussian elimination.
In contrast to implicit method, explicit methods are constrained by the CFL condition,
which assures numerical stability(Courant et al.,1928),(jvj + c)%t  %r

(2.28)

where v c %t and %r are the uid velocity, the speed of transport of information, thetimestep, and the width of the grid, respectively. Solving the moment equations of radiation on an explicit scheme would require an unacceptably small value of %t comparedto the evolutionary time scale because c in (2.28) is the speed of light in this case. Then,the moment equations of radiation (2.24) and (2.25) with the uid energy equation (2.23)
are often solved on an implicit scheme while the uid motion can be solved explicitly forsome cases. However, the CFL condition which rules the hydrodynamic equations |i.e.,equation (2.28) in which c is now replaced by the speed of sound| still imposes sucha severe restriction that explicit methods are ruled out in numerical calculations for thewhole evolution of protostellar collapse (TW). It is therefore imperative that the full setof equations (2.20)-(2.25) is solved with an implicit scheme.
The frequency-averaged absorption coecients E  F  and P and Eddington factor
fE are needed in order to solve equations (2.20)-(2.25). The weighting values E and F ,or J and H , are given by direct integration of I over . We explain the method for
calculating J and H afterward. The frequency-averaged values are obtained using thesolved spectral proles J and H as
R1
0  H d 
(2.29)F=
H
R1
0  J d 
(2.30)E=J
R1
0  B d 
(2.31)
P=
T 4=
and
R1
(2.32)
fE = 0 fJ J d 
where
Z 1
J = J d
032

2.2. THE COMPUTATIONAL METHOD IN SPHERICAL SYMMETRY
Z 1

H = H d:0
The entire procedure at each time step is carried out as follows.
1. Trial estimates of density and temperature structures ((r) and T (r)) are assumed.
2. Equations (2.11) and (2.12) are solved under the given density and temperature
distributions to calculate the rst three moments for radiation (2.17),(2.18), and(2.19).
3. The integrations over frequency (2.29)-(2.32) are carried out to obtainfE (r), E (r), F (r), and P (r).
4. Equations (2.20)-(2.25) are solved using the given Eddington factor and mean opacities.
Trial estimates of an iterative procedure are given by the solutions at the previous time
step.
In earlier versions of the numerical code we imposed the condition that the transfer
equation solver and the moment equations solver should be iterated until convergence isachieved to a desirable extent. One thus obtains the radiation hydrodynamical quantities
which are consistent with the radiation eld at each time step. While this procedure works
well during calculations for the rst collapse (MMI), convergence is found to be much worse
after the second collapse due to complicated behaviors of the Eddington factors near theaccretion shock front. Therefore convergence between the transfer equation solver andthe moment equations solver is not imposed in the present version of the code, where
the transfer equation is solved once in several time steps under the constraint that thevariations of the Eddington factors should not be too large.
In Eulerian coordinates, the basic equations (eqs.#2.20]-#2.25]) are rewritten into
@Mr + 3vr2 @Mr = 0(2.33)@t@r3r = 43 @M(2.34)
3@r@v + v @v = ; GMr ; @p + gas
F F ; vQ 
(2.35)
2
@t @r
r
@r c

@ (e + E ) + v @ (e + E ) + (p + fE ) 3 @ (r2 v) ; (3f ; 1)E v = + ; " ; 3 @ (r2 F ) 
@t

@r

 @r3
r CR Q gd
@r3
(2.36)
1 ( @F + v @F ) + 1 @ (fqE ) = ; F F (2.37)
c2 @t@r q @rc33

CHAPTER 2. NUMERICAL SCHEME@e + v @e + p 3 @ (r2 v) = (c gasE ; 4
E@t @r  @r3

gas B ) +P

CR + Q ; "gd (2.38)

Equations (2.26) and (2.27) remain unchanged. For calculations in Eulerian coordinates,
variables to be solved are Mr   v T E F and Tgr(n) (n = 1 2 3) and the independent
variable is r. The computational procedure is similar to the Lagrangian case describedabove. Our computational policy how to choose either Lagrangian or Eulerian coordinatesin calculations is mentioned in x2.4 below.

2.3 The Energy Mean and Flux Mean of Opacities
In the transfer equation solver, the mean opacities ( E , F , and P ) are calculated as wellas the Eddington factor. We note here how to obtain J (r) and H (r), which are requiredfor calculating the mean opacities in equations (2.30)-(2.31).
MM suggested that one should obtain J (r) and H (r) by solving the monochromatic
moment equations of radiation using the current estimates of T (r) and f (r) (p.484 in
MM). One can, however, obtain J (r) and H (r) also directly by integrating I (r )
and I (r ) over . The former method, however, has a serious problem for non-grey
calculations as described below.
In a non-grey calculation the Planck function, B , increases exponentially with temperature at a frequency in the Wien regime, in contrast to the grey approximation, forwhich B is only proportional to T 4. Therefore, even a moderate contrast of temperature
throughout a computational region sometimes requires an extremely wide dynamic range
for I . For instance, an evolved protostar has a central core of T > 3000K with a coldenvelope of T 10K . Such a situation involves a great range of I over a hundred ordersof magnitude at  1014Hz throughout the computational region. It exceeds the available numerical range by far for ordinary computational tools. Thus, J ,H , and K in thecold region at high frequency should equal zero numerically. In such cases the Eddington
factor is not easily obtained because f  J =K is not well-dened numerically, although
the actual Eddington factor should behave regularly however the individual values of J
and K are close to zero. If the Eddington factor is not obtained, one cannot solve themoment equations.
In contrast to the method suggested by MM, we calculate the moments of radiation
directly by integrating I (r ) over  instead of solving the monochromatic moment
equations in the transfer equation solver. Our method works well even in the case whereJ and K are very small.
34

2.4. COORDINATES FOR COMPUTATION

2.4 Coordinates for Computation
2.4.1 Lagrangian or Eulerian
In the initial collapse phase, highly concentrated density structure is developed in the center and hence Lagrangian coordinates are suited to the calculations. Lagrangian scheme,
however, is ruled out when the sharp structure like a shock wave is generated in diuse areas where grid points are sparsely distributed. By this reason, calculations are performedon Lagrangian coordinates for the initial isothermal collapse, and the Eulerian takes theplace after the adiabatic core is formed|i.e., after the accretion shock is generated. Thismethod was found to successfully resolve the accretion shock at the surface of the rst
and second cores in our calculations. In early versions of our code, second-order accurate
advection (i.e., the van Leer advection) was adopted for Eulerian calculations, supplementing the Lagrangian RHD equations (2.20)-(2.25) (MMI). This method works wellfor calculations at the rst collapse stage but cannot be applied after the second collapse
occurs. This is because the advection solver suers from the CFL condition where theow is supersonic and timesteps become unacceptably small compared to the evolutionary timescale at late stages of protostellar collapse. Therefore the discretized equations in
the Eulerian formulation (eqs. #2.33]-#2.38]) have to be directly incorporated into the linearized matrix, instead of employing the external advection terms solver. Since upstreamdierentiation of rst-order accurate is adopted, this \direct Eulerian scheme" is inferiorin accuracy to the van Leer advection scheme. Nevertheless, it has a great advantagethat large CFL numbers are permitted and computational time is reduced to a reasonable
extent.

2.4.2 Grid Conguration
For discretization 400 radial grid cells and 64 frequency bins are used. The paths along
which the transfer equation is integrated are set tangentially to each radial grid (see Figure2.1). The innermost zone has 30 additional paths so as to provide a sucient number ofangle-mesh points for innermost rings of the radial grid.
A most important and cumbersome task in numerical calculations for protostar formation is how to resolve the radiative relaxation layer just ahead the accretion shock front at
the stellar surface(WN). Mesh points should be prepared in a suciently large number toresolve the thickness within  1 in the relaxation layer. Otherwise numerical diculties
can occur as experiences have shown, because the radiative cooling rate is not estimatedcorrectly in the relaxation layer.
35

CHAPTER 2. NUMERICAL SCHEME

Mesh Sizes

Initial Stage-1

Final Stage

-1.5

-2

-2.5

-4

-3

-2

-1

0
log R [AU]

1

2

3

4

Fig. 2.2.| Distributions of local mesh sizes is shown. The ordinate representslog(log ri+1 ; log ri) as a function of ri where i denotes the radial grid number. Opencircles denotes the initial conguration while lled triangles shows the nal distribution
after Eulerian scheme takes the place of Lagrangian scheme (see x2.4.1).
36

2.4. COORDINATES FOR COMPUTATION
Figure 2.2 delineates the distributions of local mesh sizes adopted in calculations
which are to be described in Chapter 5. The initial conguration (open circles) has alocal minimum between 1 and 10 AU, which is adjusted to be settled at the accretion
shock front after Lagrangian coordinates convey the grid to the position at which they
are switched into the Eulerian (lled triangles). Some trial calculations are necessary inorder to estimate the position of the shock front since one has to know a priori where ashock front stands.
Our scheme for determining the grid conguration may seem less sophisticated thanan adapted grid method, where grid points are gathered or dispersed automatically toachieve a desirable conguration to resolve the physical structure. Nevertheless, we prefer
our method because it is rm against numerical instabilities by adopting either Lagrangianor Eulerian scheme instead of freely moving coordinates. Any adaptive grid schemes are
not used in order to assure computational results to be free from unknown numericalinstabilities (see x1.1.2).

37

CHAPTER 2. NUMERICAL SCHEME

38

Chapter 3
INPUT PHYSICS
Various types of physical processes play signicant roles in the evolution of protostar
formation. Hydrogen is initially in molecular form but dissociates into atoms if the central
temperature reaches about 2000 K. Dissociation of hydrogen molecules triggers the second
collapse, and is hence essentially important to the dynamical evolution of protostellar
collapse. The ratio of specic heats for molecular hydrogen is 7/5 for suciently hightemperatures, but it varies depending on T for lower temperatures due to discretenessof the rotational transitions of H2 molecules. The partial degeneracy of electrons begins
to dominate the gas pressure when the central density arrives at 1g  cm;3 , and theequation of state (EOS) for ideal gas should be replaced by the non-ideal EOS. Opacities
due to dust grains, molecules, and atoms should be given for solving the radiative transfer
equation from the stellar interior to the diuse envelope of infalling gas and dust.
In this chapter we introduce physical processes considered in calculations. Chemical
abundances(x3.1), specic heat for molecular hydrogen(x3.2), equation of state(x3.3), andopacities(x3.4) are described.

39

CHAPTER 3. INPUT PHYSICS

3.1 Chemical Abundances
3.1.1 General Remarks
Elements that we consider in calculations consist of H2, H, H+, He, He+ , He++, and asmall fraction of heavy elements. The mass fractions of hydrogen, helium, and heavy
elements are supposed to be X = 0:70, Y = 0:28, and Z = 0:02. Dissociation degree ofhydrogen molecules and ionization degrees of hydrogen and helium are evaluated assuming
the chemical equilibrium under given temperature and density at each timestep. This
assumption is well justied in our calculations because both dissociation and ionization
becomes important after the density is so high that the timescale of chemical reactions is
negligibly short.
Hence we evaluate dissociation and ionization degrees by solving classical Saha equations, shown in x3.1.2 below. Classical Saha equations are applicable where free particles
obey the Maxwell-Boltzmann distribution, but are ruled out where the degeneracy ofelectrons is not negligible. Furthermore, classical Saha equations should not be used also
where the density is so high ( > 1g  cm;3 ) that the pressure ionization dominates thethermal ionization. We abandon classical Saha equations when the density reaches a critical value (typically  = 10;4-10;3g  cm;3 ) and instead use the tabulated data provided
by Saumon, Chablier, & Van Horn(1995) (x3.3).
Actual molecular clouds have a small but nite fraction of electrons (typically ne=nH;
10 7 for nH  104cc;1 , where nH = n(H ) + 2n(H2)) despite extremely low temperatures( 10K) because cosmic rays ionize molecules and atoms (Nakano 1984 and references
therein). The ionization degree in molecular clouds plays a signicant role when one
considers the coupling between gas and magnetic elds, but it has no importance fordynamical evolution once magnetic elds cannot support cloud cores against the selfgravity. While the heating by the non-thermal ionization by cosmic rays is taken into
account as a uniformly embedded heat source(Chap. 2), we do not solve explicitly theionization degree by cosmic rays in molecular clouds.

3.1.2 Equations for Dissociation and Ionization Degrees
Here classical Saha equations for evaluating dissociation and ionization degrees are derived.
The equilibrium state of the following type of reactions,A*) B + C40

3.1. CHEMICAL ABUNDANCES
are described as (Saha equation)nB nC = ZB ZC exp(; )
(3.1)
nA
ZA
kT
where ZA, ZB , and ZC are the partition functions for particles A, B , and C , respectively,
and is the energy potential for the reaction.
First we consider dissociation of hydrogen molecules. The partition functions are
written down as
3=2
(3.2)
ZH = geHgNH (2 mHh3kT ) 
3=2
ZH2 = 12 geH2 gNH2 (2 2mhH3kT ) #1 ; exp(; Tvib )];1 T 
(3.3)
rot
where ge and gN represent the statistical weights for electrons in bound states and nuclear
spin of hydrogen, respectively, and vib  6100K and rot  84:3K are the energy deposits
for the lowest transitions of vibration and rotation, respectively, of molecular hydrogen(x3.2). Simple consideration leads geH = 2, gNH = 2, geH2 = 1, and gNH2 = 4.
The number densities for hydrogen atoms and molecules, nH and nH2 , are related withthe mass density as
nH = yX m HnH2 = (1 ; y)X 2m 
H
where the dissociation degree, y, is dened byy = (H)+(H)(H ) :2
Combining the above equations, one obtains Saha equation for hydrogen dissociation fromequation (3.1):

y2 = 4mH ( mHkT )3=2 #1 ; exp(; vib )] rot exp(; d )(3.4)
1 ; y Xh3
T T
kT
where the dissociation potential d is 4.48 eV. This is a quadratic equation with respect
to y and is easily solved under given T and . Equation (3.4) is equivalent to equation(A.15) in Nakano, Ohyama, & Hayashi(1968). Black & Bodenheimer(1975) derived an
approximate formulation for equation (3.4).
Saha equations for ionization equilibrium are also derived similarly. The partition
functions for hydrogen and helium ions and electrons are)3=2 
ZH+ = gNH (2 mHhkT
(3.5)
3
41

CHAPTER 3. INPUT PHYSICS
3=2) 
ZHe = geHe gNHe (2 mHhkT
3
+

+

3=2

ZHe = gNHe (2 mHh3kT ) 
++

3=2

) :
Ze = ge (2 mhe kT
3

(3.6)
(3.7)
(3.8)

The statistical weights are found to be gNHe = ge = geHe+ = 2. The ionization degrees aredened by(H+ )x1 = (H+ ) + (H) 
for H+,+)
x2 = (He+(He) + (He)for He+ ,(He++ )x3 = (He++) + (He+ )for He++, and hencenH = (1 ; x1 )yX m HnH+ = x1 yX m 
HnHe = (1 ; x2 )Y 4m 
H
nHe+ = x2 (1 ; x3 )Y 4m HnHe++ = x2 x3 Y 4m :H
Using the above relations Saha equations for ionization are derived asx1 = emH (2 mekT )3=2 exp(; 1 )(3.9)
1 ; x1h3kT
x2 (1 ; x3 ) = 4 emH (2 mekT )3=2 exp(; 2 )(3.10)
1 ; x2h3kT
x3 = emH (2 mekT )3=2 exp(; 3 )(3.11)
1 ; x3

h3
kT
where
;e 1  x1 yX + x2 (1 + x3 ) Y4 :The ionization potential for H+, He+ , and He++ are 1 = 13:60eV, 2 = 24:58eV, and3 = 54:40eV.
42

3.1. CHEMICAL ABUNDANCES
Since equations (3.9)-(3.11) are coupled with each other through electron densities,
they should be solved simultaneously. Therefore the linearized system of equations (3.9)(3.11) is solved by Gaussian elimination of a 33 matrix, and this procedure is iterated
according to Newton-Raphson method until the solutions converge.
Figure 3.1 presents solutions of Saha equations (eqs. #3.4] and #3.9]-#3.11]). Dissociation of molecular hydrogen is seen to begin near T 2000K for  10;7, whichcorrespond to the central temperature and density when the second collapse is triggered.
Ionization of H and He occurs around T 104K and He+ begins to ionize when temperature reaches about 3  104K. The temperature at which ionization begins rises as
the density increases, but the pressure ionization takes place when the density arrives at
 1g  cm;3 and atoms are enforced to ionize for higher densities regardless of temperature. Classical Saha equations are ruled out there(x3.3).

43

CHAPTER 3. INPUT PHYSICS

1

1
(a)

(b)
.8

.6

.6

.4

.4

.2

.2

0

0

y

.8

3

3.5

4

4.5

5

1

3

3.5

4

4.5

5

3.5

4

4.5

5

1
(c)

(d)

.8

.8

.6

.6

.4

.4

.2

.2

0

0
3

3.5

4

4.5

5

3

Fig. 3.1.| Solutions of Saha equations for dissociation degree of molecular hydrogen(panel (a)) and ionization degrees of hydrogen, (b), and helium, (c) and (d) for He andHe+ , respectively.
44

3.2. SPECIFIC HEAT FOR MOLECULAR HYDROGEN

3.2 Specic Heat for Molecular Hydrogen
The specic heat at constant volume is described as@ 2 (T ln Z )CV = mk T @T
(3.12)
2
H
where Z is the partition function. The specic heat consists of the vibrational androtational parts in addition to the translational part.

CV = CV tra + CV vib + CV rot:
The translational part in the specic heat is simplyk :CV tra = 23 m
H
The partition function for molecular vibration, approximated by a harmonic oscillator, isgiven by1X
Z = expf;h (n + 12 )=kT g:
n=0
Equation (3.12) has an analytical expression for the vibration:

CVH2vib = mk ( Tvib )2
H

exp( vib
T ) 

vib
#exp( T ) ; 1]2

wherevib  h
k = 6100K:
The partition function for rotational energy levels are

Zo =for ortho-H2 andfor para-H2, where

Zp =

1
X


(2j + 1) exp#;j (j + 1) rot ]
T
j =135
1X(2j + 1) exp#;j (j + 1) rot ]
T
j =0242
rot  8 h2Ik = 85:4K:
Here I is the moment of inertia of a hydrogen molecule.
A principle in quantum mechanics demands the wave functions describing the nuclear
spin and rotation of H2 to be asymmetric with respect to substitution of the nuclei. The
wave function for the rotation is symmetric for even j (i.e., para-H2) or asymmetric for odd

45CHAPTER 3. INPUT PHYSICS
j (i.e., ortho-H2), while the total nuclear spin of H2 has three independent states which issymmetric and a single state which is asymmetric. Consequently, ortho-H2 corresponds
to the \triplet" and para-H2 to the \singlet", and then the ortho:para ratio would beclose to 3:1 if they are in thermal equilibrium. Equilibrium states, however, are notachieved generally in astronomical sites. Especially, the abundance ratio of ortho- to paraH2 is highly uncertain in molecular clouds. While rotational transitions by radiation ishighly improbable for H2 in molecular clouds(Raich & Good 1964), thermal equilibrium isattained over a timescale of 106 or 107 yr in molecular clouds through the proton exchange
reaction, which indicates that actual ratio of ortho-H2 to para-H2 in a given cloud would
depend on its chemical history and neither a normal nor a thermal equilibrium ratio
may be safely assumed (Flower & Watt 1984). Observations of the ortho:para ratio in
molecular clouds are quite dicult since molecular hydrogen does not emit the dipole
radiation, so that no information is available about the physical state of H2 .
Therefore we consider a couple of extreme conditions for the abundance ratio, i.e.,in equilibrium and in completely non-equilibrium. When ortho-H2 and para-H2 are inthermal equilibrium, the rotational part in the specic heat is determined by@ 2 #T ln(3Z + Z )]:CV rot = mk T @Top2
HOn the other hand, in case that there is no interaction between ortho-H2 and para-H2,
2@k
CV rot = m T @T 2 #T (fo ln Zo + (1 ; fo ) ln Zp )]:H
Here fo is the fractional abundance of ortho-H2.
Figure 3.2 delineates the specic heat at constant volume as a function of temperaturefor each case for the equilibrium and non-equilibrium cases with those of pure ortho-H2
and para-H2 for reference. In the non-equilibrium case, where the ortho:para ratio is
uncertain as mentioned above, we assume fo = 3=4 here as an example.
Some general tendencies are observed: (1) neither the rotational nor the vibrationalenergy levels are excited at low temperatures near T = 10K (therefore CV = 3k=2mH)
but (2) the rotational levels begin to be excited as T increases, which brings CV up to
5k=2mH for T > 500K, and (3) further increase of temperature excites the vibrational
levels and CV rises again toward 7k=2mH.
Divergence of the CV curves around T 100K is due to discreteness of a few lowest
rotational energy levels of molecular hydrogen. Only the equilibrium ortho-para ratio
permits the collisional transition of J = 1 ; 0, which yields a striking peak at T 50K.
Another peak at T 170K corresponds to the J = 0 ; 2 transition and is observed only46

3.2. SPECIFIC HEAT FOR MOLECULAR HYDROGEN

4

3.5

3

2.5

in equil.

2

in non-eq.
pure ortho

1.5pure para
1
1

1.5

2

2.5
log T [K]

3

3.5

4

Fig. 3.2.| Specic heat at constant volume for molecular hydrogen is depicted as afunction of temperature. Ortho- and para-H2 are assumed to be in equilibrium (dashedline) or not in equilibrium (solid line). Cases for pure ortho- or para-H2 are shown aswell for reference (short-dashed and long-dashed lines, respectively). Ortho-para ratio
is supposed to be 3:1 as an example in the non-equilibrium case. Actual ratio is quite
uncertain in molecular clouds (see text).
47

CHAPTER 3. INPUT PHYSICS
in case of the pure para-H2. These dierences disappear when temperature becomes sohigh that the classical treatment for the rotational energy levels is valid.
In our calculations we adopt the non-equilibrium case. The non-equilibrium caseis most plausible because the thermal equilibrium is highly improbable to be achieved
within the free-fall timescale of a molecular cloud core. Possible eects on dynamical
evolutions due to the dierence of the CV curves are negligibly small because it occurs for
T < 100K, where the isothermal collapse has not yet been decelerated signicantly andthus the evolution proceeds insensitively to the ratio of specic heats.

48

3.3. EQUATION OF STATE

3.3 Equation of State
3.3.1 EOS for Ideal Gas
The internal energy, e, is the sum of molecular and atomic hydrogen, helium atoms,
dissociation and ionization energies, i.e.,

e = eH2 + eH + eHe + edis + eion :Note that eH and eHe include contributions by ions and electrons. For ideal gas, eachcomponent of the internal energy is expressed as

X C H2 T
eH2 = (1 ; y) 2mVH

eH = 32 y(1 + x1 ) mX kTH

Y kTeHe = 23 (1 + x2 + x2 x3 ) 4m
edis = yX 2mk d HHeion = #x1 yX1 + x2 f(1 ; x3 )2 + x3 3 g Y4 ] mk 

and the mean molecular weight, , is given byH

2x1y X + 1 + x2 + x2 x3 Y:
;1 = 1 + y +2
4

The dissociation and ionization degrees, y x1 x2  and x3, are discussed in x3.1 above, andthe specic heat at constant volume for molecular hydrogen, CVH2 , in x3.2 above. Equationof state (EOS) for ideal gas is written simply ask T
p = mH(3.13)

where p is the gas pressure.
For low densities (typically  < 10;4g  cm;3 ) the ideal gas approximation is valid and
we use equation (3.13) for the equation of state. However, non-ideal eects take place forhigher densities, where both the equation of sate for ideal gas (eq. #3.13]) and classical
Saha equations (eqs. #3.4], #3.9]-#3.11]) are inapplicable. The equation of state for thehigh density regime is described in the next subsection.
49

CHAPTER 3. INPUT PHYSICS

3.3.2 EOS for Non-ideal Gas
We adopt tabulated data obtained by Saumon, Chablier, & Van Horn(1995) for equationof state in the high density regime. Non-ideal eects taken into account by Saumon et al.
are summarized as follows.
1. Pressure ionization.
Pressure ionization should be taken into account when e2 =a0 < F , where a0 is
the Bohr radius and F is the zero-temperature Fermi energy of an electron. Thiscondition corresponds to  > 1g  cm;3 and hence pressure ionization dominates inthe stellar interior.
2. Degeneracy of electrons.
The pressure of degenerate electrons dominates the thermal pressure when kT < F .
This condition is satised in the center of a protostar at a later stage of the evolution,and degenerate gas cause the temperature inversion in the stellar interior as seen
later.
3. Non-ideal Coulomb eects.
When the electrostatic potential energy between two protons exceeds their kinetic
energy, i.e., kT < e2 =a, where a is the mean interparticle distance, non-ideal eects
between particles through the Coulomb interaction play an important role.
The tabulated data are calculated for pure hydrogen and pure helium, and then interpolation in composition between them is required for their mixtures. However, ionization
degrees in H and He mixture are related with each other through the electron density andno interpolation schemes would provide exact estimates for the mixture. Nevertheless,
we adopt interpolated values for the mixture as the best estimates available for non-ideal
EOS to date.

50

3.4. OPACITIES

3.4 Opacities

3.4.1 Dust Opacities
The opacities we adopt for dust grains are provided by Yorke(1998). The tabulated data
are based on the dust model of Preibisch et al.(1993). The opacity tables contain thethree components of amorphous carbon, silicate, and dirty ice-coated silicate. The dust
properties are tabulated in Table 3.1. The grain radius for each component is assumed tobe constant for brevity.
Figure 3.3 shows the frequency-dependence of the extinction opacities and albedos.
The value of  , which is dened by  /   in the submillimeter region (c.f., Beckwith
et al. 1990), is 2 for silicate and ice-coated silicate, and is 1 for amorphous carbon.
Albedos are almost zero for  < 1013 Hz but they are close to unity for higher frequencies.
This means that the eects of scattering are important where the radiation eld is dominated by near-infrared and optical photons. In order to show the temperature dependenceof the dust opacities, their Rosseland mean is depicted in Figure 3.4. The opacity jumps
at T = 125K 1500K and 2000K correspond to sublimation of the icy mantle, silicate,
and amorphous carbon, respectively.
For T > 2000K all the dust components are destroyed, and then the gas opacities (i.e.,
molecular and atomic opacities) take the place. The gas opacities are introduced in thefollowing subsection.

3.4.2 Gas Opacities
For the gas opacities we adopt Alexander & Ferguson(1994) for log T #K] < 3:75 andIglesias & Rogers(1996), or updated OPAL, for higher temperatures. The opacities asfunctions of temperature are shown in Figure 3.5 under some xed values of R, whereR  #g  cm;3 ]=(T=106#K]).
Opacity data for T  2000K, where the dust opacities by Preibisch et al. (1993) areused (x3.4.1), are abandoned in our calculations in order to avoid double count of thedust opacities. (Alexander & Ferguson also takes into account the dust opacities, whichcorrespond to at curves for T < 1500K.) The sharp drop at T 2500K is due to destruction of H2O, which makes a large contribution to the opacities for lower temperatures.
Hydrogenic absorption (principally, H; and Rayleigh scattering) dominates the opacities
for higher temperatures.
The bound-bound and bound-free absorption of hydrogen raise the opacities dramatically as temperature increases over 104K, and the opacities gradually decrease for higher
temperature where the free-free absorption dominates. For lower values of R, that is, for51

CHAPTER 3. INPUT PHYSICS

4

2

0

-2

-4

11

12

13

14

15

11

12

13

14

15

1

.8

.6

.4

.2

0

Fig. 3.3.| The dust opacities based on data by Preibisch et al.(1993) as functions offrequency are depicted. The upper panel shows the extinction opacities for each component of amorphous carbon (solid line), silicate (dotted line), and dirty ice-coated silicate(dashed line). In the lower panel the albedos are delineated (denotations for lines are thesame as the upper panel).
52

3.4. OPACITIES
lower densities at a xed temperature, the opacity is retained at a constant value whichcorresponds to the Thomson cross section for T > 104K.
Tabulated data of the gas opacities are given in the form of the Rosseland mean. Hence
we encounter a problem how to incorporate the gas opacities into the monochromaticradiative transfer equation. We choose a simple resolution for this problem that the
monochromatic opacities for gas is constant with respect to frequency, i.e., gas
 (T ) gas
R (T ). This treatment for  is physically incorrect but it would cause no errors inpractice for the problems of our present interest because of the following reason. Solutionsof the monochromatic radiative transfer equation are used to obtain fE , E , and F . The
gas opacities dominates in general in the protostellar interior, where the optical depth isso large at every frequency that fE is always 1=3 and E and F equal R independently of
the frequency dependence of opacities. Although the opacity gap between the protostellar
surface and the dust sublimation front can be optically thin, possible errors in fE , E , andF would not be large because the spectrum of the radiation eld in the opacity gap does
not signicantly depart from the Planck function of the local temperature. Furthermore,shape of SEDs would not be aected by the simplication on the gas opacities since SEDs
of YSOs are contributed mostly by the dust envelope and the blackbody radiation fromthe protostellar surface, and then SEDs are independent of the frequency-dependence ofthe gas opacity.

53

CHAPTER 3. INPUT PHYSICS

2

1

0

-1

-2

-3
1

1.5

2

2.5

3

3.5

Fig. 3.4.| The Rosseland mean dust opacity is shown. Ice mantle is vapored at T =125K, and silicate and amorphous carbon are sublimated at T = 1500K and T = 2000K,
respectively.

54

3.4. OPACITIES

6
1
4

-1

2

-3
-5

0
log R = -7
-2

-4

-6
3

3.5

4

4.5

5

5.5

6

6.5

Fig. 3.5.| The gas opacities are presented. The data are taken from Alexander andFerguson(1994) for log T #K] < 3:5 (triangles) and from Iglesias and Rogers (1996) for
higher temperatures (squares). Each curve corresponds to dierent R, where R = #g 
cm;3 ]=(T=106#K])3.
55

CHAPTER 3. INPUT PHYSICS

Table 3.1. Dust Properties.
Dust Components

Grain Radiusa

Number Densityb

Dirty ice-coated silicate
Silicate
Amorphous carbon

5:998  10;6
5:236  10;6
9:462  10;7

1:652  1012= ;=1:652  1012
2:438  1014=2:623  1014

Sublimation
Temperature
125 K
1500 K
2000 Ka Radius of a dust grain in #cm].
b Number density of dust per gram gas where the ice-mantle survives (left) or other-wise(right).

56

Chapter 4
THE FIRST COLLAPSE
The aim of this chapter is to realize typical properties of the earliest stage in star formation, i.e., the rst collapse. Physical processes put into calculations are simplied here
in order to extract essential features from computational results. Ratio of specic heatsis xed at 5/3, and dissociation of molecules and ionization of atoms are not taken into
account (and therefore evolutions after the second collapse are not considered in thischapter). We parameterize dust opacities by supposing power-law distributions. Theseassumptions may modify the results quantitatively by a small factor but they do not cause
qualitative changes. Under these simplications, we nd that the conditions for termination of isothermality and the physical properties of the rst core are in good agreement
between numerical results and analytical estimates.
In particular, the conditions for termination of isothermality derived by us have great
importance because they are found to indicate that the conventional idea that opaquenessbreaks isothermality is invalid. Appendix C is devoted to derivations of the criteria forviolation of isothermality and to discussions on problems that would arise from these
newly derived criteria.

57

CHAPTER 4. THE FIRST COLLAPSE

4.1 Assumptions

4.1.1 Initial Conditions
For the total mass of the cloud core Mcc , we chose Mcc = 1M, 0:1M,0:3M , and
3M . The initial density and the radius are adjusted so that the cloud core should beslightly more massive than the Jeans mass for each case. As is shown in Table 4.1, theinitial density is init 10;20-10;17g  cm;3 (i.e., nH2 104-107cm;3 ) and the radius is
Rcc 1000-30000AU , which is compatible with the typical scales of observed molecular
cloud cores. The initial density distribution is assumed to be homogeneous except a casein which the cloud core is initially in hydrostatic equilibrium (see x4.3.2). In every model,
the velocity of each uid element is set to zero in the initial state.
The initial temperature is also an important model parameter, Tinit =10K, 5K, and
30K. Dierent initial temperatures account for possible variations of the environments
in molecular clouds. Dependence of protostellar evolution on initial temperature will be
discussed in x4.3.3.

4.1.2 Boundary Conditions
The inner and outer boundaries are given by the center and the edge of the cloud core,respectively. The outer boundary is xed through evolution at R = Rcc shown in Table
4.1. Boundary conditions for radiation involve a brief discussion as follows.
Boundary conditions for the transfer equations (2.11) and (2.12) (see Appendix 2.2)
are given by h on the p-axis shown in Figure 2.1 and at the cloud edge. On the p-axish = 0 because of symmetry, and h at the cloud edge is determined by given incident
radiation I ;(r = Rcc ). In this chapter, two dierent boundary conditions for radiation
are considered, which provide a couple of extreme cases.

Case 1. There is no incident radiation but there is a heat source CR to compensate

for the radiative energy loss in the initial state. The heat source represents theheating by cosmic rays and photons from surrounding stars, which support thetemperature of the molecular cloud against cooling by the thermal emission of dust
grains (Hayashi & Nakano 1965 Goldsmith & Langer 1978).

Case 2. The cloud core is surrounded by isotropic blackbody radiation of T = Tinit.
Though Case 2 seems to keep the cloud core naturally at a constant initial temperature,
ambient blackbody radiation of Tinit is too idealized a condition for actual molecular
cloud cores. Furthermore, Case 1 is more favorable rather than Case 2 for the purpose58

4.1. ASSUMPTIONS
of calculating emergent SEDs because in Case 2 incident radiation modies the emergent
spectra. In particular, in the initial state for Case 2 the net ux equals zero and emergent
SEDs cannot be calculated. Therefore we adopt Case 1 for all the models except whereindicated.
Boundary conditions of radiation moment equations are given as radiation uxes at the
inner and the outer boundaries. The inner ux at the center is set equal to zero becauseof symmetry. The outer boundary condition is that the outer ux FI +1 is eliminated usingh, where h is dened by FI +1 = hcEI +1 and is given by the transfer equation solver.

4.1.3 Opacity
We assume that the opacity is due to the absorption by dust grains, neglecting the scattering in this chapter, and that the gas temperature couples closely with the dust temperature. These assumptions are known to be well justied in molecular cloud cores andin very young protostars. (Scattering by dust grains should be taken into account todetermine the temperature structure of envelopes for more evolved young stellar objects
like class I and early class II.) We have referred to Adams & Shu(1986) for frequencydependence of the dust opacity in higher frequency than far-infrared except for the linefeatures, which have been removed in our models for simplicity. Dust grains are assumed
to consist of silicate and graphite, where H2O ice is excluded since it has no signicantinuences on the dynamical evolutions. The opacity for lower frequency is characterizedwith a parameter  , where  = 0(=0 ) , taking the ambiguity of the dust opacity into
account. Three dierent models are considered with respect to  .

(a).  = 2:5 for 4:6  1011 <  < 3  1012 Hz and  = 1:5 for  < 4:6  1011Hz, taking

a recent result into account (Pollack et al. 1994). (This model has a comparable
value for P (10K) to that of the model of  = 2.)(b).  = 1:5 for  > 3  1012Hz.
(c).  = 0:1 Model(a).

In Models (a) and (b), 0 is determined to be connected at 0 = 3  1012Hz with Adams& Shu's opacity, that is, 31012 Hz = 0:3cm2  g;1 . Note that the absorption coecients
 are given with the mass opacity of dust  by  (r) =  (r). We shall examine in
x4.3.4 how the dierence among Models (a)-(c) aects the evolution.
We consider the vaporization of dust grains at Tvap  1200K . Though the opacity forhigher temperature should be taken place by the molecular opacities, we simply reducethe dust opacity to nearly (but not exactly) zero for T > Tvap . In order to assure the59CHAPTER 4. THE FIRST COLLAPSE
convergence of the Eddington factor between the transfer equation solver and the moment
equations solver, the dust grains are supposed to vanish completely at T = Tvap + %T ,where %T is taken to be 100 degrees, and to decrease linearly with increase of temperature
for Tvap < T < Tvap + %T .
The dust temperature is supposed to be equal to the gas temperature in this chapter.

4.1.4 Equation of StateWe adopted the equation of state of ideal gas for simplicity in this chapter.
e =  ;1 1 mk T
(4.1)
Hwhere  k  and mH are the ratio of specic heats, the Boltzmann constant, the mean
molecular weight and the mass of hydrogen atom, respectively. In this chapter,  is
xed at 5/3 except a case where  equals 7/5 instead, and  is assumed to be constantat 2.46, which corresponds to the same chemical abundance as in L69. More realistic
models should include the eects of rotational and vibrational degrees of freedom for H2
molecules, and of dissociation and ionization of the molecules. The former modies  from
5/3 and the latter changes the value of  as functions of the temperature. Moreover, thedissociation and ionization energies should be taken into account in the energy equation.
All of these eects will be considered in the subsequent chapter.

60

4.2. RESULTS

4.2 Results
4.2.1 Overview for a Typical Case
Results of the typical case(i.e., Model M1a in Table 4.1) are overviewed in this section.
In Model M1a, the mass and the initial temperature are taken to be 1M and 10K,
respectively, with the opacity of Model (a) in x4.1.3. First, we present a brief summaryof the results in Figures 4.1 and 4.2.
The initial state is indicated by the thick line in each panel. At early stages, thecloud core collapses isothermally, which corresponds to the rst two lines in each panel.
The density distribution at large radii develops an r;2 law reminiscent of the isothermal
self-similar solution (Larson-Penston solution, c.f., x1.1.1). This isothermal stage lasts fora free-fall time, i.e.,sR3 = 1:77  105yrt = 2 2GM(4.2)
until the central temperature departs signicantly from the initial value (dotted lines).
The rapid increase of the central temperature to T 2000K forms a central adiabaticcore|that is, the rst core. The drop in the density distribution at a few AU showsthe eective radius of the rst core, whose position gradually moves inward throughthe rst collapse stage. The super-sonically infalling gas in the envelope is decelerated
and thermalized at the surface of the rst core, where the shock front is observed as a
sharp peak in the velocity distribution. In the radiative ux distribution, a few lines justbelow the dashed line have small humps at the radius of the rst core, which indicatesthat the thermal energy generated by the shock heating is immediately radiated away.
The distribution of change-rate of the specic entropy shows the negative peak on thesurface of the rst core, where the entropy is taken out by radiation. In terms of themass distribution, the inner region has a r3 slope which corresponds to the homogeneousdensity structure inside the rst core, while the outer region is proportional to radius
which resembles the Larson-Penston solution in the infalling envelope (i.e.,  / r;2).
Therefore, the plateau between these two slopes implies the mass of the rst core, whichreads 0:05M in the nal stage of this calculation. This value is consistent with the result
by Boss & Myhill(1992). Figure 4.2 shows more directly that the mass of the rst core,
which is indicated by the peak in the velocity distribution and the drop in the density
distribution, is increasing up to 0:05M.
The dashed lines represent the stage when the dust grains begin to evaporate in thecentral region. The complicated behavior in the radiative ux for R < 1AU is caused by
the sudden decrease of the opacity due to the dust vaporization.
61

CHAPTER 4. THE FIRST COLLAPSE

(a)

(b)
-10

1000

log T [K]

-12

-14
100
-16

-18

10
-2

0
2
log R [AU]

4

-2

0
2
log R [AU]

4

0
2.5
(c)

(d)
-2

2

1.5

-4

1

-6

.5

-8

0
-10
-2

0
2
log R [AU]

4

-2

0
2
log R [AU]4

Fig. 4.1.| Evolutionary sequences for a typical case (Model M1a) as functions of radius.
Dierent lines in each panel represent the dierent elapsed times tabulated in Table 4.2.
The initial state is specied by thick solid lines. See text for the meanings of dotted anddashed lines. The gure depicts (a)temperature, (b)density, (c)infall velocity, (d)mass,(e)radiation temperature, (f)radiative ux, (g)Planck mean optical depth, and (h)rate of62change of the specic entropy.

4.2. RESULTS

1000

4

(e)

(f)

2

100

0

-2
10
-4

-6
1
-2

0
2
log R [AU]

4

-2

0
2
log R [AU]

4

0
(g)

4

(h)
-.02

2
-.04
0
-.06
-2
-.08
-4
-.1
-6

-.12

-8

-.14
-2

0
2
log R [AU]

4

-2

63

0
2
log R [AU]4CHAPTER 4. THE FIRST COLLAPSE

(a)

(b)
-10

1000

log T [K]

-12

-14
100
-16

-18

10
-8

-6

-4

-2

0

-8

-6

-4

-2

0

2.5
(c)

(d)

4

2
2

0

1.5

-2
1
-4
.5
-6

0

-8
-8

-6

-4

-2

0

-8

-6

-4

-2

0

Fig. 4.2.| Evolutionary sequences for a typical case (Model M1a) as functions of mass.
See also the caption of Fig. 4.1. The gure depicts (a)temperature, (b)density, (c)infall
velocity, and (d)Planck mean optical depth.
644.2. RESULTS

5

5
R < 1 AUR < 10 AU

0

0

-5

-5

-10
11

11.5

12

12.5

13

13.5

-10
11

5

5

0

0

-5

-5

-10
11

11.5

12

12.5

13

13.5

-10
11

11.5

12

12.5

13

13.5

11.5

12

12.5

13

13.5

11.5

12

12.5

13

13.5

5
4

2

0

0
-5

-2
-10
11

11.5

12

12.5

13

13.5

11

Fig. 4.3.| Evolution of spectral energy distributions (SEDs) for Model M1a. Panels
except the lowest right show SEDs integrated over the dierent beamsizes of an imaginary
telescope. A beamsize corresponds to a concentric circle of the radii indicated in thepanel. The lowest right panel shows the radius at which the optical depth measured froman observer reaches unity as a function of frequency.
65

CHAPTER 4. THE FIRST COLLAPSEThe evolution of the spectral energy distribution (SED) is presented in Figure 4.3.
The panels except the lowest right suppose observations with dierent beamsizes of thetelescopes. These panels show that evolutionary changes would be detected more clearlyby observations with smaller beamsizes. This is because a ner beam can concentrate theeld of view on a more restricted region in the center of the collapsing cloud core, where theevolution is proceeding most dramatically. The spectral evolution for ne beams is mainly
due to the developing density distribution, and is not due to the temperature evolution.
Indeed, the SEDs give us little information about the rapid increase of temperature in
the rst core. The lowest right panel in Figure 4.3 gives us an answer for the reason why
SEDs are not aware of the existence of the hot central core. The rst core of T 2000K
radiates mostly in the frequency range near  1013Hz, at which R( = 1) > 100AU,
while the size of the rst core is only a few AU. An observer therefore cannot see throughdown to the surface of the rst core, but only observes the cold envelope of T 20K.
Second, we examine the physical processes which determine the thermal evolution fromthe isothermal collapse to the adiabatic contraction. Figure 4.4 illustrates the thermal
evolution at the center of the collapsing cloud core for the same model as in above (i.e.,
M1a in Table 4.1). In Figure 4.4 the evolution proceeds from left to right as the central
density increases monotonically with time.
Initially the radiative cooling rate, "th  4 P B , and the cosmic ray heating rate,
CR , are balanced, while the radiative heating rate, c E E , and the gas compressional
heating rate, ;g  (p=)r  v, are considerably smaller as compared to them. Here theheating and cooling rates are evaluated per unit mass. As the cloud core collapses, ;gshows an increase proportional to the square root of the density for c > 10;17g  cm;3 ,which is accounted for in equation (C.3). The slight disturbance of the compressional
heating rate appearing in c < 10;17g  cm;3 is caused by the homogeneous and static
initial condition, which is far from the similarity solution. The radiative heating rate also
increases with the central density, reecting that the dust grains (the emitter of photons)
are being concentrated into the central region.
In principle, the increase of temperature due to the compressional heating can raise theradiative cooling rate because "th has a positive dependence on T . For c < 10;14g  cm;3,
however, the radiative cooling rate is much larger than the compressional and radiative
heating rates, and "th is a very steep function of temperature, so that a very slight increaseof temperature is enough to compensate for the compressional heating. Thus the gas stays
isothermal for c < 10;14g  cm;3 .
When ;g gets suciently close to "th at c 10;14g  cm;3 , the radiative coolingno longer keeps the cloud core isothermal against the compressional heating. Thus, the66

4.2. RESULTS

15(a)cosmic ray heatingradiative heating10

radiative coolinggas compression
5

0

-5
-18

-16

-14

-12

-10

-16

-14

-12

-10

3.5
(b)
3

2.5

2

1.5

1
-18

Fig. 4.4.| Thermal evolution at the center of a collapsing cloud core for Model M1a.
(a) The heating and cooling rates per unit mass as functions of the central density. The
evolution proceeds from left to right. (b) log T -log  plane.
67CHAPTER 4. THE FIRST COLLAPSE
central temperature Tc departs from Tinit. At c 10;12g  cm;3 the radiative heating
rate overtakes "th , when the optical depth to the thermal radiation reaches the orderof unity. (Note that the condition \the radiative heating ratethe radiative coolingrate" corresponds to  > 1 only when the radiation eld has no external contributions.)
Further evolution occurs without ecient cooling and thus enters the adiabatic phase.
The central temperature increases linearly in the log Tc -log c plane with the gradient of
2=3 for c > 10;12g  cm;3 . This behavior proves that the isothermal collapse has
been replaced by the adiabatic contraction, because in adiabatic evolutions d log T=d log equals  ; 1 with  ; 1 = 2=3 for  = 5=3. The adiabatic evolution decelerates theself-gravitational collapse and thus diminishes the value of d ln ;g =d ln c to be less than
1=2. The abrupt drops of the radiative heating and cooling rates at c 10;9 are due tothe dust vaporization.
It can be shown that our results are consistent with past works. In L69 the \initial"
mass and radius of the rst core when the central density reaches 170K are 1031g and6  1031g  cm;3 , respectively, which are roughly equal to our results at a corresponding
evolutionary stage (see Fig. 4.1). Figure 4.4 shows good agreement with Figure 3 inBoss(1988), who performed three-dimensional RHD calculations with the 3D Eddington
approximation. Moreover, his result for the lower opacity model is also quite similar to
our lower opacity model (Model M1c, see x4.3.4 and Figure 4.9 below). These imply that
our results in spherical symmetry would be helpful for investigations in general symmetry.

4.2.2 Characteristic Densities
Appendix C shows that isothermal evolution of gravitationally collapsing clouds is terminated when the central density of the cloud reaches a critical density, crit, dened
byTinit )6
=; = 1:4  10;14g  cm;3 ( 1:6  10;2Pcm2  g;1 )2( 10
(4.3)
K
or
Tinit )4=3 :dif MD = 1:2  10;13g  cm;3 ( 1:6 P10;2 );2=3 ( 10K(4.4)
Equation (4.3) is derived from the condition that the compressional heating rate of gas
overtakes the thermal cooling rate of dust grains. On the other hand, equation (4.4) comesfrom the condition that the compressional heating rate reaches the energy transport rateby the radiative diusion. If the central density of a collapsing cloud arrives at =; before
the optical depth of the cloud reaches unity, crit equals =; . Otherwise crit is equal todif MD. For reference, the central density when the optical depth of the collapsing cloud68

4.2. RESULTS
Table 4.1. Model Parameters.
Model

Mcca
M ]

Tinit

P (Tinit)

M1a
M2a
M3a
M4a
C1a
H1a
M1b
M1c

1.0
0.1
0.3
3.0
1.0
1.0
1.0
1.0

10
10
10
10
5
30
10
10

1:6  10;2
1:6  10;2
1:6  10;2
1:6  10;2
3:2  10;3
2:1  10;1
7:7  10;2
1:6  10;3

K]

cm g;1 ]2initg  cm;3 ]
1:41  10;19
1:41  10;17
1:57  10;18
1:57  10;20
1:77  10;20
3:82  10;18
1:41  10;19
1:41  10;19

Rccb

Rfc c

AU]

AU]

1:0  104
1:0  103
3:0  103
3:0  104
2:0  104
3:3  103
1:0  104
1:0  104

4.2
4.2
4.2
4.2
3.7
5.1
6.7
3.1

ad d

g  cm;3 ]2  10;13
2  10;13
2  10;13
2  10;13
1  10;13
6  10;13
5  10;14
5  10;13

Le
L ]

f
Lacc
fc
L ]

0.061
0.014
0.031
0.12
0.018
0.24
0.14
0.050

0.04
0.04
0.04
0.04
0.01
0.3
0.1
0.02

Mass of the parent cloud core.
b Radius of the parent cloud core.
c Radius of the rst core.
d Estimated from the numerical results.
e At  = 10;9 g  cm;3 .
cf Evaluated from equation (4.13), with 10 for .
a

Table 4.2. Elapsed Times Assigned to Figs. 4.1-4.3
Time #105yr]

Time/tff a

Line Typeb.02.04830
2.16677
2.171392.173012.181222.18725.01.15643
1.223311.225931.22684
1.23147
1.23488

thick

ff is the free fall time  (3 =(32Ginit)) 2 = 1:77  105 yr

at1

b For correspondence to Figures 4.1 and 4.3.

69

dotteddashedCHAPTER 4. THE FIRST COLLAPSE
reaches unity is given:

Tinit );1: 1 = 4:3  10;13 g  cm;3 ( 1:6  10;2Pcm2  g;1 );2( 10
K

(4.5)

Equations (4.3) and (4.4) correspond to equations (C.13) and (C.15) in Appendix
C. (Temperature dependence of the opacity is also parameterized in Appendix C.) One
shall nd that these critical densities provide excellent analytical predictions to numericalresults in x4.3 below. More detailed discussions will be made in Appendix C.

70

4.3. PROPERTIES OF THE FIRST CORE

4.3 Properties of the First Core
In this section we dene the mass and radius of the rst core from the numerical results.
The radius and mass of the rst core have signicance for determination of the accretion
luminosity, which can relate the theoretical results with observational data. In this chapter
we focus our investigation on estimating the size, mass and accretion luminosity of therst core, and discuss possible interpretations of observational evidence, taking SEDs as
well into consideration.

4.3.1 Denitions of the Core RadiusFirst we argue how to dene the radius of the rst core, Rfc. Three dierent (i.e., dynamical, thermodynamical, and geometrical) denitions are examined as follows.
Def. 1. Rfc is dened where the gas pressure is balanced with the ram pressure of theinfalling envelope.
Def. 2. Rfc is dened at the position of the peak of jds=dtj, where s is the specic entropy(see Figure 4.1h).
Def. 3. Rfc is considered as the position where jd ln =d ln rj has the maximum value (seeFigure 4.1b).
Figure 4.5a shows Rfc in each denition as a function of c . The lines converge afterthe rst core is formed (c > 10;11g  cm;3 ) and indicate a few AU for Rfc . We consider
that it is most desirable to adopt Def. 1 (Rfc  r(Pgas = Pram)), since it shows anintermediate value among the three and is in the best accordance with the analyticalprediction estimated in x4.4.1 below.
In the remainder of this section, the dependences of Rfc on the mass, initial temperature and opacity are argued. The model parameters are tabulated in Table 4.1. In eachmodel the calculation is performed until the central temperature reaches 2000K, because
at this temperature the rst collapse should actually have been succeeded by the second
collapse, which is not the subject in this chapter.

4.3.2 Dependence on Cloud Mass
We present here the radii of the rst core for the parent cloud cores of dierent masses(Mcc ). Models are denoted by M2a, M3a, and M4a for Mcc = 0:1M, 0:3M , and 3M,
respectively, as well as Model M1a for Mcc = 1M . The radius of the rst core in eachmodel is dened according to Def.1 in x4.3.1 and is plotted in Figure 4.5.
71

CHAPTER 4. THE FIRST COLLAPSE

1.2

(a)

1.2

1

1

.8

.8

.6

.6

.4

.4

-12

-11

-10

-9

1.2

-12

(c)

-11

-10

-9

1.2

1

1

.8

.8

.6

(b)

(d)

.6

T_init = 10K
T_init = 5K

.4
-12

T_init = 30K
-11

.4
-10

-9

-12

-11

-10

-9

Fig. 4.5.| Radius of the rst core, Rfc. (a) Various denitions for Rfc are shown. The
denition R(Pgas = Pram) is represented by solid line, R(jd ln =d ln rj = max) by dashedline, and R(jds=dtj = max) by dotted line. R(Pgas = Pram) is adopted as Rfc in thefollowing panels. (b) Rfc for dierent masses of the parent cloud core and for an initiallyinhomogeneous density distribution, (c) for dierent initial temperatures, and (d) fordierent opacities. See text for denotation of72Models a, b, and c. The analytical estimation
for Rfc assuming 10K for Tinit and 10;13 g  cm;3 for ad is shown as well (dot-dash lines).

4.3. PROPERTIES OF THE FIRST CORE
Interestingly Rfc as a function of c is not sensitive to Mcc . That is, Rfc does notdepend on the initial density and the size of the parent cloud core. Consequently, themass of the rst core also does not depend on the mass of its parent cloud core. A
qualitative explanation of the result is briey summarized as follows. A given mass of theparent cloud core determines the initial density and its radius. In the isothermal phase thedensity evolution of the collapsing cloud core \forgets" its initial condition and develops
into centrally peaked density prole (which resembles the Larson-Penston solution) in thesound crossing time over the entire cloud core (as noted also by Bodenheimer & Sweigart
1968). Therefore the evolution, regardless of mass, converges to the sequence which isdetermined uniquely by the initial temperature. The physical processes which change theisothermal phase into the adiabatic phase are aected only by the initial temperature andopacity, according to equations (4.3) and (4.4), and thus the properties of the rst core
are aected by Tinit and P but are independent on the mass of the cloud core.
For comparison, we carried out a calculation with an inhomogeneous initial density
distribution. We provided an isothermal equilibrium conguration with the central densityof 10;17g  cm;3 as the initial density distribution, for which the outer radius truncated at
104AU. The mass of the cloud is 0:87M. The resultant value of Rfc shows no signicantdierences from the other cases, as seen in Figure 4.5b, which supports the \forgetful"nature of collapsing cloud cores. However, the dierence in the initial density structurecauses a substantial change in the luminosity curve after the second collapse, as we will
see later in Chapter 5.

4.3.3 Dependence on Initial Temperature
The Taurus-Auriga molecular cloud has a typical temperature of 10K, which motivatedus to choose the initial temperature of Tinit = 10K as a typical value. However, theexistence of massive stars in the neighborhood of star-forming clouds would raise initial
temperatures for the collapse of cloud cores such as in the Orion clouds. In contrast, theinterior of cloud cores may well have a lower temperature than the typical value of theentire cloud because visual or UV photons from neighboring stars would not reach deep
inside cloud cores. Indeed, the value of 9:6  10;2erg  s;1g;1 for CR , which was foundto be required to sustain the initial cloud core at T = 10K, is larger by two orders ofmagnitude than the typical cosmic ray heating rate.
Hence, model calculations for Tinit = 5K (Model C1a) and 30K (Model H1a) are
carried out to be compared with the typical case of Tinit = 10K. The other conditions
are the same as Model M1a. A more realistic condition will be considered for the initial
temperature in Chapter 5.
73

CHAPTER 4. THE FIRST COLLAPSE

Case for Tinit = 5K
Equation (4.3) claims that =; decreases quite sensitively with decreasing initial temperature. The Planck mean opacity, moreover, has a lower value for a lower initial temperature, which makes =; diminish with Tinit more sensitively. For Tinit =5K, we have=; = 8:8  10;18g  cm;3 
noting that P (5K) is calculated to be 3:2  10;3cm2  g;1. It is found that lowering Tinit
by a factor of two yields a =; which is lower by three orders of magnitude. In contrast, 1 has a negative dependence upon the initial temperature (see eq. #4.5]). For Tinit =5K
one nds
 1 = 2:2  10;11g  cm;3 :In this case  1 exceeds =; by six orders of magnitude, in contrast to the case ofTinit =10K, where they are comparable to each other. This contrast causes a large dierence in thermal history between the two cases.
Figure 4.6 illustrates the thermal evolution for Tinit =5K. One can see that the isothermal phase ceases rather early, when c reaches =;, and that the central temperature
increases gradually during a long transitional phase before the adiabatic stage is enteredwhen c exceeds  1 . In the transitional phase, where =; < c <  1 , the compressional heating rate attempts to dominate the radiative cooling rate, but is pursued
immediately by it in response to the raised central temperature. Consequently ;g and "thare almost balanced with each other as they increase gradually with c before the opticaldepth becomes unity.
The isothermal approximation, which many theoretical works suppose for collapse
calculations, is often considered to be valid when c <  1 10;13g  cm;3 . However,under the situations where =;  1 as in this Model C1a, the transitional phasebreaks the isothermal condition far before the adiabatic stage begins, which implies thatthe isothermal approximation is not so valid even for c <  1 . Especially, this aspect iscritical for cylindrical collapses, where even a slight hardening of the isothermal equationof state can decelerate the gravitational collapse (Inutsuka & Miyama 1997). Further
discussions will be made in Appendix C.

Case for Tinit =30K
If the cloud core is surrounded by the energy source such as OB stars, the temperatureof the cloud can be larger from the beginning of the collapse. Indeed, cloud cores in theOrion clouds should be exposed to hard radiation elds, and the average temperature of74

4.3. PROPERTIES OF THE FIRST CORE

15
cosmic ray heating(a)radiative heating10

radiative coolinggas compression
5

0

-5
-20

-18

-16

-14

-12

-10

-18

-16

-14

-12

-10

(b)

3

2.5

2

1.5

1

-20

Fig. 4.6.| Thermal evolution at the center of a collapsing cloud core for Model C1a(Tinit = 5K). See also the caption of Fig. 4.4.
75

CHAPTER 4. THE FIRST COLLAPSE
the Orion clouds is considerably higher than that of dark clouds. Therefore the calculation
with higher initial temperature should be instructive for the understanding the eect ofdierent environments. Thus in Model H1a, we have removed CR and introduced incident
blackbody radiation of Tinit instead (i.e., Case 2 in x4.1.2) for Tinit = 30K. The results areshown in Figure 4.7.
For Tinit =30K and P (30K) = 2:1  10;1cm2  g;1,dif = 9:2  10;14g  cm;3 :Figure 4.7a clearly shows the turning point between the isothermal and adiabatic phasesat c = dif . In great contrast to the case of Tinit = 5K, where the transitional phaseinterrupts the isothermal stage before  1, the isothermal phase lasts even after the
optical depth exceeds unity.

4.3.4 Dependence on Opacity
We consider here the inuences caused by the ambiguity of the dust opacity for submillimeter and millimeter wavelengths. As described in x4.1.3, the three models of dierent opacities are examined. We provide Models M1b and M1c, that correspond to Model(b) and (c) in x4.1.3.
The opacity of Model (b) gives P = 7:7  10;2cm2  g;1 , which is larger by a factorof 5 than Model (a). Equations (4.3) and (4.5) for this case yield

=; = 3:2  10;13g  cm;3  1 = 1:9  10;14g  cm;3 :Considering  1 < =;, the critical density for violation of isothermality isdif = 4:2  10;14g  cm;3 :For Model (c), where the dust-to-gas ratio is decreased by a factor of ten from Model(a), P equals 1:6  10;3cm2  g;1 , which is of course a tenth of the value of Model (a).
Equations (4.3) and (4.5) give=; = 1:4  10;16g  cm;3  1 = 4:3  10;11g  cm;3 :Figure 4.9 shows that the transitional phase appears again where =; < c <  1 (seex4.3.3). It would be worth noting that the transitional phase also appears for the lower
opacity model in Figure 3 in Boss(1988).
76

4.3. PROPERTIES OF THE FIRST CORE

15
cosmic ray heating(a)radiative heating10

radiative coolinggas compression
5

0

-5

3.5

-16

-14

-12

-10

-16

-14

-12

-10

(b)

3

2.5

2

1.5

Fig. 4.7.| Thermal evolution at the center of a collapsing cloud core for Model H1a(Tinit = 30K). This is a case where CR is absent and the incident radiation eld ispresent. Note the dierence in Panel (a) from the other models. See also the caption ofFig. 4.4.
77

CHAPTER 4. THE FIRST COLLAPSE

15
cosmic ray heating(a)radiative heating10

radiative coolinggas compression
5

0

-5
-18

-16

-14

-12

-10

-16

-14

-12

-10

(b)
3

2.5

2

1.5

1
-18

Fig. 4.8.| Thermal evolution at the center of a collapsing cloud core for Model M1b(larger opacity). See also the caption of Fig. 4.4.
78

4.3. PROPERTIES OF THE FIRST CORE

15
cosmic ray heating(a)radiative heating10

radiative coolinggas compression
5

0

-5
-18

-16

-14

-12

-10

-8

-16

-14

-12

-10

-8

3.5
(b)
3

2.5

2

1.5

1
-18

Fig. 4.9.| Thermal evolution at the center of a collapsing cloud core for Model M1c(smaller opacity). See also the caption of Fig. 4.4.
79

CHAPTER 4. THE FIRST COLLAPSE

4.4 Discussions
In the previous section, we noted the importance of determining the radius and mass ofthe rst core and dened them in the numerical results. We will show in this sectionthat they can be derived from semi-analytical considerations. Possible interpretations of
the observational YSOs are also discussed, on the basis of the accretion luminosities and
SEDs in our results.

4.4.1 Analytical Derivation of the Core Radius
We approximate the rst core as a polytropic sphere in equilibrium. Though the rst
core is not strictly isentropic, one nds that this approximation provides an excellent
prediction of the radius of the rst core. The equilibrium congurations are expressed by
the Lane-Emden equation
1 d ( 2 d ) = ;n
(4.6)
 2 d d
where we follow the standard notation of Chandrasekhar (1967). In the present calculations the ratio of specic heats,  , is xed to 5/3, i.e., n =1.5. The radius of a polytropic
sphere in equilibrium, RLE , is1)K n1 ;1 12
RLE (c ) = # (n 4+ G
(4.7)
c ] 1
where 1 = 3:65375 for n = 1:5 (Chandrasekhar 1967) and K  p= , which relates withthe entropy that the system has. One immediately obtains K as follows:

K = c2s01ad; (4.8)

where cs0 is the sound speed at the initial isothermal phase and ad is dened at theintersection of extrapolated lines of the isotherm and the adiabat in the log-logT plane(as illustrated in Figure 4.10). The log-logT planes of numerical solutions are shown in
Figures 4.4, 4.6, 4.7, 4.8, and 4.9, which read the typical value of 10;13 g  cm;3 for ad astabulated in Table 4.1. Substituting equation (4.8) into equation (4.7), one obtains theanalytical description of Rfc as follows.
1
Tinit ) 21 (ad); 3 :
Rfc (c) = 5:3AU( 10;9gc cm;3 ); 61 ( 10
(4.9)
;
13
;
3
K 10 g  cm
Figure 4.5 delineates Rfc (c) in equation (4.9) (dot-dash line), with Tinit and ad set
to be 10K and 10;13g  cm;3 , respectively. One will nd that the analytical prediction isin excellent accordance with the numerical results.
80

4.4. DISCUSSIONS

log T

ρad
log

Fig. 4.10.| Denition of ad (see text).

81ρCHAPTER 4. THE FIRST COLLAPSE
Apparently ad relates with =;,  1, and/or dif . Because all of =; ,  1 anddif depend solely on Tinit and P as seen in equations (4.3) and (4.4), Rfc also depends
only on Tinit and P when c is given, which supports the numerical results that Rfc isindependent of the mass of the parent cloud core (see x4.3.2).
Interestingly, ad seems never to depart from 10;13g  cm;3 further than an orderof magnitude for any models examined here, in contrast to =; and  1 themselves,
which vary quite sensitively to the initial temperature. Moreover, considering the weak
dependences of Rfc on c , Tinit, and ad in equation (4.9), one will nd that the typical valueof 5AU for Rfc would be a reasonable estimate for any cases unless extreme conditions
for Tinit and P are supposed. Indeed, Figure 4.5 shows that the radius of the rst core
converges to several AU in every model.

4.4.2 Estimate of the Accretion Luminosity
The emergent ux of the radiation is the most important quantity when we compare ourresults with observational data. In this section we express the total luminosity obtained
from the numerical calculations, in terms of analytic expression.
So-called accretion luminosity is expressed as follows:

GMfcM_ Lacc
fc = Rfc(4.10)

where Mfc and M_ are the mass of the rst core and the mass accretion rate, respectively.
For determining Lacc
fc , Mfc and M_ are estimated as follows.
n
Mfc = 4 c Rfc3 (; 1 dd ) = 1 :(4.11)

M_ =  cGs0 :(4.12)We dene Mfc as the mass of an equilibrium polytropic sphere of the radius of Rfc inequation (4.11), where ((1= )(dn=d )) = 1 = 5:5642  10;2 for n = 1:5 (Chandrasekhar1967). In equation (4.12)  is a dimensionless coecient, which is found to be 10 bycomparison with the numerical results for Model M1a. For another example, Shu(1977)
found  = 0:975 in the expansion wave solution. This dierence arises because theexpansion wave solution assumes an initially hydrostatic singular sphere in contrast to
our homogeneous density distribution. Note that the mass accretion rate would become
smaller with time in much later stages. This aspect will be discussed in detail in thesubsequent chapters.
382

4.4. DISCUSSIONS
Substituting equations (4.11) and (4.12) into equation (4.10), we obtain
2 Tinit 5

c
ad
;2
Lacc
fc = 0:07L ( 10 )( 10;9 g  cm;3 ) 3 ( 10K ) 2 ( 10;13 g  cm;3 ) 3 :(4.13)

Table 4.1 tabulates Lacc
fc evaluated by equation (4.13) as well as the numerical result, L. By
acc
comparing Lfc with L, one can see that equation (4.13) oers a reasonable approximation
for the luminosities in most of the models.

4.4.3 Case for = 7=5
Hydrogen molecules in molecular cloud cores of 10K have the ratio of specic heats,
 , equal to 5/3 because rotational and vibrational levels of H2 are not excited at low
temperature. However, the rotational degree of freedom raises the ratio of specic heats
up to 7/5 before the temperature reaches several hundred K (Black & Bodenheimer 1975).
The calculations presented in this chapter suppose  = 5=3 throughout the evolution, forsimplicity.
To study this eect, another calculation was performed under the same conditions as
Model M1a except that  is set to 7/5 instead of 5/3. We found that Rfc does not show
apparent dierences between the two cases and Mfc is smaller by only a factor of twofor  = 7=5. These results are consistent with past works where the rotational degreeof freedom was taken into account (L69 TW). Therefore this simplied treatment for 
remains useful for this kind of studies.
However we will use a more realistic equation of state in the subsequent chapters,
where we consider the dissociation, ionization, and the rotational and vibrational degrees
of freedom for molecules, as well as the non-ideal eects of the equation of state such as
the pressure ionization and the degeneracy of electrons.

4.4.4 Comparison with Observations
Our non-grey radiation hydrodynamical calculations provide information on thefrequency-dependent radiative ux. Recently a number of SEDs of young stellar objects are available. In this section we examine possible suggestions and interpretations of
the observations from our results.
In Figure 4.3 we presented SEDs for the typical model, whose shape were found to
be quite insensitive to the thermal evolution in the center throughout the rst collapse
because the envelope is opaque. In Figure 4.11, the SEDs of the other models are shown.
The shapes of those SEDs also do not reect the rapid increase of the central temperature
except for Model M1c (unusually small opacity model). Therefore, if we rely only on the83CHAPTER 4. THE FIRST COLLAPSE
shape of the SEDs, it is dicult to distinguish cloud cores which have already formed therst (or possibly more evolved) cores from star-less cores or pre-collapse cores. However
the total luminosity should experience a considerable increase even at this early stage ofthe collapse which is hidden in the opaque envelope.
Ward-Thompson et al.(1995) presented SEDs and luminosities for four \star-less"cores, on the basis of the results of a submillimeter continuum survey. Their masses are
1-2M and are comparable with our models. The SEDs (Figure 5 in Ward-Thompson
et al.1994) can be tted by greybody spectra of T 10-20K, which are quite similar
to our results in Figures 4.3 and 4.11. Furthermore, the estimated luminosities of theobserved star-less cores, 10;1L < L < 2L, seem to be rather higher than expected
for cores which have not yet experienced dynamical collapse. Indeed, Figure 4.12 depicts
the evolutions of the luminosities, which reach at most 10;1L even at the end of therst collapse. These results show a relatively luminous ( 1L ) star-forming cloud core
might be recognized to have already formed the second core at the center. Therefore wesuggest that the objects observed by Ward-Thompson et al.(1994) are possibly identied
more naturally as young protostars which are hidden by opaque envelopes, rather thanas pre-protostellar cores as concluded by Ward-Thompson et al.
Boss and Yorke(1995) calculated SEDs at the rst collapse stage, which is \class;I" stage by their denition, by their three-dimensional numerical code. No essential
discrepancy is recognized between their SEDs and ours, although the SED of the poleon view by Boss and Yorke shows a somewhat more excess of the component of higher
temperatures. (Note that the vertical axis in Figures 1 and 2 in Boss and Yorke denotes
L , whereas F by our denition.) Thus the pole-on view of young stellar objects would
possibly reveal the dramatic evolution at the center of the cloud cores.

84

4.4. DISCUSSIONS

5

5

0

0

-5

-5

-10
11

11.5

12

12.5

13

13.5

-10
11

5

5

0

0

-5

-5

-10
11

11.5

12

12.5

13

13.5

-10
11

5

5

0

0

-5

-5

-10
11

11.5

12

12.5

13

13.5

-10
11

11.5

12

12.5

13

13.5

11.5

12

12.5

13

13.5

11.5

12

12.5

13

13.5

Fig. 4.11.| SEDs which would be observed with the beamsize of 103AU for dierent
models. The SED for Model H1a has been omitted because it supposes a surrounding
isotropic radiation eld.
85

CHAPTER 4. THE FIRST COLLAPSE

.15

.1

.05

0
-12

-11.5

-11

-10.5

-10

-9.5

.15

.15

.1

.1

.05

.05

0
2.16

2.17

2.18

2.19

0
6.48

.05

.05

.04

.04

.03

.03

.02

.02

.01

.01

0
.22

.225

.23

.235

.24

.245

0
.65

-9

6.49

.655

-8.5

6.5

.66

6.51

.665

.67

Fig. 4.12.| Evolutionary curves of the luminosity. L.vs.c (top panel) and L.vs.t (others)
for Models M1a, M2a, M3a, and M4a. The luminosity remains at a constant value since
the collapse is initiated (i.e., t = 0) until the luminosity begins to increase.
864.5. SUMMARY

4.5 Summary
In this chapter, we have investigated the properties of the rst collapse and the rst coreby means of radiation hydrodynamical numerical calculations, supplemented by analyticalconsiderations. Our ndings are summarized as follows.
1. The rst core nally has the mass of 0:05M and the radius of 5 AU for
Tinit = 10K and the typical opacity. If the initial temperature and opacity are given,
these values for mass and radius are uniquely determined and are independent bothof the mass of the parent cloud core and of the initial density prole.
2. We have found the analytic expression for the radius of the rst core (eq. #4.9])
and the accretion luminosity (eq. #4.13]) as functions of the central density. They
are weakly dependent on the initial temperature and opacity. Those expressions areuseful for understanding the numerical results.
3. The shape of the spectral energy distributions (SEDs) does not reect the central
thermal evolution of collapsing cloud cores throughout the rst collapse. This result,of course, does not exclude the possibility that the rst cores might be observed on
the pole-on view, where the envelope material would actually be more depleted than
predicted by spherically symmetrical calculations. The luminosity increases up to
0:1L at the end of the rst collapse phase.
4. Dense cloud cores with luminosities higher than 0:1L possibly have young protostars inside the opaque envelopes, even if observations show no evidence for protostellar activities.
In the subsequent chapter, we will discuss the evolution after the rst collapse, inwhich the dissociation of hydrogen molecules causes the second collapse. Main subjects
there are to examine the properties of the main accretion phase and to answer the question
when a protostar emerges before us in the near-infrared.

87CHAPTER 4. THE FIRST COLLAPSE

88

Chapter 5
THE WHOLE EVOLUTION OF
PROTOSTAR FORMATIONIn Chapter 4 computational results were presented under simplied treatments for theequation of state and opacities in order to clarify the physical properties at the rst
collapse stage. We there neglected dissociation of molecular hydrogen, and thus further
evolution initiated by the second collapse was beyond the scope of Chapter 4.
The whole evolution from a star-less molecular cloud core to a protostar is pursued
in this chapter, where all the physical processes described in Chapter 3 are taken into
account. After introducing the initial and boundary conditions in x5.1, we overview thecomputational results in x5.2 and the behavior of the luminosity curve is discussed in x5.3.
The evolution of SEDs is described in x5.4, and summary and conclusions are remarked
in x5.5.

89

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

5.1 Initial and Boundary Conditions
5.1.1 Initial Density Structure

An appropriate initial condition for gravitational collapse of a cloud core is related withmechanism for the initiation of star formation, which involves a controversial problem.
Star formation may be triggered by external factors such as shock compression of gas due
to supernovae or stellar winds in the neighborhood of star formation regions. Otherwise,a star formation may occur spontaneously because of the plasma drift (or the ambipolardiusion) of magnetic elds (e.g., Nakano 1979 Shu, Adams, & Lizano 1987 Shu 1995),
or the dissipation of turbulence which reduces the eective sound velocity in cloud cores(Nakano 1998). The \spontaneous star formation" process would produce nearly hydrostatic cloud cores as the initial condition for the dynamical collapse, while the \triggered
star formation" would not necessarily. Although some latest observations resolved spatial
structures of molecular cloud cores, the initial condition based on observational evidence
is still dicult to determine.
L69 and WN supposed a homogeneous distribution as the initial density structure in
their calculations for protostar formation. A homogeneous cloud core is not in hydrostatic
equilibrium and hence the gravitational collapse is initiated with a nite acceleration.
The dynamical evolution is known to be insensitive to the initial density conguration(Chap. 4), unless the gravitational energy is signicantly larger than the thermal energy
at the initial stage. Nevertheless, dierent initial conditions for the density distribution
may cause a dierence by a factor in the mass accretion rate, which possibly hastens orretards the growth of a protostar. Furthermore, line spectral proles reect the variationsof the initial condition more sensitively.
In this paper, two limiting cases are considered for the initial density distribution: thehomogeneous and the hydrostatic. The homogeneous cloud core has 1M for the total
mass and is bounded by the xed boundary at R = 104AU, which yields the densityof 1:415  10;19g  cm;3 . This condition is equivalent to those by L69 and WN. The
hydrostatic cloud core is obtained by solving the equation of hydrostatic equilibrium for
an isothermal sphere:
1 d ( r2 d ) = ; 4 G 
(5.1)
r2 dr  drc2s
where cs is the isothermal sound speed. We numerically integrate equation (5.1) outward
from the center with the same central density as in the homogeneous model (i.e., c 
1:415  10;19). For the initially hydrostatic models the outer boundary is xed at R =4  104AU and the total mass is then found to be 3:852M. We adopt a larger value ofthe total mass than in the homogeneous model because the mass should exceed a critical90

5.1. INITIAL AND BOUNDARY CONDITIONS
value so that an isothermal hydrostatic core is gravitationally unstable (Bonner 1956
Ebert 1955).

5.1.2 Initial Temperature Structure
For the thermal energy balance at the initial stage, we consider a more realistic condition
than in Chapter 4 and than in any other works concerned with protostar formation. Previous studies assumed isotropic blackbody radiation of 10K surrounding the cloud core
(L69 WN) or heat sources embedded homogeneously with no incident radiation (MMI
Chap. 4). In contrast, the initial temperature prole here is taken as the thermal equilibrium state between the thermal cooling of dust grains and the heating by the externalradiation and cosmic rays (Hayashi & Nakano 1965). The external radiation eld considered here is the sum of the background radiation of 2.7K and ambient stellar radiation.
Stellar photons are given as diluted blackbody radiation of 6000K with the intensity integrated over frequency of 2  10;4erg  s;1 cm;2 for the initially homogeneous model and1  10;4erg  s;1 cm;2 for the initially hydrostatic model. We obtain the equilibrium state
by solving numerically the radiation moment equations and the energy equation underthe xed density distribution. The stellar radiation sustains the cloud core to be nearly
10K with the cosmic ray heating rate of CR = 6:4  10;28n(H2)erg  s;1 cm;3 (Goldsmith& Langer 1978). The initial cloud core has the optical depth close to unity to opticalphotons. Incident optical photons, therefore, are attenuated as passing through the cloudcore and the temperature in deep interior is slightly lower than the temperature near the
cloud edge (see Figure 5.1 below).
In order to consider the thermal balance between the thermal cooling and externalradiation, a frequency-dependent treatment for radiation transport is imperative. The
dust opacity for optical wavelengths is larger by four orders of magnitude than the value
for the sub-millimeter (see Fig. 3.3) and then gray approximations are totally invalid.

91CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

5.2 Overview of Results
We present computational results of protostellar collapse of an initially homogeneous cloudcore. Results are found to be qualitatively the same in the initially hydrostatic model
except for the behavior of the luminosity curve. Hence we describe the evolutionary
properties only for the initially homogeneous model in this section, and the luminosity
curve will be discussed for each case in x5.3.
Figure 5.1 shows the evolution before the main accretion phase. Evolution at theearly stage, i.e., the rst collapse stage, proceeds along the line described in Chapter
4, where the density distribution develops almost isothermally according to the LP-type
solution. A major dierence arises in the thermal evolution because of the more realistic
assumption for the external radiation eld (see x5.1). In contrast to the previous models
where the temperature never descends below the initial value, temperature in the central
region once decreases to 5K and then rises again. The reason for the decrease of T
is that increasing column density of the collapsing cloud core prevents optical photons
from penetrating into the interior. The cloud core then cools down to 5K, at which thecosmic ray heating support the temperature against the thermal cooling. However, thecompressional heating rate (;g ) immediately overtakes the thermal cooling rate (4 P B )and the central temperature begins to increase (Fig. 5.2). This picture is supported by
the result in x4.3.3, where the isothermal evolution of a cloud core of T = 5K was foundto cease very early at c =; = 8:8  10;18g  cm;3 . Cooling of gas due to collisions
with dust grains ("gd in Fig. 5.2a) is closely coupled with the radiative cooling rate andthe compressional heating rate for c > 10;17g  cm;3 .
Although the central temperature begins to rise far before the optical depth arrives
at unity, the evolution becomes completely adiabatic when   1 at c  10;12g  cm;3 ,at which the radiative heating rate ( E cE ) overtakes the radiative cooling rate. For
adiabatic evolutions the eective adiabatic exponent, e  d ln T=d ln  + 1, equals theratio of specic heats (e = 7=5 for T > 100K) as one can see in Figure 5.2b. The
gravitational collapse is decelerated by the increasing gas pressure, which leads to theformation of the rst core. The physical properties of the rst collapse and rst core werediscussed in detail in Chapter 4. Typical scales of radius and mass of the rst core agreewith the values obtained there. We do not repeat further discussions on the rst collapse
here.
The second collapse takes the place (curves 6-8 in Fig. 5.1) when the central density
and the central temperature reach c 10;7g  cm;3 and Tc 2000K, where dissociationof molecular hydrogen reduces the eective adiabatic exponent to 1:1. Dissociationof H2 ceases when the dissociation degree approaches unity in the center at c 10;292

5.2. OVERVIEW OF RESULTS

0
9

9

(a)

(b)
7

4
-5

log T [K]

7

3

5
-10

5

2
-15
3
1

1

3

1

-20
-4

-2

0
log R [AU]

2

4

-4

-2

0
log R [AU]

2

4

0
(c)
20

(d)
-1

9

-2

9

15

-3

6

3

1

10
8
-4
5
7
-5

0

1
-4

-2

0
log R [AU]

-6
2

4

-4

-2

0
log R [AU]2

4

Fig. 5.1.| Numerical results for the initially homogeneous model before the main accretion phase are shown. The succeeding evolution is depicted in Figure 5.4. Distributionsof temperature (a), density (b), infall velocity (c), and mass (d) are presented. Thick line
denotes the initial condition. Elapsed times are tabulated in Table 5.1.
93

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

Table 5.1. Elapsed Times for the Initially Homogeneous Modelt ; t0 #yr]a;1:7522  105;2:9237  104;4:1894  103;6:5192  102;4:3844  102;1:2656  100;3:1250  10;2
;1:5625  10;2
1:5625  10;2
1:8984  101
4:9265  103
2:2958  104
1:3788  105

Label
1
2
3
4
5
6
7
8
9
10
11
12
13

a The oset of time, t  1:7526  105 yr, represents the instance when the second collapse
0

begins.

94

5.2. OVERVIEW OF RESULTS

20

(a)

10

0

5

-15

-10

-5

0

-15

-10

-5

0

(b)

4

3

2

1

Fig. 5.2.| Thermal evolution at the center of the cloud core is depicted for the initiallyhomogeneous model. The upper panel shows the energy exchange rate between gas anddust ("gd), the radiative heating rate( E cE ), the radiative cooling rate(4 P B ), and thecompressional heating rate of gas(;g ) as functions of the central density. The lower panel
shows the central temperature as a function of the central density.
95CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION
and Tc 104K. This is the instance when the second core is formed. The dynamical
time of the second collapse is of the same order as the free-fall time corresponding to = 10;7g  cm;3 , which is 0.1 yr. The second collapse therefore begins and ends almost
instantaneously compared to the rst collapse, which lasts for about t 105 years (eq.
#4.2]). The adiabat recovers its steepness after the second collapse is halted, but e does
not reach 5/3 as expected for monoatomic gas because the partial degeneracy of electrons
begins to dominates the pressure in the center.
All the component of dust grains evaporate at c 3  10;7g  cm;3 and the radiative
heating and cooling rates drop suddenly. The increasing density and temperature, however, again raise the radiative heating and cooling rates. The compressional heating rateof gas increases with c until the collapse is almost halted at c 1g  cm;3 , where thecompressional heating rate drops by more than 6 orders of magnitude and the quasi-static
contraction of the protostar takes the place of the gravitational collapse.
L69 and WN reported rebound of the second core just after its formation. In our
results, indeed, a shock wave propagates outward as if the second core rebounded (see Fig.
5.1c). The rebound, however, does not actually occur because the velocity always remains
negative everywhere. The feature that the infalling gas is decelerated (but not reversed)
by an outgoing shock wave is observed more clearly in Figure 5.3, which illustrates the
trajectories of mass elements during the second collapse. Here the oset of time, t0 
1:7526  105 yr, is taken as the instance when the second collapse begins. As expected,
t0 is very close to tff dened by equation (4.2). A shock wave propagates outward as
settling the infalling material onto the second core.
The very short second collapse phase is followed by the main accretion phase (curves
9-13 in Fig. 5.4). The duration time of the main accretion phase is characterized by
the accretion time, tacc , which is dened by the time required for a mass element in theinfalling envelope to be accreted by the central star:tacc ' vR 
(5.2)

q

where R is the radius of the infalling envelope and v  2GM=R is the free-fall velocity.
Rearranging equation (5.2), we obtain
s

R 
(5.3)
tacc ' 2GM
which gives essentially the same value as the free-fall time (eq. #4.2]) except for numericalfactors of order of unity.
Figure 5.4 shows the evolution succeeding to Figure 5.1. The second core, which has
963

5.2. OVERVIEW OF RESULTS

2

1

R [AU]

0

-1

-2

-3

-4
-3

-2.5

-2

-1.5

-1

-.5

0

Fig. 5.3.| Trajectories of mass elements of the cloud core are delineated just after thesecond collapse. The innermost and uppermost mass shells correspond to Mr = 10;3M
and Mr = 2  10;2M, respectively. Dierence in mass between the neighboring shells isconstant in logarithmic scale, where % log Mr = 5:42  10;2. See legend of Table 5.1 forthe denition of t0.
97.5

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

6

12

13

13

0

(a)

(b)

9
10

5

-5

9
4
log T [K]

9
10

-10
3

2

9

12 13
10

13

-15

1
-20
-4

-2

0
log R [AU]

2

4

-4

-2

0
log R [AU]

2

4

2

4

0
13
250

(c)

(d)
-1

11

200

10
-2
9

150
-3
100

13
-4
11

50
-5

10
9
0

-6
-4

-2

0
log R [AU]

2

4

-4

-2

0
log R [AU]Fig. 5.4.| Numerical results for the initially homogeneous model in the main accretion
phase are shown. The preceding evolution is depicted in Figure 5.1. Elapsed times aretabulated in Table 5.1.
98

5.2. OVERVIEW OF RESULTS
now become a protostar, grows monotonically in mass as accreting the infalling material(Fig. 5.4d). The radius of the second core is settled at 2  10;2AU or 4R during themain accretion phase. This value of the protostellar radius is comparable with WN and
SST.
The internal structure of the protostar shows an o-centered peak in the temperature distribution, as was also found by WN and SST. The o-centered maximum of
temperature arises because of partially degenerate electrons in the central region. Since
the pressure of degenerate gas depends only weakly upon temperature, the gravitational
contraction of the protostar does not raise temperature eectively where degenerate gas
dominates the pressure. In the o-centered region where the density is not high enoughand the thermal energy dominates the Fermi energy, the gas pressure is more sensitive totemperature and the o-centered temperature becomes higher than in the center.
The density structure in the infalling envelope is approximately  / r;1:5 in the
inner region and  / r;2 in the outer, and the boundary between the both regions shifts
outward with time. Although this feature resembles the Expansion-Wave collapse solution
by Shu(1977), it represented the continuous counterpart of the LP solution obtained by
Hunter(1977) rather than the Expansion-Wave collapse solution or another member ofa family of the solutions by Shu. Shu's criticism of the unphysical aspect of the LP
solution that the infall velocity is nite ( 3:3cs) at the innite radius is inappropriate
here because the velocity is xed at zero at a nite radius in our calculations.
Plateau in the mass distribution (Fig. 5.4d) indicates the current mass of the central protostar. The protostellar mass monotonically reaches 0:73M at the end of thecalculation.
Figure 5.6 depicts spatial distributions of the Eddington factor in the main accretion
phase. The Eddington factor is kept at 1/3 in the protostellar interior (R < 10;2AU)
where the optical depth is quite large ( 1015). At the accretion shock front on the protostellar surface the Eddington factor sharply drops below 1/3, which occurs because theoptically thin radiative relaxation layer is lled with radiation hotter than in the preshockgas and the post-relaxation layer (Fig. 5.5, see also Balluch 1988). The Eddington factor
rapidly rises toward 1 in the opacity gap between the protostellar surface and the dust
sublimation front, but it decreases again onto 1/3 near the dust sublimation front because
the opaque dusty envelope reprocesses the outwardly-peaking radiation from the center
into a nearly isotropic eld of the thermal radiation emitted by dust grains. The Eddington factor increases with radius once again as the dusty envelope becomes diuse andthe outwardly-peaking radiation dominates the isotropic eld. At the nal stage (curve13) the Eddington factor departs from 1/3 even near the dust sublimation front. This is99

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION
f=1

µ
r

1/3 < f < 1µ
r

f = 1/3µ
r

0 < f < 1/3

µ
r

f=0

µ
r

Fig. 5.5.| A schematic picture for the Eddington factor is illustrated. Arrows represent
the angular distribution of the specic intensity. A special case where the Eddington
factor drops below 1/3 as well is shown.
100

5.2. OVERVIEW OF RESULTS
because the dusty envelope becomes optically thin for its own thermal emission due tothe depletion of the infalling material.
Figure 5.6 also shows the distribution of the mass accretion rate, M_ . On the protostellar surface M_ (r) has a sharp spike where the infalled material through the diuse
envelope settles onto the dense protostar. The mass accretion rate in the infalling envelope (R 10;2AU) is 10;5M yr;1 independently of radius except near the outeredge in the early main accretion phase (curves 10-12). This value of M_ is consistent with
SST, who adopted M_  10;5Myr;1 c3s0=G as a free parameter, where cs0 is the soundvelocity in the initial cloud core, However, M_ decreases to 10;6M yr;1 (curve 13)
at later stages in the main accretion phase. Our calculation supposes no incident mass
ux through the outer boundary and then the mass accretion rate should decline after asignicant fraction of the total mass is accreted by the protostar.

101

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

1

Eddington factor

.8

13

.6
11 12

.4
10

-4

-2

0
log R [AU]

2

4

0

2

4

-4

10
11

-5

12

13
-6

-4

-2

Fig. 5.6.| The upper panel depicts the Eddington factor as a function of radius in themain accretion phase. Elapsed times are tabulated in Table 5.1. The lower panel shows
the mass accretion rate, i.e., M_ (r)  4 r2(r)(;v(r)).

102

5.3. LUMINOSITY CURVES

5.3 Luminosity Curves
While the luminosity remains < 0:01L during the rst collapse as was shown in Chapter
4, it rapidly rises when the second core, or a protostar, is born. The luminosity reaches
the maximum value of  25L and then decreases gradually in the initially homogeneousmodel (Fig. 5.7). The luminosity of a protostar surrounded by an infalling envelope isdominated by the accretion luminosity:_ GM_ Z t
GM
M

_ 0
Mdt
(5.4)
Lacc = R = R t0
where M is the protostellar mass, R is the protostellar radius, and M_ is the massaccretion rate onto the protostar. As the protostellar radius is settled at 4R during themain accretion phase (x5.2), the shape of the luminosity curve reects the time evolutionof M and M_ , which closely relate with each other. The initial rise in the luminosity
curve is produced by the rapid growth of the protostellar mass. The luminosity declines
after then because of the decreasing mass accretion rate, which is seen as a declining
solid curve in the lower panel in Figure 5.7. The average mass accretion rate reads
5:7  10;6M yr;1 from the gradient of the M -t curve, while M_ just after the birth of a
protostar is 1:3  10;5Myr;1.
A homogeneous cloud core begins the gravitational collapse with inward acceleration
because it has no pressure gradient and the force balance is not satised initially against
the self-gravity. On the contrary, collapse of a hydrostatic cloud core is initiated without
any acceleration, which yields a lower mass accretion rate than in the initially homogeneous case. Although the dynamical evolution \forgets" the initial density prole inthe sound-crossing time (Chap. 4), small surviving inuence on the mass accretion ratecauses a qualitative dierence in the behavior of the luminosity curve as seen in Figure
5.7.
In contrast to the initially homogeneous model, the initially hydrostatic model produces neither a rapid rise at t ; t0 = 0 nor a decrease after the maximum in the luminosity
curve, and instead shows a gradual and monotonic increase of the luminosity. The M -tdiagram presents a linear growth of the protostellar mass for this model, which indicatesthat the mass accretion rate remains nearly constant. Lack of the rapid growth of M
disables the abrupt rise in the luminosity curve at an early stage, and prevents immediate
exhaustion of the infalling material which is responsible for the decrease of the luminosity.
The average mass accretion rate is 5:3  10;6M yr;1, which is very close to the value inthe initially homogeneous model.
SST found that L increases monotonically with time. Evidently this increase of L
shown by SST originates their \steady state approximation" in which M_ is xed at
103

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

1.5

1

.5

0

-.5

-1
0

50000

0

50000

1

.8

.6

.4

.2

0

Fig. 5.7.| The luminosity curve (upper) and the protostellar mass (lower) are shown
as functions of time after a protostar is formed. t0 is dened in Table 5.1. Solid anddotted lines denote the initially homogeneous model and the initially hydrostatic model,
respectively.
1045.3. LUMINOSITY CURVES
10;5Myr;1 throughout the calculation, where Lacc inevitably rises as the central protostar grows in mass. The steady state approximation by SST is valid for our initiallyhydrostatic model but is invalid for the initially homogeneous model. WN adopted ahomogeneous cloud core as the initial condition, and then found the same tendency as in
our initially homogeneous model.
Even for the initially homogeneous model where the luminosity declines with time,the nal value of the luminosity ( 10L) is still too high for the typical value for class
I sources ( 1L, Kenyon et al. 1990). This \luminosity problem" for class I protostars,
which was pointed out by Kenyon et al.(1990), indeed arose in our spherically symmetriccalculation. The luminosity problem, however, can be resolved by considering a possibledierence in the mass accretion rates between the circumstellar disk and the infallingenvelope (Kenyon, Calvet, & Hartmann 1993), which is not taken into account in ourspherically symmetric calculations.
Exhaustion of the infalling material in our calculations is simply due to the condition
that no mass ux should enter through the boundary. Physical mechanism in which themass accretion onto a protostar is halted is more complicated for actual YSOs. Nakano,
Hasegawa, & Norman(1995) showed that the mass outows from YSOs nally halt further
mass accretion by blowing o the residual envelope matter until the envelope becomes nolonger gravitationally bound. This eect would suppress M_ at late stages of the evolution,
which may also contribute to resolve the luminosity problem.
Which of the decreasing or increasing luminosity accounts for observations? In order
to discuss observational counterparts of our results, we present results of SEDs in the next
section.

105

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

5.4 Spectral Energy Distributions
The computational result for the emergent SED of the initially homogeneous model isshown in Figure 5.8. The SEDs are calculated from the monochromatic outgoing radiation
through the cloud edge, convolved by the Gaussian beam pattern with a given FWHM.
In Figure 5.8 an imaginary telescope is assumed to observe the object at a distance of 150
pc with the resolution of 1000 (FWHM). The interstellar extinction is not considered.
At the initial stage the SED reects optically thin thermal emissions by dust grains ofnearly 10K (the curve denoted by '1'). The SED once becomes colder (curves 2 and 3),
due to a slight decrease of temperature mentioned in x5.2, and then gradually becomes
hotter. However, the rst collapse stage (curves 1-6) shows no apparent changes in SEDs
because the optically thick envelope reprocesses hot photons from the central region into
cold emissions in the outer envelope as was mentioned in Chapter 4.
A deep absorption feature at 13.5 Hz (or 10m) corresponds to a striking feature insilicate (see Fig. 3.3). The absorption is extremely deep just before near infrared emission
appears in the SED. The deepness reects the contrast in temperature between the cold
region in the outer envelope, which is traced by the 10m feature, and hotter regions
close to the center where near infrared photons are emitted. The absorption feature
becomes shallower as the material in the envelope depletes because of the accretion onto
the protostar.
The object is detectable in near infrared wavelengths at the last stages (curves 12 and13). The SED labeled by '12' shows typical features for class I protostars, i.e., a large
excess of the far infrared component and undetectability in optical light. The SED labeledby '13' corresponds to a more evolved YSO like a class II object, for which the enhancementof luminosity in the far infrared has disappeared and the peak in F has been shifted to
the near infrared. The SED of '13' shows a at spectrum in the infrared and hence thisobject might correspond to a at spectrum T Tauri star if the optical component appears
in the SED due to bipolar cavities (Chap. 7). We stopped the calculation at this stage
because the further evolution of SED involves a non-spherically-symmetric treatment,
where emissions from the circumstellar accretion disk and photons that escape throughthe bipolar cavities aect the SED in the near infrared and the optical.
The bolometric temperature, Tbol (Myers & Ladd 1993), dened by equation (1.3)
calculated from the obtained SED is plotted in Figure 5.9. The bolometric temperaturemonotonically increases with time. While Tbol remains at 10K before a protostar is formed(t ; t0 < 0), it abruptly rises up to 40K after t ; t0 exceeds 0 and then gradually increases
to 410K. Chen et al. (1995) found that observed class 0 sources show Tbol less than 70K,
which corresponds to t ; t0  2  104 yr in our results. This is smaller by an order of106

5.4. SPECTRAL ENERGY DISTRIBUTIONS

4

2
13
10 11

12

0
5
3
1

-2

4

2

-4
11

11.5

12

12.5

13

13.5

14

14.5

15

11.5

12

12.5

13

13.5

14

14.5

15

-6

-8

-10

-12

-14
11

Fig. 5.8.| Evolution of the spectral energy distribution for the initially homogeneousmodel is shown. Elapsed time of each curve is listed in Table 5.1. The upper paneldelineates the emergent ux, F , in the unit of #Jy], while the lower depicts F in the
cgs unit. An observer is assumed to be at a distance of 150pc with a telescope of theresolution of 1000.
107

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION

100

10

1
-50000050000-50000

0

50000

100

10

1

Fig. 5.9.| Evolution of the bolometric temperature, Tbol , is shown. Solid line delineates
Tbol estimated by the original denition for  given by equation (1.2). The other lines
denote Tbol estimated using nu
 with a upper cuto frequency, fc, in the integral (see text),
where fc is 3  1012Hz for dotted line, 5  1012Hz for short-dashed line, and 1  1013Hz forlong-dashed line. The upper panel is for the initially homogeneous model, and the lower
for the initially hydrostatic model.
1085.4. SPECTRAL ENERGY DISTRIBUTIONS
magnitude than the typical age of class I sources. We call these very young protostars
as \genuine" class 0 objects. We nd no qualitative and quantitative dierences in theevolution of Tbol between the initially homogeneous model and the initially hydrostatic
model, in contrast to the luminosity curve (x5.3).
Computational results for the luminosity and Tbol implies that genuine class 0 sources
should be more luminous than class I sources for the initially homogeneous model, andless luminous for the initially hydrostatic model. Observations have revealed that the
typical luminosities are 1L for class I sources (Kenyon et al. 1990) and 5:7L for class
0 sources (Gregersen et al. 1997). Obviously protostars in a younger phase tend to be
more luminous, which supports the initially homogeneous model. Note that the \homogeneity" considered here is just for convenience and any other initial density structures
which are not in hydrostatic equilibrium would trace a similar dynamical evolution to theinitially homogeneous model rather than to the initially hydrostatic model, as long as theinitial cloud core is (marginally) Jeans unstable. A cloud core that is not in hydrostatic
equilibrium obtains a relatively large M_ because of the initial inward acceleration, whichleads to large luminosity just after the birth of a protostar (x5.3). In contrast, the initiallyhydrostatic model, where a lower M_ admits neither the rapid growth of a protostar nor
the immediate exhaustion of the infalling material, fails to reproduce the qualitative tendency in the luminosity of observed YSOs and instead the luminosity increases graduallyand monotonically. Consequently, star formation generally tends to be triggered under
initial conditions that are not in hydrostatic equilibrium.
Our calculations suppose spherical symmetry. Actually, the same YSO may be observed as either class I or class 0 depending on its inclination angle to an observer because
an edge-on view suers from such a high column density that the object would be seen as
if it were less evolved. Considering this point, we also show in Figure 5.9 the bolometrictemperatures which would be estimated if emissions were decient in high frequencies.
We also calculate the mean frequency, , imposing an upper cuto frequency, fc, that is,
R fc

  R0 f F d 
(5.5)
0 F dinstead of equation (1.2). Figure 5.9 indicates that Tbol is sensitive to fc for 3  1012 <fc < 1  1013Hz and Tbol can be signicantly colder even after t 105 when near infraredcomponents are hidden by locally enhanced column density as in the disk plane. A
fraction of edge-on class I, or \class 0-like class I", would not be large relative to the total
number of class I sources, but it can be substantially large for genuine class 0 sources,
which are less abundant by an order of magnitude than class I sources. Although our
results indicates that protostars younger than 104yr would be observed always as classc109

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATION
0 sources, contamination of older (and thus more abundant) protostars is not negligible
for YSOs with cold Tbol .
A clue for discriminating genuine class 0 objects from edge-on class I objects lies in
the luminosity of individual sources since our initially homogeneous model shows that theluminosity reaches the maximum value of 25L at 2  104yr, which is coincident with thetypical age of genuine class 0, and the luminosity declines as the evolution proceeds. YSOs
with cold Tbol which have luminosities considerably larger than 1L would be genuine
class 0 objects. On the other hand, cold YSOs with comparable luminosities with class I
sources are probably more evolved protostars than genuine class 0 objects. The luminosity
problem would not arise for the protostars substantially younger than tacc 105yr, where
non-spherically-symmetric eects are less important because a signicant fraction of theangular momentum of the rotating cloud core still remains far away from the central
protostar.
Unfortunately, the bolometric luminosity estimated from observations suers froma large uncertainty for such young protostars. A major fraction of the luminosity isemitted in the submillimeter and only an upper-limit value is obtained for the infrared
ux. Furthermore, the luminosity of a YSO on an edge-on view would be underestimatedbecause only a small fraction of the total emission leaks out through the disk plane. Directcomparison of the observed luminosities with our numerical results require a great care.

1105.5. SUMMARY

5.5 Summary
Our ndings in this chapter are summarized as follow.
1. An RHD numerical calculation is carried out to reveal physical processes in theformation of a 1M protostar. The whole evolution are pursued from the beginning of the rst collapse to the end of the main accretion phase. All the physicalphenomena that appear in protostellar formation, (e.g., the dynamical collapse of acloud core, the formation of the pressure-supported core, the accretion shock front
on the protostellar surface, and the supersonically infalling gas) are incorporated ina self-consistent manner in the calculation.
2. For the initially homogeneous model, the accretion luminosity rapidly rises to the
maximum value of 25L just after the formation of a protostar, and declines gradually as the mass accretion rate decreases. In contrast, the luminosity increases
monotonically with time for the initially hydrostatic model. This discrepancy arises
because the mass accretion rate varies depending inward acceleration is zero or
non-zero in the initial condition, which aects the luminosity curve.
3. We conrm that the SED evolves from a 10K greybody spectrum to hotter spectra typical for class I and II sources. The bolometric temperature derived by theobtained SED shows that cold sources with Tbol less than 70K are younger than2  104yr. This indicates that class 0 objects have younger ages by an order ofmagnitude than class I objects. However, Tbol can be signicantly colder even aftert 105 when near infrared components are hidden by locally enhanced column density as in the disk plane. Observed class 0 sources should be the compound of the\genuine" class 0 that is as young as 104 yr and more evolved protostars on edge-on
view. The contamination of older protostars are not negligible because they are
intrinsically abundant than genuine class 0 objects.
4. Since observations indicate that class 0 sources is typically more luminous than class
I sources, we exclude the initially hydrostatic model where the luminosity increases
monotonically with time. The initially homogeneous model, in contrast, is found toshow the tendencies consistent with observations. Star formation therefore tends tobe triggered under initial conditions that are not in hydrostatic equilibrium. The
decreasing luminosity with time implies that signicantly luminous and cold YSOs
( 1L and Tbol < 70K) should be genuine class 0 objects. Less luminous YSOs
are more evolved protostars on edge-on view even if they have cold Tbol.
111

CHAPTER 5. THE WHOLE EVOLUTION OF PROTOSTAR FORMATIONIn the present version of the numerical code, the convective energy transport and thenuclear energy production are not taken into account. The deuterium burning, which is
ignited at T 106K, should occur in the protostar at the nal stage of our calculation.
These eects might change the structure of the protostar but a possible change would
not be large because the available energy from the deuterium burning is not larger than
the released energy by the quasi-static contraction of a protostar. Indeed, SST consideredboth the convection and deuterium burning and they found the protostellar radius of
3 4R, which accords with our result quite well. Furthermore, a slight dierence in the
protostellar structure would give no inuence of the structure of the infalling envelope.
Hence the SEDs and molecular line proles, which are introduced in Chapter 6, would
remain unchanged by incorporating the convection and the nuclear reactions.
We did not encounter the dynamical instabilities that found by Tscharnuter(1987) and
Balluch(1988) (x1.1.1). Therefore, the updated version of input physics is not responsible
for the instability. It is most likely that the large oscillation of protostars is an artifact
by a numerical technique, possibly an adaptive grid scheme.

112Chapter 6
MOLECULAR LINE SPECTRAAs mentioned in x1.2.3, molecular line spectra of protostellar objects often show a doublepeaked prole with a stronger blue peak, which is considered as evidence of infall motions.
We here carry out model calculations for line spectral synthesis on the basis of the resultsof radiation hydrodynamic calculations introduced in Chapter 5.

113

CHAPTER 6. MOLECULAR LINE SPECTRA

6.1 Numerical Scheme
6.1.1 Beyond LVG and Microturbulence Approximations
The rst trial for modeling the double-peaked features in line spectra of collapsing clouds
was done by Snell & Loren(1977). They solved radiative transfer problems supposing
power-law distributions for velocity, temperature, and density of a cloud under the LVG(Large Velocity Gradient) approximation. The LVG or Sobolev approximation reduces
the complexity in radiative transfer problems for moving media to a tractable level byassuming that the intrinsic line width is negligibly small compared to the local Doppler
displacement due to the systematic velocity gradient (Sobolev 1960). Photons are absorbed just at the place where they emitted or escape freely under this assumption, andthen the problem is localized. Following Sobolev's method, Castor(1970) derived moregeneral formulations for the photon escape probability,  , and he found that it leads toa simple function of optical depth  , that is,  = #1 ; exp(; )]= , in case of v(r) / r.
Since Goldreich & Kwan(1974) adopted this formulation for  , the velocity distribution
proportional to radius has been a favorite model in the LVG approximation.
Leung & Brown(1977) pointed out that application of the LVG approximation to line
prole synthesis, as done by Snell & Loren(1977), can lead to serious errors when the
intrinsic line width is no longer smaller than the systematic velocity of the cloud. In thelimit opposite to the LVG regime, one can consider the situation where the velocity eld
is dominated by random velocity such as thermal and/or turbulent motions rather thansystematic velocity. Leung & Liszt(1976) carried out radiative transfer calculations in
the presence of \microturbulence", where the correlation length of the velocity eld isassumed to be small compared to the photon mean free path and hence velocity elds aretreat as thermal Doppler shifts with the eective sound velocity.
For protostellar collapse systematic velocity can be either higher or lower than thethermal (or turbulent) velocity width and neither the LVG nor the microturbulence regime
is valid to calculate line spectra. Therefore we have developed a numerical code for multilevel line transfer calculations which is applicable to any physical structures of density,
temperature, and systematic velocity with an arbitrary line prole function. Our scheme
thus demands neither the LVG nor microturbulence approximation. Non-LTE eects
become important near the cloud edge and hence the level population is solved selfconsistently with the radiation eld. Basic equations and computational procedure aredescribed below.
1146.1. NUMERICAL SCHEME

6.1.2 Basic Equations
The radiative transfer equation is solved under the given absorption coecient,  , andthe source function, S ,2
 @I@r + 1 ;r  @I@ = ; (I  ; S )
(6.1)
where
 = h' (  )(nJ BJJ +1 ; nJ +1BJ +1:J )
(6.2)
4 J +1J
S = (n B nJ +1;AnJ +1JB ) :(6.3)
J JJ +1
J +1 J +1:J
Here  is direction cosine (see Fig. 2.1 in x2.2). The Einstein coecients for the dipole
radiation emitted by rotational transitions between the quantum numbers of J and J + 1
are
4 2
BJ +1J = 323h2c d 2JJ ++13 
3
BJJ +1 = 22JJ +B + 1 J +1J
and3
3
AJ +1J = 16chB
2 (J + 1) BJ +1J 
where d is the electric dipole moment and B is the rotation constant. Only the transitions
between %J = 1 are permitted for the dipole radiation. The rotation constant is denedbyB  8 h2I 
where I is the moment of inertia
of a molecule. The energy level corresponds to the J -thqangular momentum LJ = J (J + 1)h is given by

EJ = BJ (J + 1)hHere we consider only linear molecules.
In order to evaluate  and S , populations at every energy level are solved.
@nJ = n #A
J +1 J +1J + JJ +1J BJ +1J ] + nJ ;1 JJJ ;1BJ ;1J
@t;nJ #AJJ ;1 + JJJ ;1BJJ ;1] ; nJ JJ +1J BJJ +1
X
Xc ]
+NH2 # Jc J nJ ; nJ JJ
where

J

0

0

J

0

0

0

1
1
JJ +1J = 21 d dI  'J +1J (  )
;1
0
Z

Z

115

(6.4)

(6.5)

CHAPTER 6. MOLECULAR LINE SPECTRAand Jc J is the collisional transition rate of J 0 ! J . In our calculations the LHS of theabove equation is set to be zero (i.e., @=@t  0.)
The line prole function at rest frame, '0J +1J ( ), is assumed here to be the thermal
Doppler broadening.
0'0J +1J ( ) = % 1p exp#;(  ;%J +1J )2]
D

D

where

(6.6)

s

%D = J +1J 2kTcm
The prole function is shifted in the frequency space due to bulk motion with the radial
velocity vr .
J'J +1J (  ) = '0J +1J ( ; J +1(6.7)
c vr )6.1.3 Computational Procedure
The computational procedure follows an ordinary "-iteration method. First, the transfer
equation (6.1) is integrated under trial values for nJ in equations (6.2) and (6.3). The
Boltzmann distribution is used as a trial level population. Second, equation (6.4) is solved
to by using JJ +1J which is estimated by the solution of the transfer equation. Newly
obtained nJ provides improved values of the absorption coecient and source function,with which equation (6.1) is solved again. The above procedure is iterated until the selfconsistent solution between the source function and radiation eld is obtained everywherein the collapsing cloud. The criterion for convergence is that the maximum of relative
errors in  and S is below 1 %. Convergence are found to be good enough although anacceleration technique for "-iteration is not introduced (e.g, Rybicki & Hummer 1991).
Together with the radiative transfer equations for N transitions from J =0 to N ,equation (6.4) is solved for each J between 0 and N ; 1 along with the normalization
equation:N
X
nJ = cXnH2 
(6.8)
J =0

where X is the fractional abundance of the molecular tracer relative to molecular hydrogen. A dimensionless factor c is unity when most molecules occupy lower rotational
levels than J = N . Otherwise, c should be adjusted appropriately to give the correct
normalization. Correction by c becomes important when T is large enough to pump up
the populations to higher levels than J = N . In our cases such hot regions are conned
near the central protostar, where the density is so high that the Boltzmann distribution is116

6.1. NUMERICAL SCHEME
achieved. We therefore determine c assuming that nJ obeys the Boltzmann distribution:

c =

1
X
(2J + 1) exp(; EJ )= (2J + 1) exp(; EJ )kT J =0
kT
J =0
N
X

(6.9)

Emergent spectra obtained from the converged solutions are presented in terms of theantenna temperature, TA, versus the Doppler velocity. First, we dene the brightness
temperature as
2
TB = 2c2 k I 
(6.10)
where I is the calculated spectral prole. The antenna temperature is evaluated fromTB :RP ())d) 
TA = sourceR TBP())(6.11)4  ())d)
where P represents a Gaussian beam pattern with a given FWHM, which is given as
the wavelength divided by the diameter of an imaginary telescope. Our numerical code
permits an o-centered beam convolution in addition to the centered beam. The aperture eciency and the beam eciency of the telescope are assumed to be unity in our
calculations.
The Doppler velocity is simply dened by
v = f ; f0 
(6.12)cf0
where c is the speed of light, f is the frequency, and f0 is the frequency at the line center
measured in the rest frame.

117

CHAPTER 6. MOLECULAR LINE SPECTRA

6.2 Calculations and Results
6.2.1 Molecular Tracers

We suppose CS, HCO+ , and H13 CO+ as molecular tracers for collapsing clouds. Rotational transitions of these molecules have the high critical densities, > 105-106cm;3 in nH2 ,at which the collisional de-excitation rate is balanced with the Einstein's A coecient.
They therefore trace dense regions in molecular clouds, which are candidates for star
forming cloud cores. Typical abundances are 10;9-10;8 relative to molecular hydrogenfor CS and HCO+ and smaller by nearly a hundredth for H13 CO+ . For probable valuesof column density of star forming cores, HCO+ lines are typically opaque in contrast to
H13 CO+ lines, which are optically thin because of their low abundance. Molecules are
assumed to be destructed for T > 2000K. Spectral proles do not suer from arbitrariness of the destruction temperature unless it is so low that the region where molecules
are destroyed occupies a substantial volume in the cloud core.
Tables 6.1 and 6.2 show the parameters used in calculations for the molecular tracers.
The collisional de-excitation rates, Jc J (J 0 > J ), are taken from tabulated data by Green& Chapman(1978) for CS, and data from Monteiro(1985) (for J = 0 to J = 3) andGreen(1975) (for J = 4 and J = 5, scaled appropriately) for HCO+ and H13 CO+. The
c (J 0 > J ), is calculated using  c from detailed balance:collisional excitation rate, JJ
J J
0

0

0

c
JJ
2J 0 + 1 exp(; EJ J ):=Jc J 2J + 1
kT
0

0

0

(6.13)

Equation (6.13) assures that the level population obeys the Boltzmann distribution when
the radiative transitions are negligible compared to the collisional transitions.
Since data of the collisional transition rates are given only for low temperatures (e.g.,
30K in Monteirop1985), we extrapolate the de-excitation rates for higher temperatures
proportionally to T assuming that the cross section would not vary drastically withtemperature in the transition rate, h vi. Although this assumption is rather crude, itwould cause no misestimation of the level population because the density is high enough forthe population to achieve the Boltzmann distribution wherever the temperature exceeds
30K in our present problems.

6.2.2 Results
Figure 6.1 shows molecular line proles of HCO+ and H13 CO+ (J =3-2 and 4-3) for the
initially homogeneous model (Chap. 5). Four stages of the evolution were chosen: at theinitial (labeled by '1'), just after the rst core is formed ('6'), early ('12') and late ('13')118

6.2. CALCULATIONS AND RESULTS
Table 6.1. Molecular Parameters 1.
MoleculeCS
HCO+
H13 CO+Ba
24.4956
44.592943.3772

d b
1.958
3.3
3.3

Xc
4  10;91  10;8
1  10;10

Nd
10
6
6

a Rotation constant in GHz.
b Electric dipole moment in Debye.
c Fractional abundance relative to H .
2d Levels of J = 0-N are considered in calculations.

Table 6.2. Molecular Parameters 2.
Molecule

Transition

Frequency#GHz]

Resolutiona#00]

Resolutionb
#AU]

CS

J =1-0J =2-1J =3-2
J =4-3J =1-0J =2-1J =3-2
J =4-3J =1-0J =2-1J =3-2
J =4-348.99115
97.98230
146.97345
195.96460
89.1858
178.3716
267.5574
356.7432
86.7544173.5088
260.2632
347.0176

60
60
60
60
69
35
23
17
71
36
24
19

9000
9000
9000
9000
10350
5250
3450
2550
10650
5400
3600
2850

HCO+
H13 CO+a A 10m telescope is supposed except for CS lines where a constant resolution is given forcomparison with Zhou(1992).
b The object is supposed at the distance of 150pc from an observer , i.e., 100 = 150AU.
119

CHAPTER 6. MOLECULAR LINE SPECTRA

14
12
12
1
10
6
8

12

10

8

13
6

6

4

4

2

2

0

0
-2

0

2

-2

0

2

-2

0

2

1.4
1

1.2

1

.8

.8
.6
.6
.4
.4
.2
.2

0

0
-2

0

2

Fig. 6.1.| Computational results of spectral line synthesis for HCO+ and H13 CO+ (J =32 and 4-3) are shown. Adopted temperature, density, and velocity distributions are basedon the initially homogeneous model described in Chapter 5. Line proles are given for the
initial condition(labeled by '1'), just after the rst core is formed('6'), early ('12') and late('13') stages of the main accretion phase. Refer to Table 5.1 in x5.2 for label numbers.
120

6.2. CALCULATIONS AND RESULTS
stages of the main accretion phase. On the basis of the discussions in x5.4, curve '12'
corresponds to the class 0 or early class I phase and curve '13' to more evolved stages.
Refer to Table 5.1 in x5.2 for label numbers.
While the kinetic temperature is almost constant throughout the cloud core, the excitation temperature, Tex , decreases with radius to 2.7K because the density near the outerboundary is too low to collisionally excite levels against the de-excitations by spontaneous emission. A static cloud core with outwardly decreasing Tex produces a symmetricdouble-peaked prole for an optically thick line, which appears in HCO+ J =3-2 at theinitial stage. For the other lines shown in Figure 6.1 the initial proles are resolved out.
As the evolution proceeds, double-peaked proles in a blue asymmetry appear in
optically thick lines (HCO+) as also found by past studies. The most intensive feature isseen in an early stage of the main accretion phase, when the maximum value is achievedin the luminosity (Fig 5.7 in Chap. 5). At the end of the main accretion phase, the peakintensity decreases again because of the depletion of material in the infalling envelope. On
the other hand, line widths continue growing throughout the evolution. Wings extend tov = 2km/s in the main accretion phase. Wing components survive even after the peakintensity is diminished by the depletion of the infalling gas, because the infall velocity
becomes larger as the gravitational potential is deepened by growth of the protostellar
mass.
The optically thin lines (H13 CO+ ) show single-peaked proles with a slight blue asymmetry. The line width increases monotonically with time similarly to the optically thicklines (Table 6.3 below). An optically thin line is helpful to discriminate the infall signature from other possible interpretations of the double-peaked proles. Cloud complex
which has a couple of separate velocity components would also show a double-peaked line
prole, but the double peaks should appear in both optically thick and thin lines in sucha case. A set of double-peaked opaque lines and single-peaked transparent lines with ablue asymmetry is strong evidence for the infall motions.
We also present spectral proles convolved by an o-centered Gaussian beam by4000(6000AU) in Figure 6.2. The line proles are totally dierent from those convolvedby the centered beam. The peak intensity is weak even at the most active phase (Curve12). The lines remain narrow throughout the evolution and wings never come out. These
drastic changes in the line properties indicate that the broadening of lines and the appearance of wings are produced in a conned region near the center, and thus they are diluted
away in an o-centered beam. The blue asymmetry becomes weaker in an o-centered
beam for HCO+J =3-2 line, and it completely disappears from the other lines shown in
Figure 6.2. In general, an o-centered beam suppresses the infall signatures.
121

CHAPTER 6. MOLECULAR LINE SPECTRA

1
.2
1
.8
6
.15
12
.6
13
.1
.4

.05
.2

0

0
-2

0

2

.014

-2

0

2

-2

0

2

.0012

.012

.001

.01
.0008
.008
.0006
.006
.0004
.004
.0002

.002

0

0
-2

0

2

Fig. 6.2.| Spectral line proles convolved by an o-centered beam by 4000 for HCO+ andH13 CO+ (J =3-2 and 4-3) are shown. See also Figure 6.1.
122

6.2. CALCULATIONS AND RESULTS
We dene a line width, %V , as follows according to Zhou(1992).
%V = 2#2 ln 2(hV 2 i ; hV i2 )]1=2 (6.14)

where the angle bracket indicates an average, hqi  P(qTA)= P TA. A Gaussian line
prole yields %V equal to FWHM by this denition. Table 6.3 tabulates %V for eachline prole. Line widths increase monotonically with time for every transition, reecting
the increasing infall velocity.
For comparison with the results by Zhou(1992), we perform spectral synthesis calculations for CS lines using the same parameters as adopted by Zhou for the fractional
abundance (X = 4  10;9) and the resolution (6000 at 140 pc). We found that the largest
line widths throughout the evolution are 1.45, 2.01, and 2.58 km/s for CS J=1-0, 2-1, and
3-2 lines, respectively, at the end of the main accretion phase. As the dynamical evolution
resembles the LP solution rather than the Shu's Expansion-Wave collapse solution(x5.2),
our results should be compared with the \LP ow" in Zhou(1992). Zhou obtained %V =4.03, 4.24, and 4.43 km/s for CS J = 1-0, 2-1, and 3-2, respectively, at a corresponding
evolutionary stage (t = 1:5  105yr) for the \LP ow". The line widths obtained by
Zhou are signicantly larger than our results. Zhou excluded the LP ow because of theoverestimation of the line width, but our results do not suer from such unrealisticallybroad line proles. Indeed, our results for CS line proles (Fig. 6.3) are totally dierent
from Zhou's \LP ow" proles, which showed an extremely broad absorption at the line
center. As predicted in x1.2.3, very large %V found by Zhou(1992) is due to a large valueof the infall velocity (3.3cs) at large radii in the \naive" LP solution, while the velocity
declines to zero near the outer boundary in a \realistic" LP solution that we obtained,
which yields acceptable values for %V .
Another important result of our calculation is the existence of wings extending tov = 2km/s, in contrast to some past studies where wings could not be produced byinfall models (e.g., Gregersen et al. 1997). A main reason why the past calculations failed
to reproduce wings in observed line spectra is attributed to the use of the Expansion-Wavecollapse solution obtained by Shu(1977). The Expansion-Wave collapse solution consistsof an inner free-falling region and a surrounding hydrostatic envelope, where the latter
gives no contribution to the high velocity wing. Moreover, the Expansion-Wave collapse
solution has the minimum mass accretion rate (M_ = 0:975c3s =G at the origin) among thefamily of the isothermal self-similar solutions that include the singular origin. Hence thegrowth rate of the central point-mass (or a protostar) is so slow that the gravitational
potential grows only gradually high velocity wings are not produced.
As a consequence, we found that a radiation hydrodynamical model reproduces typicalfeatures in line spectra of class 0 and I sources more naturally than the isothermal self123

CHAPTER 6. MOLECULAR LINE SPECTRA

Table 6.3. Velocity Widths.
Molecule

TransitionCS

J =1-0J =2-1J =3-2
J =4-3J =1-0J =2-1J =3-2
J =4-3J =1-0J =2-1J =3-2
J =4-3HCO+
H13 CO+1%V a
6#km/s]
12

13

0.770
0.911
0.854
0.707
0.645
0.725
0.672
0.485
0.384
0.382
0.360
0.358

1.07
1.23
1.28
1.30
0.885
1.06
1.06
1.04
0.615
0.704
0.812
1.027

1.41
1.77
2.00
2.23
1.22
1.69
1.93
2.10
0.810
1.03
1.26
1.53

1.45
2.01
2.58
3.38
1.22
1.85
2.27
2.57
0.885
1.44
2.04
2.54

a The velocity width, %V , is dened by equation (6.14). Refer to Table 5.1 in x5.2 for label

numbers. See also Figures 6.1 and 6.3.

124

6.2. CALCULATIONS AND RESULTS

1.4

1.4

1.2

1

1

6

1.2

1

12
.8

.8
13

.6

.6

.4

.4

.2

.2

0

0
-4

-2

0

2

4

1.2

-4

-2

0

2

4

-4

-2

0

2

4

.8

1
.6
.8

.6

.4

.4
.2
.2

0

0
-4

-2

0

2

4

Fig. 6.3.| Computational results of spectral line synthesis for CS (J =1-0, 2-1, 3-2, and4-3) are shown. Adopted temperature, density, and velocity distributions are based on
the initially homogeneous model described in Chapter 5. Line proles are given for the
initial condition(labeled by '1'), just after the rst core is formed('6'), early ('12') and late('13') stages of the main accretion phase. Refer to Table 5.1 in x5.2 for label numbers.
125

CHAPTER 6. MOLECULAR LINE SPECTRA
similar solutions. Any of the self-similar solutions are found to be too idealized for detailed
line spectral synthesis.

1266.3. SUMMARY

6.3 Summary
We summarize the results obtained in this chapter.
1. We have developed a numerical code for non-LTE line transfer problems and applied it to protostar formation. The LVG and micro-turbulence approximations are
not used in our formulation, and arbitrary distributions for density, temperature,and velocity are applicable. For dynamical models, results by the radiation hydrodynamical calculations are used, where the density, temperature, and velocity
distributions are obtained in a self-consistent manner.
2. We conrm that the computational results show double-peaked proles in a stronger
blue peak for optically thick molecular lines. Optically thin lines show single-peakedproles with a slight blue asymmetry. These qualitative features are consistent withpast studies.
3. The peak is most intensive at an early stage of the evolution, which corresponds to
the genuine class 0 phase. While the peak intensity is lowered in later stages, line
widths increase monotonically with time. An o-centered beam is found to suppress
the infall signatures. On the contrary to the remarks by Zhou(1992), we did notnd overestimation of line widths despite the dynamical evolution resembles the LP
solution rather than Shu's Expansion-Wave collapse solution. Modest values of line
widths are obtained under a realistic boundary condition.
4. The infall motion also produces wings extending to v = 2km/s in line spectra.
Wings could not be produced in previous works where the Expansion-Wave collapse
solution was adopted as the infall model. Our results indicate that the presence ofwings in line spectra does not always assure the existence of outows.
Our results do not exclude possible contamination of outows in observed spectra.
Since 87% of class 0 and I sources show evidence of driving an outow (Mardones et al.
1995), the contamination would be inevitable. However, there is no theoretical investigation for contributions of outows to line spectral proles. Quantitative studies on thecontamination of outows are postponed until a reliable model for molecular outows isdeveloped.

127

CHAPTER 6. MOLECULAR LINE SPECTRA

128

Chapter 7
AN EVOLUTIONARY PICTURE
OF PROTOSTAR FORMATIONIn the preceding chapters we established a theoretical model for protostar formation
which accounts for recent observations. The theoretical treatment in this model is fully
consistent in spherical symmetry, and thus the computational results including the age,luminosity, and bolometric temperature are quantitatively reliable. On the other hand,the physical structure of actual YSOs is not in spherical symmetry because of the eectsby the angular momentum, the magnetic elds, and the bipolar outows. In this chapter
we discuss possible contributions of these non-spherically-symmetric eects to our results.

129

CHAPTER 7. AN EVOLUTIONARY PICTURE OF PROTOSTAR FORMATION
We pointed out in Chapter 5 that edge-on class I objects would be observed as class
0 sources (\class 0-like class I"), and they would contaminate the \genuine" class 0 objects which have the age of 2  104yr and are substantially younger than class I sources.
Similarly, evolved class I objects would be seen as at spectrum T Tauri stars dependingon the inclination angle to an observer.
Many theoretical and observational evidence support the idea that at spectrum T
Tauri stars are associated with circumstellar envelopes. Natta(1993) reproduced the observed at spectra by considering the backwarming of the disk by the scattering of stellar
photons in the envelope. Calvet et al.(1994) showed that the thermal emission from theextended envelope with a polar cavity forms the at spectra in addition to the emission from the central star and disk. Kikuchi, Nakamoto, & Ogochi(1999) performed
two-dimensional radiative transfer calculations in a fully self-consistent manner and successfully produced the at spectra by the backwarming of the disk by the scattering andthermal reprocessing of stellar photons in the envelope. Observations revealed the presence of extended envelopes surrounding HL Tau (Hayashi, Ohashi, & Miyama 1993) andT Tauri (Momose et al. 1996), which are typical at spectrum T Tauri stars. T Tauristars with at spectra tend to have larger extinction than typical T Tauri stars (Osterloh& Beckwith 1995). A recent high resolution observation (Close et al. 1997) succeeded indirect imaging of HL Tau in the near infrared. They conrmed the existence of bipolar
cavities, through which optical photons from the central star escape out.
Considering comprehensively our ndings and other theoretical and observational
evidence, we illustrate an evolutionary picture of protostar formation in Figure 7.1.
Protostellar evolution proceeds from top to bottom. Before t reaches t0 , where t0
2  105yr indicates the instance of the birth of a protostar (see Table 5.1 in x5.2), a cloudcore experiences gravitational collapse without any signicant changes in SED (x5.4).
Such cloud cores would be observed as star-less cores even after the formation of the rst
core (Chap. 4). After a protostar is born (t > t0 ), the appearance of a YSO depends
not only upon the evolutionary time but on the inclination angle to an observer, which isindicated by the horizontal axis. A pole-on view enables an observer to see more deeply
inside the protostellar envelope because the infalling material along the polar axis is swept
up by the centrifugal force and the bipolar outows. On the contrary, an edge-on viewsuers from enhanced column density in the disk plane, which cause signicant underestimation of the bolometric temperature (x5.4). A YSO at a certain evolutionary positionis therefore observed as if it were more evolved on a pole-on view and less evolved on anedge-on view. This tendency is schematically illustrated in Figure 7.1. Observed class
0 sources consist of genuine class 0 objects and edge-on class I, or class 0-like class I,
130t - t 0[105 yr]

Inclinationedge-onpole-on-2-1

M[M * ]
0.25

Starless Cores

0
0.2

Genuine Class 0
Molecular
Outflows(?)

Class 0-like

0.5Class I
0.7

1Class I

0.8

Flat Spectrum
T Tauri Stars
2

Fig. 7.1.| A Schematic picture of the evolution of protostar formation is illustrated.
The vertical axis represents the evolutionary time in the unit of 105yr, where the oset
value, t0 , indicates the instance of the birth of a protostar (see Table 5.1 in x5.2). The
protostellar mass, M , is shown based on Figure 5.7 in Chapter 5. The epoch in whichoutows associate with YSOs is also indicated. The horizontal axis corresponds to theinclination angle of the YSO to an observer. The area occupied by each type of YSOs
roughly corresponds to the number of the objects expected to be observed.
131

CHAPTER 7. AN EVOLUTIONARY PICTURE OF PROTOSTAR FORMATION
objects. A at spectrum T Tauri star would be evolved to a typical classical T Tauri staras the circumstellar envelope is dispersed or exhausted.
The area occupied by each type of YSOs in Figure 7.1 roughly corresponds to thenumber of the objects expected to be observed. On the basis of this fact, one can givea constraint on the epoch in which molecular outows associate with YSOs. Mardones
et al.(1997) listed 23 class 0 sources and 24 class I, where 23 sources of class 0 (100%)
and 18 sources of class I (75%) show evidence of outows. Therefore, the outow epochshould cover the whole region of genuine class 0 and class 0-like class I, and 3/4 of class
I in Figure 7.1. If an outow begins to be driven when the main accretion phase starts(i.e., t = t0 ) as predicted by Tomisaka(1998), an outow should continue being driven
until t ; t0 1:3  105yr to naturally reproduce the observations. This estimation is, ofcourse, very crude because it assumes that outows are driven universally and steadily inthe same epoch regardless of the individuality of each object. More elaborate estimation
requires a further progress in both observational and theoretical studies in the driving
mechanism of outows.
Class 0 sources are observed typically on edge-on views (Hirano et al. 1997). This fact
is attributed to the contamination of edge-on class I (or class 0-like class I) in our picture.
The inclination angle would be randomly distributed for genuine class 0 objects, but thecontribution of edge-on class I shifts the average inclination angle to a higher (i.e., moreedge-on) value.
A genuine class 0 object on a pole-on view may be seen as a class I source. Such an object would be very luminous (25L, see x5.3) and have a small protostellar mass(0:1M ).
Momose et al.(1998) derived the protostellar mass as 0:1M for L1551 IRS5, whose
luminosity is 30L (Cohen et al. 1984). Hence L1551 IRS5 may be a genuine class 0
object on a pole-on view, or with a polar cavity of a large opening angle, which enables
infrared photons to escape out so that IRS5 is observed as a class I object despite its
extreme youth.
A more quantitative picture than Figure 7.1 should involve numerical calculations
based on multi-dimensional radiation hydrodynamics. Such calculations require verylarge computational costs and is then not available at present (x1.1.2), but it will bea challenging task in future.

132

Chapter 8
SUMMARY AND CONCLUSIONSOur ndings throughout this paper are summarized as follows.
1. In order to construct a theoretical model for protostar formation which is capable
to account for recent observations, we have developed a numerical code for radiation hydrodynamic calculations which includes all the necessary physical processes,
such as dissociation of molecules, ionization of atoms, non-ideal equation of state,adequate opacities for gas and dust, and so forth. The numerical code is designed
to yield the evolution of spectral energy distributions as well as the dynamical evolution. We established in this work a theoretical model for the whole evolution of
protostar formation in a consistent scheme which is accountable for observations.
No such studies have been done to date.
2. First, we performed numerical calculations with simplied input physics for clarifying the evolution before the central protostar is formed. The rst core nally has
the mass of 0:05M and the radius of 5 AU for Tinit = 10K and the typical
opacity. If the initial temperature and opacity are given, these values for mass andradius are uniquely determined and are independent both of the mass of the parent
cloud core and of the initial density prole.
3. We have found the analytic expression for the radius of the rst core and the accretion luminosity as functions of the central density. They are weakly dependent on
the initial temperature and opacity. Those expressions are useful for understanding
the numerical results.
4. The shape of the spectral energy distributions (SEDs) does not reect the central
thermal evolution of collapsing cloud cores throughout the rst collapse. This result,of course, does not exclude the possibility that the rst cores might be observed on
the pole-on view, where the envelope material would actually be more depleted than133

CHAPTER 8. SUMMARY AND CONCLUSIONS
predicted by spherically symmetrical calculations. The luminosity increases up to
0:1L at the end of the rst collapse phase.
5. Dense cloud cores with luminosities higher than 0:1L possibly have young protostars inside the opaque envelopes, even if observations show no evidence for protostellar activities.
6. Following the study on the rst collapse, an RHD numerical calculation is carriedout to reveal physical processes in the formation of a 1M protostar. The wholeevolution are pursued from the beginning of the rst collapse to the end of the mainaccretion phase. All the physical phenomena that appear in protostellar formation,(e.g., the dynamical collapse of a cloud core, the formation of the pressure-supported
core, the accretion shock front on the protostellar surface, and the supersonicallyinfalling gas) are incorporated in a self-consistent manner in the calculation.
7. For the initially homogeneous model, the accretion luminosity rapidly rises to the
maximum value of 25L just after the formation of a protostar, and declines gradually as the mass accretion rate decreases. In contrast, the luminosity increases
monotonically with time for the initially hydrostatic model. This discrepancy arises
because the mass accretion rate varies depending inward acceleration is zero or
non-zero in the initial condition, which aects the luminosity curve.
8. We conrm that the SED evolves from a 10K greybody spectrum to hotter spectra typical for class I and II sources. The bolometric temperature derived by theobtained SED shows that cold sources with Tbol less than 70K are younger than2  104yr. This indicates that class 0 objects have younger ages by an order ofmagnitude than class I objects. However, Tbol can be signicantly colder even aftert 105 when near infrared components are hidden by locally enhanced column density as in the disk plane. Observed class 0 sources should be the compound of the\genuine" class 0 that is as young as 104 yr and more evolved protostars on edge-on
view. The contamination of older protostars are not negligible because they are
intrinsically abundant than genuine class 0 objects.
9. Since observations indicate that class 0 sources is typically more luminous thanclass I sources, we exclude the initially hydrostatic model where the luminosity
increases monotonically with time. The initially homogeneous model is found toshow the tendencies consistent with observations. Star formation therefore tends tobe triggered under initial conditions that are not in hydrostatic equilibrium. The
decreasing luminosity with time implies that signicantly luminous and cold YSOs
134

( 1L and Tbol < 70K) should be genuine class 0 objects. Less luminous YSOs
are more evolved protostars on edge-on view even if they have cold Tbol.
10. We have also developed a numerical code for non-LTE line transfer problems andapplied it to protostar formation. The LVG and micro-turbulence approximations
are not used in our formulation, and arbitrary distributions for density, temperature, and velocity are applicable. For dynamical models, results by the radiationhydrodynamical calculations are used, where the density, temperature, and velocity
distributions are obtained in a self-consistent manner.
11. We conrm that the computational results show double-peaked proles in a stronger
blue peak for optically thick molecular lines. Optically thin lines show single-peakedproles with a slight blue asymmetry. These qualitative features are consistent withpast studies.
12. The peak is most intensive at an early stage of the evolution, which corresponds to
the genuine class 0 phase. While the peak intensity is lowered in later stages, line
widths increase monotonically with time. An o-centered beam is found to suppress
the infall signatures. On the contrary to the remarks by Zhou(1992), we did notnd overestimation of line widths despite the dynamical evolution resembles the LP
solution rather than Shu's Expansion-Wave collapse solution. Modest values of line
widths are obtained under a realistic boundary condition.
13. The infall motion also produces wings extending to v = 2km/s in line spectra.
Wings could not be produced in previous works where the Expansion-Wave collapse
solution was adopted as the infall model. Our results indicate that the presence ofwings in line spectra does not always assure the existence of outows.
14. Finally, we illustrated an evolutionary picture of protostar formation considering ourndings in this paper and other theoretical and observational evidence. In termsof the evolutionary time and the inclination to an observer, protostellar objects arefound to be successfully categorized.

135

CHAPTER 8. SUMMARY AND CONCLUSIONS

136

Appendix A
Test Calculations for Radiation
To check the abilities of the moment equations solver, the transfer equation solver, and theconvergence between them, test calculations for the numerical code are presented brieyin this section. Because the moment equations of radiation are dicult to solve analytically when they are coupled with hydrodynamics, eects of radiation are investigatedseparately from hydrodynamics here (e.g., Ensman 1994). Test calculations for hydrodynamics are not presented in this paper, but we have checked that numerical results ofsome hydrodynamical tests show good agreements with analytical predictions.
For radiation tests, material is assumed to be static, with the RHS of the equationof motion set to be zero. For simplicity, the density distribution is assumed to be homogeneous and radiation to be grey. Then, suppose small uctuations of temperaturein radiative equilibrium. The uctuation should be smoothed out by radiation quasistatically if the radiation energy density is suciently small compared to the internalenergy density of material. This mode is called the \thermal relaxation mode" of radiation (Chap. 8 in MM). The relaxation rate n(k) for a perturbation with wavenumber kis (Spiegel,1957)
wheren(k) =  #1 ; k0 cot;1 k0 ]
=

16

3
0 T0

(A.1)

(A.2)
cv 
where cv is the specic heat for material and subscript 0 represents unperturbed quantities.
Spiegel(1957) derived equation (A.1) for a perturbation of plane wave, i.e., eik x.
However a spherical wave solution is required for the test of our numerical scheme. We
have found that the 0-th order spherical Bessel function, j0 (r)  sin kr=kr, is the eigenfunction in the spherical case, and it yields the exactly same dispersion relation as equation(A.1). See Appendix B for derivation.
137

APPENDIX A. TEST CALCULATIONS FOR RADIATION
Figure A.1 shows relaxation of a uctuation of the spherical wave j0(kr) in the case of
the total optical depth of 10. The unperturbed temperature distribution is homogeneousat 10K . A temperature disturbance (Figure A.1a) damps toward the unperturbed state.
The radiation temperature Trad  (Erad=a) 41 (Figure A.1b) damps correspondingly tothe relaxation of the material temperature, but slightly more rapidly. For more opaquematerial, of course, the dierence between T and Trad would be smaller.
The dispersion relations are shown in Figure A.2. The ordinate represents n(k)t ,where t  (c );1 is the mean ight time of a photon. The ordinate therefore means how
much a temperature disturbance decays while a photon travels over a mean free path.
First, to test the moment equations solver solely, we calculate the thermal relaxation
mode in the Eddington approximation, xing fE  1=3 instead of solving the transfer
equation solver. In this case equation (A.1) reduces to the dispersion relation under the
Eddington approximation (Unno & Spiegel, 1966),

nE (k) = 1 + 3( 0 )2 :
k

(A.3)

The numerical results show good agreement with equation (A.3) in Figure A.2a. Therefore
the moment equations are found to be solved correctly.
Second, accuracy in the transfer equation solver is tested. In Figure A.2b, the numerical result which was obtained using the Eddington factors calculated by the transfer
equation solver is plotted with the analytical dispersion relation (eq. #A.1]). They agree,
again, with each other within a relative error of one percent. These results guarantee theaccuracy of the method over a wide range in optical depth.

138

T [K]

(a)

(b)

10.1

10.1

10.05

10.05

10

10

0

2000

4000 6000
R [AU]

.0002

8000

10000

0

2000

4000 6000
R [AU]

8000

10000

(c)

.00015

.0001

0

0

2000

4000 6000
R [AU]

8000

10000

Fig. A.1.| Results of test calculations for the thermal relaxation mode. (a) Relaxation ofa small disturbance of temperature. The given perturbation (the curve of the largest amplitude) decays to a homogeneous distribution. (b) Relaxation of radiation temperature.
(c) Relaxation of radiative ux.
139

APPENDIX A. TEST CALCULATIONS FOR RADIATION

-5

-6

-7

-8

-9
-2

-1

0

1

2

-2

-1

0

1

2

-5

-6

-7

-8

-9

Fig. A.2.| The dispersion relation for the thermal relaxation mode of radiation. See
text for the physical meanings of the axises. Upper: The dispersion relation under the
Eddington approximation. Solid line is the analytical solution and open circles denotethe numerical results with f  1=3. Lower: The exact dispersion relation without the
Eddington approximation. Dotted line corresponds to the exact solution and open circles
140 the variable Eddington factor.
denote the numerical results obtained by using

Appendix B
Thermal Relaxation Mode in
Spherical Symmetry
The thermal relaxation mode of a radiating uid was investigated rstly by Spiegel(1957).
He derived the dispersion relation (A.1) for a perturbation of a plane wave, which inevitably requires that the temperature distribution must be in slab symmetry. We willshow below that a temperature disturbance in spherical symmetry decays in the samerate as given by equation (A.1).
We assume that there are no uid motions and no heat sources, and the radiation eldis grey. The uid energy equation (2.3), then, reduces to @e@t = 4 (J ; B )(B.1)

where J  (1=4 ) Id) = (c=4 )E is the mean intensity and is a grey opacity per unit
volume. The source function in LTE, B , can be replaced by ( = )T 4.
The specic intensity I is determined by the the transfer equation,
H@I = ; (I ; B )(B.2)
@swhere we have neglected the transit time of radiation.
For linearizing the above two equations, we write the temperature as follows, neglecting
the terms of higher order uctuation.
T (r t) = T0(r) + T1 (r t):(B.3)

The other quantities are written as0
B (r t) = B0 (r) + @B
@T T1  B0 + B1 

141

(B.4)

APPENDIX B. THERMAL RELAXATION MODE IN SPHERICAL SYMMETRY(r t) = 0 (r) + @ 0 T1  0 + 1 :(B.5)@TUsing equations (B.3),(B.4), and (B.5), equations (B.1) and (B.2) are reduced to thefollowing equation, which describes the relaxation of the temperature disturbance (Chap.8
in MM).
@T1 (r) = ; #T (r) ; 1 I d) Z 1 dsT (r ; ns)e; 0s]
(B.6)
1
0
1
@t
4
0
where n is the unit direction vector along a light ray and  is dened in equation (A.2).
The terms including 1 have vanished because of the assumptions that the unpurterbed
radiation eld is isotropic and in radiative equilibrium. All the assumptions introduced
above are the same as those by Spiegel(1957).
Equation (B.6) admits separable solutions of the following form for spherical symmetriccases,
T1(r t) = (k t)j0(kr)
(B.7)
where j0 (x)  sin x=x is the 0-th order spherical Bessel function.
Substituting equations (B.7) and r0  jr ; nsj into (B.6) rewritten in the cylindricalcoordinates (r ), where  is the direction cosine,@ = ; #1 ; 1 Z 1 d Z 1 de; j0(kr0) ]
(B.8)@t
2 ;1 0j0(kr)
where   0 s.
For the purpose of replacing r0 with  and , let the path of a light ray s written down
explicitly (see Figure B.1),

s = (r02 ; r2(1 ; 2)) 12 + r:(B.9)
The negative sign in the rst term in RHS is required only for r0 < r. After rearranging
equation (B.9) and replacing s with  , we have
(B.10)
  kr0 = k ( 2 ; 2 r + ( r)2) 12 
0

0

and

0

d = ; rk0 2 d

;1   ( = ;1) = k ( + 0 r) 1   ( = 1) = k j ; 0rj:
00After integrating over  instead of , the last term in equation (B.8) equals
Z 1e; #cos  ; cos  ]=j (kr)
0
d
1
;1 0
2rk2 0

142

=

0

Z 1

k 0

sin( k0  )
;

de
:



The last integral is reduced to cot;1( 0 =k). Hence, equation (B.8) becomes@ = ; #1 ; 0 cot;1 ( 0 )]
@t
kk
The thermal relaxation rate n(k) is dened as@ = ;n(k)
@tComparing equations (B.11) and (B.12), we have equation (A.1).

143(B.11)(B.12)

APPENDIX B. THERMAL RELAXATION MODE IN SPHERICAL SYMMETRY

r’ > rr’ < r

s
s+
s-

s
θ

θ
r

r’

r’

Fig. B.1.| An illustration of the path of integration, s.

144rAppendix C
Conditions for Violation of
Isothermality
The assumption that the early stage of protostellar collapse proceeds isothermally is
widely accepted in studies of star formation. This isothermality is justied because thethermal emission by dust grains is quite eective in the early evolution and the released
gravitational energy can be immediately radiated away. The isothermal evolution, however, ceases when thermal radiation no longer cools the cloud against the compressional
heating. Detailed studies in thermal processes which violate isothermality involve complicated problems including radiative transfer, and hence quantitative investigations havenot been done to date.
On the other hand, recent progress in computational facilities has enabled us to handle
such complicated problems rather easily. In our previous work (MMI) numerical calculations for protostellar collapse were carried out with an exact treatment for radiativetransfer. MMI showed that the condition with which isothermality is violated is classiedinto three dierent criteria, as reviewed in xC.1 below. MMI found, both by numericalresults and analytical estimates, that small dierences in the cloud temperature and opacity cause drastic changes in the critical density, crit, which is the central density of thecollapsing cloud when the isothermal evolution is terminated.
We generalize the analysis of MMI in the present paper and reach a conclusion that thecondition of \  1" never terminates isothermality in possible ranges of parameters for
actual molecular clouds. This result contrasts with the familiar idea that isothermality isviolated when the cloud becomes opaque to its thermal radiation at crit 10;13g  cm;3 .
It is a great necessity to critically reexamine this commonly believed idea, which has been
often stated in the literature (L69 Appenzellar & Tscharnuter 1974 WN).
In the context of star formation, it is important to determine when the isothermal evolution ceases. In spherical collapse of a preprotostellar core, the violation of isothermality
145

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY
leads the formation of a central adiabatic core (the so-called rst core), whose size andmass are characterized by crit(MMI). This critical density has more direct importance in
cylindrical collapse because an isothermally collapsing lament is expected to fragment
due to hardening of the equation of state (Inutsuka & Miyama 1992, 1997, hereafter IM92
and IM97, respectively). IM97 derived crit and the mass of a fragmented core as functions of the initial temperature and opacity. An important implication by IM97 is thatthe minimum Jeans mass for fragmentation depends sensitively upon Tinit and , on thecontrary to the conclusions by previous authors (Low & Lynden-Bell 1976 Rees 1976
Silk 1977 Boss 1988). We will revisit this topic in more detail in xC.2.
Multi-dimensional numerical calculations are carried out commonly under the assumption of an isothermal (or polytropic) equation of state. This is because the radiativetransfer equation requires unacceptable computational eort to solve exactly in multidimensions. Some numerical studies have shown that an isothermally collapsing cloud insome cases yields a barred structure and/or fragments in the central region (Truelove etal. 1998 and references therein). These results, however, may depend on whether theisothermal approximation remains valid or not throughout the evolution, especially in thecentral region where the density can be very high. Furthermore, Truelove et al.(1997,1998) pointed out that articial fragments are observed in hydrodynamic calculations
when cell size is insucient to resolve the local Jeans length, which implies that one
should consider as well the possibility that physical fragment may not actually occur before isothermality is broken down. Our criteria for violation of isothermality are useful forchecking validity of the isothermal approximation. Moreover, one can harden the equationof state (EOS) when the density reaches a certain critical value in numerical calculations
for collapsing clouds (e.g., Burkert, Bate, & Bodenheimer 1997). Our results would be
helpful to such a type of calculations, in which the isothermal EOS is employed where < crit and a more harder EOS takes the place elsewhere.
In xC.1, we review crit dened by MMI and apply it to more general values of thecloud temperature and opacity. In xC.2 the minimum Jeans mass for fragmentation isreexamined on the basis of crit obtained in xC.1. Summary and future tasks are described
in xC.3. Contents of this section are to be published in Masunaga & Inutsuka(1999).

146

C.1. DERIVATIONS FOR THE CRITICAL DENSITIES

C.1 Derivations for the Critical Densities
In general, molecular clouds are under thermal balance described by the following relation.
;g + ;ext = "th

(C.1)

where ;g , ;ext , and "th are the compressional heating rate of gas, the heating rate dueto the external sources such as cosmic rays and visual or UV photons from surroundingstars, and the radiative cooling rate, respectively. The isothermal collapse proceeds during;g "th and thus ;ext  "th, but isothermality is broken down when the compressionof gas becomes so eective as to heat up the cloud suciently against the cooling.
In this section we discuss on the critical central density crit when the isothermal
evolution is terminated in gravitational collapse, following the manner developed by IM97
and MMI. MMI found that the condition with which isothermality is broken down isclassied into three dierent criteria as follows.
Case 1. Isothermality is violated when ;g = "th, if ;g overwhelms "th before the opticaldepth of the collapsing cloud core,  , reaches unity. This critical density is denotedby =; .
Case 2. For  > 1, the cooling rate is given by the energy transport rate due to radiativediusion, "dif , instead of "th. If "th > ;g during  < 1 and "dif < ;g when   1,
the isothermal evolution ceases when  arrives at unity. The central density when  1 is denoted by  1.
Case 3. If "dif > ;g when   1, isothermality survives even after  exceeds unity until
the central density reaches a critical value dif , which denes the central density
when ;g becomes comparable to "dif .
Now we evaluate ;g , "th, and "dif to derive crit . Note that the heating and cooling rates are dened per unit mass. For brevity we suppose the local thermodynamical
equilibrium (LTE), which admits "th to be described simply as follows.
4 "th = 4 (Tinit) Tinit(C.2)

where is the frequency-averaged opacity 1 per unit mass, which is independent of density
due to LTE, and denotes Stefan-Boltzmann constant. The temperature is kept at Tinit
1 In equation (C.2) the frequency-averaged opacity  should be replaced by the Planck mean opacity.

In contrast, the Rosseland mean opacity may be preferred in equations (C.10) and (C.11) below. We,
however, do not discriminate them for brevity in this paper because dierence in values of these two
opacities is small.

147

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITYat the initial isothermal stage. The compressional heating rate for gravitational collapseis (see x5.2. in MMI)q
(C.3)
;g = Ac2s 4 G
where cs is isothermal sonic speed and G is the gravitational constant. A numericalconstant A is found to be order of unity and nearly constant through the evolution (c.f.,
IM97 and MMI).
The optical depth,  , of the collapsing cloud is dened by

=

Z 1

0

dr:

(C.4)

In order to evaluate the optical depth, one must suppose a density structure of the collapsing cloud. For a spherical cloud, we employ the isothermal similarity solution obtainedby L69 and Penston(1969). The density distribution in their solution is approximatedby  =constant for r < Rb and  / r;2 for r > Rb, where the boundary radius, Rb , isapproximately the Jeans length corresponding to the current central density, c . Hence
Rb is expressed as
cs 
Rb = C p42 G
(C.5)
cwhere a dimensionless constant C is 0:75 by comparison with numerical results(MMI).
Equation (C.4) is then reduced to
Z Rb
Z 1

cdr +c ( Rr );2dr = 2 cRb :
(C.6)
0
Rb
b
Combining equations (C.5) and (C.6),
cs :
(C.7)
 = C c p44 G
c
For a lamentary cloud we use equation (C.8) following the formulation by IM97,
where the isothermal equilibrium conguration is adopted for the density distribution.

 = 4 c

s

2c2s
Gc

(C.8)

(eq. #28] in IM97). This approximation is justied because an isothermal lament inequilibrium can collapse preserving the initial conguration, i.e., undergo homologous
collapse. Moreover, a lament formed by fragmentation of sheet-like clouds has a line mass
only twice as large as the equilibrium line mass (Miyama, Narita, & Hayashi 1987a,b),
and therefore the equilibrium structure is applicable also to collapsing clouds. One cansee that equations (C.7) and (C.8) accord with each other except for numerical factors.
148

C.1. DERIVATIONS FOR THE CRITICAL DENSITIESThe opacity per unit mass is assumed to be constant throughout the collapsing cloud.
This assumption is well justied when considering the dust continuum opacity, which is
independent of velocity structure of clouds in contrast to the atomic or molecular line
opacity.
We dene the r-parameter, which is the ratio of the radiation energy density to thematerial internal energy, for  = crit as follows,4
E = 4( ; 1) Tinit
rcrit  ec c2 :
crit s(C.9)

The energy transport rate due to radiative diusion is(MMI)4
E = 4 (Tinit) Tinit(C.10)
tdif
2when the internal energy of material dominates the radiation energy density (i.e., rcrit <1), and
2 (T )cec"dif RD= s init 2(C.11)t( ; 1)

"dif MD

difwhen the radiation energy density dominates (i.e., rcrit > 1), or collectively"dif = min("dif MD "dif RD)

(C.12)

where e and E are the internal energy density of uid per unit mass and the radiationenergy density per unit volume, respectively, and  denotes ratio of the specic heats (weset  = 5=3 in this paper). The radiative diusion time tdif is dened by  2 p=c, wherep  1=  is the mean free path of a photon and c is the speed of light.
Combining equations (C.2)-(C.10), the critical densities for the violation of isothermality are derived. For spherical symmetry,

Tinit )6+2 =; = 4:7  10;15g  cm;3 ( 0:01cm02  g;1 )2 ( 10K

(C.13)

Tinit );1;2  1 = 2:6  10;13g  cm;3 ( 0:01cm02  g;1 );2( 10K(C.14)

Tinit )(4;2 )=3 dif MD = 6:9  10;14g  cm;3 ( 0:01cm02  g;1 );2=3 ( 10K(C.15)

which correspond to equations (14), (17), and (20) in MMI. Here the opacity is approximated as(Tinit) = 0 ( Tinit ) :(C.16)
10K
149

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY
Note that in the typical range of temperature for molecular clouds  in equation (C.16)
coincides  , which is the power law index in frequency (i.e.,  /   ) for the dust opacity(e.g., Beckwith et al. 1990) and is typically 1-2 for interstellar dust grains. The
normalizing factor, 0 , reects metallicity and is 0:01cm2  g;1 for nearby molecular
clouds.
MMI found that these three critical densities are useful to predict when the isothermal evolution is terminated, comparing them with radiation hydrodynamic numericalcalculations with an exact treatment for radiative transfer in spherical symmetry.
In case of rcrit > 1, equation (C.15) above is replaced by the following equation, which
is derived using equation (C.11).
Tinit );2;2 :dif RD = 1:7  10;2g  cm;3 ( 0:01cm02  g;1 );2 ( 10K(C.17)
Note that spherical collapse is decelerated when the equation of state is harder than = 4=3, where  is the adiabatic component, i.e., d ln P=d ln . Hence the violation ofisothermality does not necessarily cause immediate inuence on the evolution(see MMI).
In contrast, cylindrical collapse should be decelerated even by slight hardening of theisothermal equation of state. Furthermore, IM92 and IM97 have shown that a collapsinglament is expected to fragment if the radial collapse is decelerated. Therefore crit forcylindrical collapse, which is considered below, has more signicance than the sphericalcase.
The critical densities for cylindrical (lamentary) collapse were obtained by IM97 as
Tinit )6+2 (C.18)
=; = 4:7  10;15g  cm;3 ( 0:01cm02  g;1 )2 ( 10K
Tinit );1;2 : 1 = 4:7  10;12g  cm;3 ( 0:01cm02  g;1 );2( 10(C.19)
K
Another critical density dif , which was not considered in IM97, is estimated as
Tinit )(4;2 )=3(C.20)
dif MD = 4:7  10;13 g  cm;3 ( 0:01cm02  g;1 );2=3 ( 10K
andTinit );2;2 :dif RD = 5:4g  cm;3 ( 0:01cm02  g;1 );2( 10K(C.21)
Dierence in crit between spherical and cylindrical cases appears only in numerical factors.
In the above equations we assumed that the gas constant k=mH  c2s =T is 3:6  107 erg=g K.
Equations (C.13)-(C.15), (C.17), and (C.18)-(C.21) show that the critical densities are
determined solely by Tinit and . Thus the criteria for the violation of isothermality can be150

C.1. DERIVATIONS FOR THE CRITICAL DENSITIES
classied by Tinit and . Figure C.1 depicts the critical density (shown by contour lines)
as a function of Tinit and . We have chosen the range of Tinit to be less than 1000K, below
which dust grains survive. Although this range is still wider than the possible temperature
range for actual molecular clouds, it is helpful to realize behavior of crit more clearly on
the Tinit - plain. The vertical axis, , is a product of 0 and the temperature dependence
according to equation (C.16). Under a given value of Tinit, the vertical axis is an indicatorof 0 , and hence, of metallicity.
For rcrit < 1, boundaries between Case 1-3 (=; =  1 and dif =  1) are found to
fall onto an equivalent condition, which is delineated by a dashed line in Figure C.1. (Thereason is described in the next paragraph.) The critical density is given by =; below thedashed line and crit = dif above. The region for crit = dif is divided into two domains:crit = dif RD and crit = dif MD. The condition of crit = dif RD, however, requires
too large Tinit and compared to the typical values for molecular clouds and therefore
the domain for crit = dif RD does not appear in Figure C.1. The shaded domain in
the lower-left corner, where the isothermality is broken before  reaches unity, indicatesaT 4 > e (a  4 =c), which, however, does not necessarily mean that the radiationenergy dominates the uid internal energy. Unless the external radiation eld is intensive
enough, E is much less than aT 4 in optically thin media. The model parameters adopted
by L69 and WN are indicated by stellated and triangular points, respectively, in Figure
C.1. The both points read crit 10;13g  cm;3 , which accords with the central density
when the isothermal evolution terminates in their numerical calculations. These points
accidentally locate near the dashed line, on which the isothermal evolution ceases when
 = 1. This seems a main reason why many authors have been misled by the idea that
opaqueness terminates isothermality.
Equations (C.2) and (C.10) shows that "th and "dif MD are smoothly connected witheach other at  = 1. In other words, =; =  1 and dif MD =  1 degenerate into an
identical line in the Tinit- plane (the dashed line in Figure C.1). Figure C.2a schematically shows the thermal evolution of a gravitationally collapsing cloud. If rcrit < 1, theisothermal evolution ceases when ;g reaches "th (which corresponds to Case 1 above), or
when ;g reaches "dif (Case 3). Although neither "th nor "dif provides exact estimates forintermediate optical depth, the critical density of  1 (Case 2) clearly plays no essential
roles when rcrit < 1.
In contrast to "dif MD, "dif RD are not ensured to be connected smoothly with "th near  1, as found by comparing equations (C.2) and (C.11). Therefore, =; =  1 anddif RD =  1 should be indicated by two separate lines in the Tinit - plane, and thedomain for crit =  1 would appear between these two lines. Figure C.2b draws the151

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY

1

-14
0

-12
-1

-14

-2

-10

-18

-3

-

-22

-4

-5
0

.5

1

1.5

2

2.5

3

Fig. C.1.| The critical density for the violation of isothermality (eqs. #C.18]-#C.21]) is
shown as contour curves. Note that the vertical axis represents not 0 but (Tinit) (see eq.
#C.16]). Contour levels begin with log10 crit #g  cm;3 ] = ;24 in the lower-left corner, andincrease as ;22 ;20 ;18   . The dashed line delineates =; =  1 and, equivalently,dif =  1 (see text for detail). The critical density is given by =; below the dashedline and by dif above. Shaded domains indicate rcrit > 1, where rcrit is the ratio of theradiation energy density to the material internal energy density, and rcrit < 1 elsewhere.
Stellated and triangular points correspond to the model parameters adopted by L69 and
WN, respectively.
152

C.1. DERIVATIONS FOR THE CRITICAL DENSITIESa) r crit < 1
Case 1.
Λ,Γ

τ<<1
Λth

Case 3.
Γ gas

Λdif
τ>>1

Γ gas
Λth = Γ gasΛdif= Γ gastb) r crit > 1
Case 1.
Λ,ΓCase 2.

τ<<1
Λth
Γ gas

Γ gas

Case 3.

Λdif

Γ gas

τ>>1
τ∼1

Λth = Γ gasΛdif = Γ gastFig. C.2.| A schematic picture which illustrates the classication for crit is shown. The
heating and cooling rates (;g , "th , and "dif ) are depicted as a function of time. The opticaldepth for thermal radiation,  , is small enough initially but increases monotonically as
the collapse proceeds. The upper panel is for the case of rcrit < 1, and the lower panel for
rcrit > 1. See text for more detail.
153

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY
thermal evolution for rcrit > 1, where Case 2 appears as well as Case 1 and 3. However,dif RD =  1 corresponds to Tinit = 1:1  1013 K(!) independently of , and hence Case
2 is highly implausible to occur.
As a conclusion, the critical density for the violation of isothermality is determinedby =; or dif MD, depending upon Tinit and , for gravitationally collapsing clouds. The
critical density of  1 has no signicance in practice. This fact has been overlooked
in all previous works on gravitational collapse of protostellar clouds including our own
papers(IM97 and MMI).

154

C.2. THE MINIMUM JEANS MASS FOR FRAGMENTATION

C.2 The Minimum Jeans Mass for Fragmentation
This section is devoted to discussions on the minimum Jeans mass of a fragment ofclouds. Previous formulations of the minimum Jeans mass are reconsidered and we showthat they should be modied on the basis of our newly derived criteria for the violationof isothermality.

C.2.1 Comments on the Previous Studies
The minimum Jeans mass for fragmentation, MF , is determined by the condition thatthe equation of state becomes suciently hard to stop further fragmentation. In order toobtain MF, Low & Lynden-Bell(1976), Silk(1977), and Kanjilal & Basu(1992) imposed two
conditions of ;g = "th (eq. #C.3]) and  = 1 simultaneously. As a result the temperatureof clouds that denes MF is determined uniquely as well as MF itself. However, as IM97
noted, Tinit as well as must be left as a free parameter because the cloud temperature
is actually given by environmental factors including the intensity of interstellar radiationand the rate of cosmic ray heating. Hence we estimate the minimum Jeans mass as afunction of both Tinit and , on the basis of the discussion on crit in xC above. In thismanner we propose a more generalized formulation for MF in xC.2.2.
MMI shows that the temperature increases only gradually in Case 1 after ;g reaches"th until the temperature eventually begins to rise rapidly when  exceeds unity. Therefore
one might wonder whether the condition of ;g = "th can actually trigger the fragmentation of spherical clouds.
Tohline(1980) raised serious objections to the previous models for fragmentation of
spherically collapsing clouds. His results indicate that the fragmentation occurs after acloud is collapsed into a highly attened sheet. A collapsed sheet-like cloud is subject
to fragmentation into laments (Miyama, Narita, & Hayashi 1987a,b), and therefore theproblem of cloud fragmentation becomes a problem of the evolution of laments. As
mentioned in the previous section, a collapsing lament is expected to fragment by aslight hardening of the isothermal equation of state in contrast to spherical collapse, andhence even a gradual increase of temperature can trigger fragmentation(IM92, IM97).
Rees(1976) considered the condition for determining MF as

GMF2 4 fR2 T 4 (C.22)
init
Rt
where t is the free fall time and f is a non-dimensional factor less than unity, which
depends on the detailed physics of the cooling and opacity. In contrast to the otherauthors mentioned above, Rees(1976) considered the thermal balance in collapsing cloud155

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY
after  exceeds unity, and thus his criterion (eq. #C.22]) corresponds to our Case 3.
This criterion, however, assumes that the released gravitational energy is instantaneously
radiated away from the cloud surface. In other words, equation (C.22) neglects the time
which is required to convey the released energy to the cloud surface by radiative diusion.
However, the energy transport rate due to radiative diusion is essentially important inCase 3 as described in the previous section. Therefore RHS of equation (C.22) should bemodied into
Z
4
2 T4
4
T
4
R
init
init 
3"dif dV
R(C.23)
2
where we used equation (C.10) and replaced R with  . Equation (C.23) claims that
RHS in equation (C.22) overestimates the cooling rate by a factor of  . It may be possible
to substitute 1= for f in equation (C.22), but in this case f cannot be a constant factor
since  varies with  (or t ) and R. Furthermore,  is a function also of both opacity andtemperature (because Tinit characterizes the similarity solution for density structure), andtherefore the original formulation by Rees (eq. #C.22]) requires a great care to use.

C.2.2 A Revised Denition of the Minimum Jeans Mass
Following IM97,qwe evaluate MF as the mean clump mass of a lament with the eective
radius of Hf = 2c2s = Gcrit, setting the mean separation between clumps to be 8  Hf :
s

16c3s
2 :
(C.24)MF ' crit f=f
G Gcrit
Eliminating crit in equation (C.24) by equations (C.18) and (C.20), we have
Tinit );(3+2 )=2 MFth = 3:7  10;2M ( 0:01cm02  g;1 );1 ( 10(C.25)
K
for Case 1,
Tinit )(5+2 )=6 (C.26)MFdif MD = 3:7  10;3M ( 0:01cm02  g;1 )1=3 ( 10K
for Case 3 with rcrit < 1, andTinit )(5+2 )=2 MFdif RD = 1:1  10;9M( 0:01cm02  g;1 )( 10(C.27)
K
for Case 3 with rcrit > 1.
Figure C.3 illustrates MF by contour lines in the Tinit - plane. The dashed line, tracing
the bottom of a ravine drawn by contour, represents MFth = MF 1 and MFdif MD =MF 1 similarly to Figure C.1. Then, MF equals MFth below the dashed line and MF =H 2  8H

156

C.2. THE MINIMUM JEANS MASS FOR FRAGMENTATION
MFdif MD above according to Figure C.1. Recalling that Low & Lynden-Bell(1976) andSilk(1977) considered crit = =; and crit =  1 simultaneously, one nds that theirsolution for MF is constrained on the dashed line in Figure C.3. Using equations (C.25)
and (C.26), one obtains MF dened on the dashed line:
MF = 8:4  10;3M( 0:01cm02  g;1 );1=(7+4 ) (C.28)

which indeed recovers MF dened by Low & Lynden-Bell(1976) and Silk(1977). Low &Lynden-Bell emphasized the weak dependence of MF on 0 as seen in equation (C.28).
On the other hand, if assuming that Tinit does not vary signicantly as is the case inactual molecular clouds, equation (C.25) shows that the dependence on 0 is no longer
weak. For the typical temperature of molecular clouds Tinit 10-20K and for the typicaldust opacity ' 0:01cm2  g;1, Figure C.3 shows MF = MFth . Hence we choose equation(C.25), which reads MF ' 10;2M. Now we consider how MF varies with the decreasing
dust opacity. As mentioned in xC, lowering under a xed Tinit corresponds to lowering
metallicity of the clouds. Equation (C.25) tells us that MF can be signicantly larger
than 10;2M inversely proportionally to 0 for low-metal clouds. Figure C.4 delineates
MF as a function of under xed initial temperatures, Tinit =10K, 20K, and 30K. For
instance, a cloud with Tinit = 10K and = 10;4cm2  g;1 (a hundredth as small as the
typical metallicity) provides MF = 3:7M. This fact might be responsible for the cuto,
if it exists, at the low-mass end of the initial mass function (IMF).

157

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY

1

0

-2

-1

-1.5

-1

-2
-2
-1
0

-3
1
-4

-5
0

.5

1

1.5

2

2.5

3

Fig. C.3.| The minimum Jeans mass for fragmentation (eqs. #C.25]-#C.27]) is shown as
contour curves. Note that the vertical axis represents not 0 but (Tinit) (see eq. #C.16]).
Contour levels begin with log10 MF=M = ;2 in the upper-left corner, and increase
as ;1:5 ;1 ;0:5   . The dashed line delineates MFth = MF 1 and, equivalently,MFdif = MF 1 according to Figure C.1. The minimum Jeans mass is given by MFthbelow the dashed line and by MFdif above.
158

C.2. THE MINIMUM JEANS MASS FOR FRAGMENTATION

1

0

-1

-2

-5

-4

-3

-2

-1

0

1

Fig. C.4.| The minimum Jeans masses for fragmentation under given temperatures.
Solid, dotted, and dashed curves correspond to Tinit =10K, 20K, and 30K, respectively.

159

APPENDIX C. CONDITIONS FOR VIOLATION OF ISOTHERMALITY

C.3 SUMMARY
Our ndings in this chapter are summarized as follows.
1. We derived three criteria for the central density crit when isothermality is violated ingravitationally collapsing clouds. These criteria are supported by numerical resultsin our previous work. We classied crit in the Tinit- plane, and found that =; ordif MD determines crit in plausible ranges of parameters for actual molecular clouds.
Another criterion, crit =  1, actually has no importance.
2. This result indicates that the condition of \  1" never violates isothermality.
Instead, non-isothermal evolutions begin either earlier or later, depending on theinitial temperature and opacity, than the optical depth reaches unity. The conventional idea that the isothermal approximation is valid during the central density isless than 10;13 g  cm;3 is not supported physically.
3. The minimum Jeans mass for fragmentation, MF, is characterized by crit , andhence MF is a function of both Tinit and as well as crit. The typical value of MF
is 10;2M in accordance with the commonly believed value, but MF can increase
substantially depending on Tinit and . In particular, molecular clouds with lower
metallicity yield larger MF inversely proportionally to metallicity. For Tinit = 10K,
a cloud with = 10;4cm2  g;1 , which is smaller by a hundred than the typicalvalue for nearby molecular clouds, yields MF = 3:7M.
The critical densities derived here would be helpful to numerical calculations in which
either isothermal or a more harder EOS (e.g.,  = 5=3) is chosen correspondingly towhether  < crit or  > crit .
In the present investigation, we suppose that the opacity per unit mass is independentof the density (i.e., assuming LTE) and is homogeneous throughout the cloud. This
assumption is justied when the cloud density is so high that interstellar dust is thedominant coolant, but it is inapplicable for lower densities where the cooling by molecularlines mainly contributes to the opacity. For the molecular line cooling, LTE is not satisedin general and the assumption that the opacity per unit mass is homogeneous is violated
because the line feature depends on the velocity structure of clouds. Our future works will
be intended to generalize the present results to be applicable to the molecular or atomic
line opacity.

160

Bibliography
Adams, F. C., Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788
Adams, F. C. & Shu, F. H. 1986, ApJ, 308, 836
Andre, P., Ward-Thompson, D, & Barsony, M. 1993, ApJ, 406, 122
Alexander, D. R. & Ferguson, J. W. 1994, ApJ, 437, 879
Appenzeller, I. & Tscharnuter, W. 1974, A&A, 30, 423
Appenzeller, I. & Tscharnuter, W. 1975, A&A, 40, 397
Auer, L. H. 1971, J.Quant.Spectrosc.Radiat.Transfer, 11, 573
Balluch, M. 1988, A&A, 200, 58
Balluch, M. 1991a, A&A, 243, 168
Balluch, M. 1991b, A&A, 243, 187
Balluch, M. 1991c, A&A, 243, 205
Barsony, M. 1994, in Clouds, Cores, and Low Mass Stars, ed. D. P. Clements & R.
Barvainis (Provo:ASP), 197
Beckwith, S. V. W., Sargent, A. I., Chini, R. S, & G*usten, R. 1990, AJ, 99, 924
Black, D. C. & Bodenheimer, P. 1975, ApJ, 199, 619
Bodenheimer, P. & Sweigart, A. 1968, ApJ, 152, 515
Bonner, W. B. 1956, MNRAS, 116, 351
Boss, A. P. 1984, ApJ, 277, 768
Boss, A. P. 1988, ApJ, 331, 370
161

BIBLIOGRAPHY
Boss, A. P. & Myhill, E. A. 1992, ApJS, 83, 311
Boss, A. P. & Yorke, H. W. 1995, ApJ, 439, L55
Buchler, J. R. 1979, J.Quant.Spectrosc.Radiat.Transfer, 22, 293
Burkert, A., Bate, M. R., & Bodenheimer, P., 1997, MNRAS, 289, 497
Calvet, N., Hartmann, L., Kenyon, S. J., & Whitney, B. A. 1994, ApJ, 434, 330
Castor, J. I. 1970, MNRAS, 149, 111
Castor, J. I. 1972, ApJ, 178, 779
Chandrasekhar, S. 1967, An Introduction to the Study of Stellar Structure, (Dover ed.)Chen, H., Myers, P. C., Ladd, E. F., Wood, D. O. S. 1995, ApJ, 445, 377
Chick, K. M., Pollack, J. B., & Cassen, P. 1996, ApJ, 461, 956
Choi, M., Evans II, N. J., Gregersen, E. M., & Wang, Y. 1995, ApJ, 448, 742
Close, L. M., Roddier, F., Northcott, M. J., Roddier, C., & Graves, J. E. 1997, 478, 766
Cohen, M., Harvey, P. M., Schwartz, R. D., & Wilking, B. A. 1984, ApJ, 278, 671
Courant, R., Friedrichs, K. O., & Levy, H. 1928, Math.Ann., 100, 32
Dor, E. A. & Drury, L. O'C. 1987, J.Comp.Phys., 69, 175
Ebert, R. 1955, Zs. Ap., 217
Efstathiou, A. & Rowan-Robinson, M. 1990, MNRAS, 245, 275
Efstathiou, A. & Rowan-Robinson, M. 1991, MNRAS, 252, 528
Flower, D. R. & Watt, G. D. 1984, MNRAS, 209, 25
Ensman, L. 1994, ApJ, 424, 275
Goldsmith, P. F. & Langer, W.D. 1978, ApJ, 222, 881
Goldreich, P. & Kwan, J. 1974, ApJ., 189, 441
Green, S. 1975, ApJ, 201, 366
162

BIBLIOGRAPHY
Green, S. & Chapman, S. 1978, ApJS, 37, 169
Gregersen, E. M., Evans II, N. J., Zhou, S., & Choi, M. 1997, ApJ, 484, 256
Hayashi, C. 1970, Proceedings of the 15th Liege International Astrophysical Symposium,
p.127
Hayashi, C. & Nakano, T. 1965, Prog.Theor.Phys., 34, 754
Hayashi, M., Ohashi, N., & Miyama, S. M. 1993, ApJ, 418, L71
Hirano, N., Kameya, O, Mikami, H., Saito, S., Umemoto, T., & Yamamoto, S. 1997,
ApJ, 478, 631
Hunter, C. 1977, ApJ, 218, 834
Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
Inutsuka, S. & Miyama, S. M. 1992, ApJ, 388, 392 (IM92)
Inutsuka, S. & Miyama, S. M. 1997, ApJ, 480, 681 (IM97)
Kanjilal, T. & Basu, B. 1992, Ap&SS, 193, 17
Kenyon, S. J., Hartmann, L., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869
Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676
Kikuchi, N., Nakamoto, T., & Ogochi, K. 1999, in preparation
Lada, C. J. 1987, in Starforming Regions: IAU Symp.115, eds. Peimbert & Jugaku(Reidel), 1
Larson, R. B. 1969, MNRAS, 145, 271 (L69)
Leung, C. M. 1975, ApJ, 199, 340
Leung, C. M. & Liszt, H. S. 1976, ApJ, 208, 732
Leung, C. M. & Brown, R. L. 1977, ApJ, 214, L73
Low, C. & Lynden-Bell, D. 1976, MNRAS, 176, 367
Mardones, D., Myers, P. C., Tafalla, M., Wilner, D. J., Bachiller, R., & Garay, G. 1997,
ApJ, 489, 719
163

BIBLIOGRAPHY
Masunaga, H. & Inutsuka, S. 1999, ApJ, Vol. 510, the January 10 issue, in press
Masunaga, H. , Miyama, S. M. , & Inutsuka, S. 1998, ApJ, 495, 346 (MMI)
Mihalas, D. & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (NewYork: Oxford Univ. Press) (MM)Miyama, S. M., Narita, S., & Hayashi, C. 1987a, Prog.Theor.Phys., 78, 1051
Miyama, S. M., Narita, S., & Hayashi, C. 1987b, Prog.Theor.Phys., 78, 1273
Mizuno, A., Onishi, T., Hayashi, M., Ohashi, N., Sunada, K., Hasegawa, T., & Fukui,
Y. 1994, Nature, 368, 719
Momose, M, Ohashi, N., Kawabe, R., Hayashi, M., & Nakano, T. 1996, ApJ, 470, 1001
Momose, M, Ohashi, N., Kawabe, R., Nakano, T., & Hayashi, M. 1998, ApJ, 504, 314
Monteiro, T. S. 1985, MNRAS, 214, 419
Myers, P. C. & Ladd, E. F. 1993, ApJ, 413, L47
Nakano, T. 1979, PASJ, 31, 697
Nakano, T. 1984, Fund.Cosm.Phys., 9, 139
Nakano, T. 1998, ApJ, 494, 587
Nakano, T., Hasegawa, T., & Norman, C. 1995, ApJ, 450, 183
Nakano, T., Ohyama, N, & Hayashi, C. 1968, Prog,Theor,Phys., 39, 1448
Narita, S., Nakano, T., & Hayashi, C. 1970, Prog.Theor.Phys., 43, 942
Natta, A. 1993, ApJ, 412, 761
Ohashi, N., Hayashi, M., Ho, P. T. P., & Momose, M. 1997, ApJ, 475, 211
Onishi, T., Mizuno, A., Kawamura, A., Ogawa, H., & Fukui, Y. 1998, ApJ, 502, 296
Osterloh, M. & Beckwith, S. V. W. 1995, ApJ, 439, 288
Penston, M. V. 1969, MNRAS, 144, 425
Pollack, J. B., Hollenbach, D., Beckwith, S., Simonelli, D. P., Roush, T., & Fong, W.
1994, ApJ, 421, 615
164

BIBLIOGRAPHY
Preibisch, Th., Ossenkopf, V., Yorke, H. W., & Henning, Th. 1993, A&A, 279, 577
Preibisch, T., Sonnhalter, C., & Yorke, H. W. 1995, A&A, 299, 144
Raich, J. C. & Good, Jr., R. H. 1964, ApJ, 139, 1004
Rees, M. F. 1976, MNRAS, 176, 483
Rybicki, G. B. & Hummer, D. G. 1991, A&A, 245, 171
Saito, M., Kawabe, R., Kitamura, Y., & Sunada, K. 1996, ApJ, 473, 464
Saumon, D., Chablier, G., & Van Horn, H. M. 1995, ApJS, 99, 713
Shu, F. H. 1977, ApJ, 214, 488
Shu, F. H. 1995, in Molecular Clouds and Star Formation, eds. C. Yuan & J. H. You(Singapore:World Scientic), 97
Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23
Silk, J. 1977, ApJ, 214, 152
Snell, R. L., & Loren, R. B. 1977, ApJ, 211, 122
Sobolev, V. V. 1960, Moving Envelopes of Stars, (Harvard Univ. Press)
Spiegel, E. A. 1957, ApJ, 126, 202
Stahler, S. W.,Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637 (SST)
Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 (SMN)
Tohline, J. E. 1980, ApJ, 239, 417
Tomisaka, K. 1998, ApJ, 502, L163
Truelove, J. K., Klein, R. I., Mckee, C. F., Holliman II, J. H., Howell, L. H. & Greenough,J. A. 1997, ApJ, 489, L179
Truelove, J. K., Klein, R. I., Mckee, C. F., Holliman II, J. H., Howell, L. H., Greenough,J. A., & Woods, D. T. 1998, ApJ, 495, 821
Tscharnuter, W. M. 1987, Models of Star Formation in Lecture Notes in Physics : Physical Processes in Comets, Stars and Active Galaxies, eds. Meyer-Hofmeister, E., Thomas,
H. C., & Hillebrandt, W., Springer Verlag Berlin, p. 96
165

BIBLIOGRAPHY
Tscharnuter, W. M. & Winkler, K.-H. 1979, Comp.Phys.Comm., 18, 171 (TW)
Unno, W. & Spiegel, E. A. 1966, PASJ,18, 85
von Neumann, J. & Richtmyer, R. D. 1950, Appl.Phys., 21, 232
Ward-Thompson, D., Scott, P. F., Hills, R. E., & Andre, P. 1994, MNRAS, 268, 276
Ward-Thompson, D., Chini, R., Kr*ugel, E., Ande, P., & Bontemps, S. 1995, MNRAS,
274, 1219
Westbrook, C. K. & Tarter, C. B. 1975, ApJ, 200, 48
Winkler, K.-H. A. & Newman, M. J. 1980a, ApJ, 236, 201 (WN)
Winkler, K.-H. A. & Newman, M. J. 1980b, ApJ, 238, 311 (WN)
Yorke, H. W. 1979, A&A, 80, 308
Yorke, H. W. 1980, A&A, 85, 215
Yorke, H. W., Bodenheimer, P., & Laughlin, G. 1993, ApJ, 411, 274
Yorke, H. W., Bodenheimer, P., & Laughlin, G. 1995, ApJ, 443, 199
Yorke, H. W. 1998, in private communication
Zhou, S. 1992, ApJ, 394, 204166