THE ASTROPHYSICAL JOURNAL, 495 : 346È369, 1998 March 1( 1998. The American Astronomical Society. All rights reserved. Printed in U.S.A.

A RADIATION HYDRODYNAMIC MODEL FOR PROTOSTELLAR COLLAPSE. I. THE FIRST COLLAPSEHIROHIKO MASUNAGA,1 SHOKEN M. MIYAMA, AND SHU-ICHIRO INUTSUKA
Division of Theoretical Astrophysics, National Astronomical Observatory, Mitaka, Tokyo 181, Japan ; masunaga=milano.mtk.nao.ac.jp,
miyama=yso.mtk.nao.ac.jp, inutsuka=yso.mtk.nao.ac.jp
Received 1997 July 14 ; accepted 1997 October 10

ABSTRACT
Dynamical collapse of a molecular cloud core and the formation of a star are investigated by performing radiation hydrodynamic calculations in spherical symmetry. The angle-dependent andfrequency-dependent radiative transfer equation is solved without any di†usion approximations, and theevolution of the spectral energy distribution (SED) is examined.
In the present paper, as the Ðrst step in a series of our work, evolutions before hydrogen molecules
begin to dissociate (the so-called Ðrst collapse) are examined for di†erent masses and initial temperaturesof the parent cloud cores and for di†erent opacities. Numerical results for a typical case [T \ 10 K
init D0.05
and i (10 K) D 0.01 cm2 g~1] show that the radius and mass of the Ðrst core are D5 AU andP
M , respectively. These values are independent both of the mass of the parent cloud core and of the_ density proÐle. The analytical expressions for the radius, mass, and accretion luminosity of theinitialÐrst core are also obtained. The SED contains only cold components of a few times 10 K throughout
the Ðrst collapse phase, because the opaque envelope veils the Ðrst core from observers. We suggest thatthe molecular cloud cores with luminosities higher than D0.1 L should contain young protostars deepin the center, even if they show no evidence for the existence of _central stars in near-infrared and opticalobservations.
Subject headings : hydrodynamics È ISM : clouds È methods : numerical È radiative transfer Èstars : formation
1. INTRODUCTION

shock at the surface of the star or the disk. The thermalradiation 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 sufficiently transparent to the radiation from the central objects. Bipolar outÑows are supposed to be associated with this phase of evolution.
When the accreting envelope becomes transparent to
optical light, the star and protoplanetary disk appear as a T
Tauri star. The dynamical evolution of a protostar is followed by the quasi-static contraction of the central star andthe hypothetical viscous evolution of the protoplanetary
disk.
Theoretical investigations of the above scenario involve
accurate solutions of two-dimensional or three-dimensionalradiation magnetohydrodynamic equations, which inevitably require numerical calculations. At present, however, thecomputational facilities do not have sufficient power for thistask, and therefore a number of simpliÐcations are supposed in such investigations. In the past studies that dealt
with the radiative energy transport, the di†usion (orEddington) approximations and/or the gray approximation
are often introduced for simplicity. For instance, both wereassumed by Larson (1969) and Appenzeller & Tscharnuter(1974). In recent works, Boss & Myhill (1992) introduced
the gray three-dimensional Eddington approximation (Boss
1984) for incorporating the radiative energy transport into
their radiation hydrodynamic code. The di†usion or
Eddington approximation, however, is justiÐed only in optically thick regions such as stellar interiors. Considering thispoint, Narita, Nakano, & Hayashi (1970) adopted the di†usion scheme only for the inner opaque region and solved theangle-dependent radiative transfer in the tenuous andextended envelope. Besides the classical di†usion approximation, the Ñux-limited di†usion (FLD) approximation is
introduced by some recent works (Yorke, Bodenheimer, &Protostellar evolution has been investigated actively
during the past two decades. A standard evolutionary scenario for the formation of isolated low-mass stars can bebrieÑy summarized as follows (e.g., Larson 1969).
The scenario begins with the gravitational collapse of amolecular cloud core. The collapse is initiated by the ambipolar di†usion of magnetic Ðelds that supported the cloudcore against gravitational collapse (e.g., Shu, Adams, &Lizano 1987 ; Shu 1995).
The collapsing cloud core is initially optically thin to thethermal emission from dust grains, and the compressional
heating rate by the collapse is much smaller than thecooling rate by the thermal radiation. The cloud therefore
collapses isothermally.
The isothermal condition is broken when the opticaldepth of the entire cloud core becomes of the order of unity,or when the compressional heating rate overwhelms theradiative cooling rate. The collapsing cloud thus forms an
adiabatic core at its center, which is the so-called ““ Ðrstcore.ÏÏ The evolution at this phase is called the ““ Ðrst collapse.ÏÏ
When the central temperature reaches about 2000 K,
hydrogen molecules begin to dissociate into atoms. Because
released gravitational energy is consumed by the dissociation, gas pressure cannot increase rapidly enough to
support the Ðrst core against its self-gravity, and then the““ second collapse ÏÏ begins. After the dissociation is completed, the ““ second core ÏÏ is born.
The surrounding envelope continues to infall through the““ main accretion phase.ÏÏ The gas cloud with signiÐcant
angular momentum forms a protoplanetary disk aroundthe central star. The infalling gas is thermalized by the
1 Also at the Department of Astronomy, University of Tokyo.

346

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.
Laughlin 1993, 1995). The FLD scheme saves the radiativedi†usion equation from yielding unphysical solutions in
optically thin regions, but still it is not assured to give thecorrect solutions.
The exact treatment of the radiation Ðelds in radiationhydrodynamics (RHD) might require some simpliÐcations
on the dynamics in order to reduce the computationale†ort. Indeed, the past works that dealt with the radiationÐelds can be classiÐed into two kinds of categories, i.e., (1)two-dimensional radiative transfer without hydrodynamics,or (2) one-dimensional RHD. The former includes Efstathiou & Rowan-Robinson (1990, 1991) and Chick, Pollack,& Cassen (1996). They calculated the temperature distributions that are consistent with the radiation Ðelds in axisymmetry. In their calculations they Ðxed the density
distributions as snapshots of the dynamical evolution,
instead of solving the hydrodynamic equations. Moreover,the density distributions assumed by them were solutions ofisothermal collapse, and therefore the temperature anddensity distributions are not self-consistent in their results.
The latter gives the self-consistent solutions with respect
to all the radiation hydrodynamic quantities. A couple ofcomprehensive works that belong to this category were
published in the same year : Stahler, Shu, & Taam (1980)
and Winkler & Newman (1980a, 1980b ; hereafter WNa and
WNb). Stahler et al. (1980) assumed steady state inÑow inthe envelope and divided the problem into several parts
that were solved individually, and the solutions in each partwere connected at the interfaces. In contrast, WNa andWNb pursued the dynamical evolution by solving the radiation hydrodynamic equations straightforwardly using thevariable Eddington factor method, which is employed also
in this paper. Both of them yielded consistent results witheach other, but cases other than a 1 M cloud with a single_ for future investiinitial condition were left to be examinedgators. Furthermore, both of them introduced the grayapproximation, in which they replaced the energy mean andÑux mean opacities with the Rosseland mean opacity
instead of solving the frequency-dependent radiative transfer equation. The gray approximation is valid where thetemperature varies spatially in a scale much larger than amean free path of thermal photons. This is, however, not thecase at the latter stage of the protostellar evolution (i.e., in
the main accretion phase). Indeed, Preibisch, Sonnhalter, &Yorke (1995) showed that the gray approximation cancause considerable errors in estimating the temperature distribution and the radiative acceleration where an opticallythin envelope is heated up by the central hot star.
Frequency-dependent calculations for radiation, without
the gray approximation, were performed by Yorke (1979,1980) for massive star formation. He solved the hydrodynamic equations simultaneously to see the structural evolution of a circumstellar envelope ; he did not resolve thestructure of the central star, but put a central radiative
source there instead.
Our aim in the present investigation is to construct theoretical models that are based on the exact treatment ofradiation Ðelds by solving the angle-dependent andfrequency-dependent transfer equation. Using the results ofcalculations, we discuss (1) conditions that change the isothermal phase into the adiabatic, (2) the radius and accretion luminosity of the Ðrst core, and (3) possible
interpretations of observations. Masses and temperatures of
initial molecular cloud cores are given as model parameters,

347

in contrast to Stahler et al. (1980) and WNa and WNb. The
dust opacity is also a model parameter because its dependence on frequency is uncertain. This is because of the lackof knowledge of the chemical composition, size distribution,and geometrical conÐguration of dust grains. We adopt
spherical symmetry, which inevitably requires that thee†ects of rotation and magnetic Ðelds should be neglected.
Both the e†ects, of course, actually have important roles inprotostellar evolution. Nevertheless, the spherically symmetric calculation still has signiÐcance for the following two
reasons : (1) the exact treatment of radiation Ðeld in general
geometry is quite difficult, and (2) moving coordinates areavailable only for one-dimensional calculations. A molecular cloud core of o D 10~19 g cm~3 and R D 0.1 pc contracts to a star of o6 D 1 g cm~3 and R D 1 R , while the_density increases by 20 orders of magnitude and the radiusdecreases by 6 orders. Any multidimensional numerical calculations to date that include adapted mesh reÐnement
technique and smoothed particle hydrodynamics have notsucceeded in covering such a wide dynamic range in a selfconsistent way.
In the present paper, as the Ðrst step toward our Ðnalaim, we discuss the properties of the Ðrst collapse, whichcorresponds to the evolution prior to the dissociation ofhydrogen molecules. After discussing possible contributions
to observational studies in the following section, the formulation of the basic equations is presented in ° 3. The initialand boundary conditions and other comments on calculations are given in ° 4. The result for a typical case isoverviewed in ° 5. In ° 6 we discuss the dependence of theradius of the Ðrst core R on model parameters such as themass, initial temperature,fc and opacity. We derive analyticalexpressions for R and the accretion luminosity L
in ° 7.
acc
We also discuss fcpossible interpretations of observational
data with our results. In ° 8 we summarize the results. The
computational scheme is described in Appendices A, B, andC with some useful test problems. The subsequent evolutions triggered by the second collapse will be discussed in
our next paper. The formation of intermediate-mass stars(M [ 3 M ) is also discussed in the subsequent papers.
_
2. POSSIBLE CONTRIBUTIONS TO OBSERVATIONS

Our frequency-dependent RHD calculations provide theemergent spectral energy distribution (SED), which has, inprinciple, information on the internal structure of theobject. The well-known classiÐcation of the young stellar
objects (YSOs)Èclass I, II, and IIIÈis based on the shapesof SEDs in the near-infrared and is supposed to represent
the evolutionary sequence of YSOs (e.g., see Lada 1987).
The SED of class I, which is the youngest among the three,
can be explained in terms of absorption and reemission (andscattering) of thermal radiation from the central objects (i.e.,the star and protoplanetary disk) by the circumstellar
envelope (Kenyon, Calvet, & Hartmann 1993). The deÐnition of classes IÈIII requires that the objects should be seenin near-infrared light, which excludes much younger objects
from the classiÐcation.
On the other hand, recent observations have discoveredstar-forming cloud cores without any optical or infrared
sources (i.e., without IRAS sources). Some of them haveevidence of protostellar activity such as bipolar outÑows
(e.g., VLA 1623 observed by Andre, Ward-Thompson, &Barsony 1993). Andre et al. (1993) considered VLA 1623 asan object younger than class I and classiÐed it into ““ class

348

MASUNAGA, MIYAMA, & INUTSUKA

0.ÏÏ So-called starless cores and class 0 objects have SEDs
that can be Ðtted with graybody spectra of T approximately
a few times 10 K. These cold SEDs seem to indicate thatthey are substantially younger than class I protostars.
Moreover, some of the starless cores are interpreted bysome authors as the cloud cores contracting quasi-staticallyby the ambipolar di†usion of magnetic Ðelds (Mizuno et al.
1994 ; Ward-Thompson et al. 1994). However, a protostar
surrounded by a dense envelope might look as if it werestarless even after the central stellar core was born.
The recent progress in observations demands theoreticalworks that model the early stage of protostellar evolution.
For example, Boss & Yorke (1995) investigated early protostellar evolution and discussed the possibility that very
young protostars are detected in the present and future
observational projects. Our investigation, starting with thepresent paper, tries to model the SEDs of such very young
protostellar candidates. Our suggestions about the observations will be presented brieÑy in ° 7.4.
3. THE BASIC EQUATIONS OF RADIATIONHYDRODYNAMICS

The basic equations of radiation hydrodynamics are
introduced in this section. Hydrodynamic equations that
account for the interactions with radiation are as follows.
Do] o$ Æ ¿ \ 0 ,
Dt

(1)

D¿
s
\ [o$' [ $p ] F0 F ,
Dt
c 0

(2)

De
] p$ Æ ¿ \ ov ] (cs E [ 4ns B) ,
CR
E0 0
P0
Dt

(3)

+2' \ 4nGo ,(4)(c [ 1)oe \ p 4 p(o, T ) ,

(5)

o
o

where o, ¿, e, and T are the Ñuid mass density, the velocity,
the internal energy density per unit mass, and the temperature, respectively, and s, v , and c are the absorption
CR rate per unit mass, andcoefficient, the cosmic ray heating
the ratio of speciÐc heats, respectively. The Planck function
is denoted by B 4 pT 4/n. E and F , where subscript zero
0 frequency-integrated
denotes comoving variables,0 are theradiation energy density and radiation Ñux, which aredeÐned as

P Q
P Q

E \
0
F \
0

=

0
=

dl

0

dX I (l , n ) ,0 0 0 0dX I (l , n )n ,dl0
0 0 0 0 00where I (l , n ) denotes the speciÐc intensity of radiation at
0 0l along0frequency
the direction vector n in the comoving
0
0 coefficients areframe. The frequency-averaged absorptiondeÐned as follows.
/= s (l )F (l ) dl
\ 0 0 0 0 0 0,F0/= s (l )E (l ) dl
s \ 0 0 0 0 0 0,
E0
E
0
s

F0

Vol. 495/= s (l )B(l ) dl
0 0.
\ 0 0 0P0Bs

Equations (1)È(5) correspond to the equation of continuity, the equation of motion, the Ðrst law of thermodynamics for material, PoissonÏs equation, and the equation ofstate, respectively. They are the ordinary hydrodynamic
equations except for the last terms on the right-hand side ofequations (2) and (3), which describe interactions between
material and radiation. To account for the radiation Ðeldexactly, the transfer equation of radiation (eq. [6]) must be
solved.
1 LI(l, r, n)] n Æ $I(l, r, n)cLt\[ s(l, r)[I(l, r, n) [ S(l, r)] ,(6)

where S(l, r) is the source function. Note that equation (6) iswritten in the lab frame.
Because the speciÐc intensity of radiation depends onfrequency and on direction angle as well as on location, it
requires considerable computational e†ort to solve
equation (6) directly. Moreover, equation (6) must be transformed into the comoving frame to couple with the hydrodynamic equations, which would modify the equation into aquite complicated form. In contrast, it reduces the computational e†ort to an acceptable level to solve the moment
equations of radiation in the comoving frame (eqs. [7] and[8]) 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 theirstudies were succeeded by Mihalas & Mihalas (1984, hereafter MM) in a comprehensive work 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 that are O(v/c) in all physicalregimes. Therefore the adopted equations can account for
all the dynamical e†ects of radiation in a nonrelativisticÑuid (see SMN).

A B

o

A B

o D F
s F
0 ] $ Æ P \ [ F0 0 ,0c2 Dt oc(7)D E
0 ] $¿ : P ] $ Æ F \ [(cs E [ 4ns B) ,
0
0
E0 0
P0
Dt o
(8)

where P is the radiation stress tensor :0=dldXI (l , n )n n .
P \0
0 0 0 0 000Equations (7) and (8) correspond to the momentum equation and the energy equation of radiation, respectively.
In our scheme, equation (3) has been replaced with thefollowing one.

P Q

o

A

B

E
D
e ] 0 ] $¿ : P ] $ Æ F ] p$ Æ ¿ \ ov , (9)
0
0
CR
o
Dt

which is given by taking the sum of equations (3) and (8).
Equation (9) ensures that the numerical scheme will be
more accurate than equation (3) when the Ñuid is near radiative equilibrium, because the interaction terms between
material and radiation have canceled (SMN).

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

Equations (7) and (8) are not closed in themselves since
they have been integrated over direction angle. Therefore, aclosure relation that compensates the lack of knowledge forthe angular distributions of the speciÐc intensity must be
given. As a closure relation, we introduce the tensor variable
Eddington factor f, which is deÐned by
P\fE .

(10)

The radiation stress tensor P is eliminated using the tensor
variable Eddington factor, which is a measure of the degreeof anisotropy of the radiation Ðeld, to close the system of
the RHD equations.
When we calculate f, the radiation Ðeld at each time stepis assumed as a static snapshot under the instantaneousdensity and temperature distributions. In other words, wesolve equation (6), neglecting the term including the time
derivative. This assumption is justiÐed only when the timescale of radiative di†usion or propagation is sufficiently
short compared to the evolutionary timescale of interest.
Fortunately, such a condition is satisÐed in many astrophysical problems, including protostellar collapse. Even if
the condition is not satisÐed, the static transfer equation isstill applicable for most cases in astrophysics, because thesolutions 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 ofO(v/c), which would appear when transforming equation (6)into the comoving frame, have been dropped because therelativistic e†ects such as aberration or Doppler shift hardlya†ect the Eddington factors for the problems of our presentinterest. Thus, the transfer equation is reduced to a fairly
simple form and then no longer bothers us for the unacceptable computational e†ort.
Once I(l, r, n) is obtained, the monochromatic tensor
variable Eddington factor is calculated using the relation :f (l, r) \P(l, r).
E(l, r)

The frequency-integrated Eddington factor is deÐned by/= f (l)E(l)dlf\ 0.
EThe computational method, which is called the variable
Eddington factor method, was established as a techniquefor astrophysical radiation hydrodynamics Ðrst by

349

Tscharnuter & Winkler (1979). There is no essential di†erence between the computational scheme of Tscharnuter &Winkler (1979) and ours, except that they supposed the grayapproximation while we solve the frequency-dependent
radiative transfer equation. This method can be divided into
two parts. In the Ðrst part, the transfer equation (6) with
L/Lt 4 0 is solved under the density and temperature distributions assumed to be known. Thus this part is called thetransfer equation solver. The tensor variable Eddington
factor at each spatial point is obtained from the solutions
I(l, r, n). In the second part, the full set of the RHD equations (1)È(4), (7), and (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 for the Ðrst part, whichis employed again. The two parts are iterated in this mannerfor consistency. In Appendix A we will show the details ofeach part for the spherical symmetrical case. Test calculations are presented in Appendix B.
4. CALCULATIONS

4.1. Initial Conditions
For the total mass of the cloud core M , we choseM \ 1 M , 0.1 M , 0.3 M , and 3 M cc
. The initialcc__ cloud coredensityand _the radius_ are adjusted
so that theshould be slightly more massive than the Jeans mass foreach case. As shown in Table 1, the initial density is o D
initthe10~20 to 10~17 g cm~3 (i.e., n D 104È107 cm~3) and2
H
radius is R D 1000È30,000 AU, which is compatible withthe typical ccscales of observed molecular cloud cores. The
initial density distribution is assumed to be homogeneous
except for a case in which the cloud core is initially inhydrostatic equilibrium (see ° 6.2). In every model, thevelocity of each Ñuid element is set to zero in the initialstate.
The initial temperature is also an important model
parameter, T \ 10, 5, and 30 K. Di†erent initial teminit
peratures accountfor possible variations of the environments in molecular clouds. Dependence of protostellar
evolution on initial temperature will be discussed in ° 6.3.
4.2. Boundary Conditions
The inner and outer boundaries are given by the centerand the edge of the cloud core, respectively. The outerboundary is Ðxed through evolution at R \ R shown in
cc

TABLE 1
MODEL PARAMETERSModel
M1a . . . . . .
M2a . . . . . .
M3a . . . . . .
M4a . . . . . .
C1a . . . . . .
H1a . . . . . .
M1b . . . . . .
M1c . . . . . .

M a(Mcc )
_
1.0
0.1
0.3
3.0
1.0
1.0
1.0
1.0

T
init
(K)

i (T )
P g~1)
init
(cm2

o
init
(g cm~3)

R b
cc
(AU)

R c
fc
(AU)

o d
ad
(g cm~3)

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

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

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

2 ] 10~13
2 ] 10~13
2 ] 10~13
2 ] 10~13
1 ] 10~13
6 ] 10~13
5 ] 10~14
5 ] 10~13

a 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 o \ 10~9 g cm~3.
c
f Evaluatedfrom eq. (28), with 10 for a.

Le
(L )
_
0.061
0.014
0.031
0.12
0.018
0.24
0.14
0.050

L acc f(Lfc )
_
0.04
0.04
0.04
0.04
0.01
0.3
0.1
0.02

350

MASUNAGA, MIYAMA, & INUTSUKA

Table 1. Boundary conditions for radiation involve a brief
discussion as follows.
Boundary conditions for the transfer equations (A1) and(A2) (see Appendix A) are given by h on the p-axis shown inlFigure 13 and at the cloud edge. On the p-axis h \ 0lbecause of symmetry, and h at the cloud edge is determinedlby given incident radiation I~(r \ R ). In the present papercc
two di†erent boundary conditions for radiation thatprovide a couple of extreme cases are considered.
Case 1.ÈThere is no incident radiation but there is a heatsource v to compensate for the radiative energy loss in theCRinitial state. The heat source represents the heating bycosmic rays and photons from surrounding stars, which
support the temperature of the molecular cloud against
cooling by the thermal emission of dust grains (Goldsmith& Langer 1978).
Case 2.ÈThe cloud core is surrounded by isotropic
blackbody radiation of T \ T .
init the cloud core naturally at
Although case 2 seems to keepa constant initial temperature, ambient blackbody radiationof T
is too idealized a condition for actual molecular
initcores. Furthermore, case 1 is more favorable thancloudcase 2 for the purpose 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(eqs. [A12] and [A13]) are given as radiation Ñuxes at the
inner and the outer boundaries. The inner Ñux at the centeris set equal to zero because of symmetry. The outer boundary condition di†ers between the two cases. In case 1, theouter Ñux F
in equation (A12) is eliminated using h,
I`1 by F
where h is deÐned
\ hcE
and is given by theI`1 outer boundary contransfer equation solver. I`1
In case 2, thedition for radiation is given simply by E
4 aT 4 ,I`3@2
init
where a 4 4p/c.
4.3. Opacity
We assume that the opacity is due to the absorption by
dust grains, neglecting the scattering in this paper, and thatthe gas temperature couples closely with the dust temperature. These assumptions are known to be well justiÐed
in molecular cloud cores and in 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 frequency
dependence of the dust opacity in higher frequency than
far-infrared except for the line features, which have been
removed in our models for simplicity. Dust grains are
assumed to consist of silicate and graphite, where H O ice is
excluded since it has no signiÐcant inÑuences 2on the
dynamical evolutions. The opacity for lower frequency ischaracterized with a parameter b, where i \ i (l/l )b,l
0 0
taking the ambiguity of the dust opacity into
account.
Three di†erent models are considered with respect to b.
a) b \ 2.5 for 4.6 ] 1011 \ l \ 3 ] 1012 Hz and b \ 1.5
for l \ 4.6 ] 1011 Hz, taking a recent result into account(Pollack et al. 1994). [This model has a comparable value
for i (10 K) to that of the model of b \ 2.]b) Pb \ 1.5 for l [ 3 ] 1012 Hz.
c) i \ 0.1imodel a.
llVol. 495

In models a and b, i is determined to be connected at0
l \ 3 ] 1012 Hz with the opacity of Adams & Shu (1986),
0that is, i
12 Hz \ 0.3 cm2 g~1. Note that the absorption
3C10
coefficientss are given with the mass opacity of dust i bylls (r) \ i o(r). We shall examine in ° 6.4 how the di†erencellamong models aÈc a†ects the evolution.
We consider the vaporization of dust grains at T 4
vap1200 K. Although the molecular opacities should take theplace of the opacity for higher temperature, we simply
reduce the dust opacity to nearly (but not exactly) zero for
T [ T . In order to assure the convergence of the Eddingvap
ton factor between the transfer equation solver and themoment equations solver, the dust grains are supposed tovanish completely at T \ T ] *T , where *T is taken tovap
be 100 K, and to decrease linearly with increase of temperature for T \ T \ T ] *T .
vap
vap
4.4. Equation of StateWe adopted the equation of state of ideal gas for simplicity :
1
k
T ,
(11)
c [ 1 kmH
where c, k, k, and m are the ratio of speciÐc heats, theBoltzmann constant, Hthe mean molecular weight, and themass of hydrogen atom, respectively. In the present paper, c
is Ðxed at 5/3 except in a case where c equals 7/5 instead,and k is assumed to be constant at 2.46, which corresponds
to the same chemical abundance as in Larson (1969). Morerealistic models should include the e†ects of rotational andvibrational degrees of freedom for H molecules, and of2
dissociation and ionization of the molecules.
The formermodiÐes c from 5/3 and the latter changes the value of k as afunction of the temperature. Moreover, the dissociation andionization energies should be taken into account in theenergy equation. All of these e†ects will be considered in thesubsequent papers.
e\

5. RESULTS

5.1. Overview for a T ypical Case
Results of the typical case (i.e., model M1a in Table 1) areoverviewed in this section. In model M1a, the mass and theinitial temperature are taken to be 1 M and 10 K, respec_ First, we present
tively, with the opacity of model a in ° 4.3.
a brief summary of the results in Figures 1 and 2 with Table
2.
The initial state is indicated by the thick line in eachpanel. At early stages, the cloud core collapses isothermally,
which corresponds to the Ðrst two lines in each panel. The
TABLE 2
ELAPSED TIMES ASSIGNED TO FIGURES 1È3
Time(105 yr)
0.0 . . . . . . . . . . .
2.04830 . . . . . .
2.16677 . . . . . .
2.17139 . . . . . .
2.17301 . . . . . .
2.18122 . . . . . .
2.18725 . . . . . .

Time/t aff
0.0
1.15643
1.223311.225931.22684
1.23147
1.23488

Line Typeb
Thick
...
Dotted...
...
Dashed...

a Thefree-fall
time
is
t 4
ff
[3n/(32Go )]1@2 \ 1.77 ] 105 yr.
init
b For correspondence
to Figs. 1 and 3.

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

351

FIG. 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 2. The initial state is speciÐed by thick solid lines. See text for the meanings of dotted and dashed 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 of change of the speciÐc entropy.

density distribution at large radii develops an r~2 law reminiscent of the isothermal self-similar solution by Larson(1969) and Penston (1969 ; Larson-Penston solution). This
isothermal stage lasts about one free-fall time until thecentral temperature departs signiÐcantly from the initialvalue (dotted lines). The rapid increase of the central temperature to T D 2000 K forms a central adiabatic coreÈ
that is, the Ðrst core. The drop in the density distribution at
a few AU shows the e†ective radius of the Ðrst core, whose
position gradually moves inward through the Ðrst collapse
stage. The supersonically infalling gas in the envelope isdecelerated and thermalized at the surface of the Ðrst core,
where the shock front is observed as a sharp peak in thevelocity distribution. In the radiative Ñux distribution, a few
lines just below the dashed line have small humps at theradius of the Ðrst core, which indicates that the thermalenergy generated by the shock heating is immediately radiated away. The distribution of change rate of the speciÐc
entropy shows the negative peak on the surface of the Ðrstcore, where the entropy is taken out by radiation. In termsof the mass distribution, the inner region has a r3 slope that
corresponds to the homogeneous density structure insidethe Ðrst core, while the outer region is proportional to theradius, which resembles the Larson-Penston solution in theinfalling envelope (i.e., o P r~2). Therefore, the plateau
between these two slopes implies the mass of the Ðrst core,352

MASUNAGA, MIYAMA, & INUTSUKA

Vol. 495

FIG. 1ÈContinued

which reads 0.05 M in the Ðnal stage of this calculation.
_ with the result by Boss & Myhill
This value is consistent(1992). Figure 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, isincreasing up to D0.05 M .
_ the stage when the dust grains
The dashed lines represent
begin to evaporate in the central region. The complicatedbehavior in the radiative Ñux for R \ 1 AU is caused by thesudden decrease of the opacity due to the dust vaporization.
The evolution of the spectral energy distribution (SED) ispresented in Figure 3. The panels except the lower rightpanel suppose observations with di†erent beam sizes of thetelescopes. These panels show that evolutionary changes

would be detected more clearly by observations withsmaller beam sizes. This is because a Ðner beam can concentrate the Ðeld of view on a more restricted region in thecenter of the collapsing cloud core, where the evolution isproceeding 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 increaseof temperature in the Ðrst core. The lower right panel in
Figure 3 gives us an answer to the question of why SEDs
are not aware of the existence of the hot central core. The
Ðrst core of T D 2000 K radiates mostly in the frequency
range near l D 1013 Hz, at which R(q \ 1) [ 100 AU, while
the size of the Ðrst core is only a few AU. An observer

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

353

FIG. 2.ÈEvolutionary sequences for a typical case (model M1a) as functions of mass. See also the legend of Fig. 1. The Ðgure depicts (a) temperature,(b) density, (c) infall velocity, and (d) Planck mean optical depth.

therefore cannot see through down to the surface of the Ðrstcore, but only observes the cold envelope of T D 20 K.
Second, we examine the physical processes that determine the thermal evolution from the isothermal collapse to
the adiabatic contraction. Figure 4 illustrates the thermal
evolution at the center of the collapsing cloud core for thesame model as above (i.e., M1a in Table 1). In Figure 4 theevolution proceeds from left to right as the central density
increases monotonically with time.
Initially the radiative cooling rate, " 4 4ni B, and theth
P the radicosmic ray heating rate, v , are balanced,
whileCR
ative heating rate, ci E, and the gas compressional heating
rate, ! 4 (p/o)$ Æ ¿, Eare considerably smaller as comparedgto them. Here the heating and cooling rates are evaluatedper unit mass. As the cloud core collapses, ! shows anincrease proportional to the square root of thegdensity foro Z 10~17 g cm~3, which is accounted for in equation (13).
c slight disturbance of the compressional heating rate
Theappearing in o [ 10~17 g cm~3 is caused by the homogec initial condition, which is far from theneous and static
similarity solution. The radiative heating rate also increases
with the central density, reÑecting that the dust grains (theemitter of photons) are being concentrated into the central
region.
In principle, the increase of temperature due to the compressional heating can raise the radiative cooling rate

354

MASUNAGA, MIYAMA, & INUTSUKA

Vol. 495

FIG. 3.ÈEvolution of spectral energy distributions (SEDs) for model M1a. Panels except the lower right show SEDs integrated over the di†erent beam
sizes of an imaginary telescope. A beam size corresponds to a concentric circle of the radii indicated in the panel. The lower right panel shows the radius atwhich the optical depth measured from an observer reaches unity as a function of frequency.

because " has a positive dependence on T . For o \th however, the radiative cooling rate is muchc
10~14 g cm~3,larger than the compressional and radiative heating rates,and " is a very steep function of temperature, so that ath increase of temperature is enough to compensate
very slightfor the compressional heating. Thus the gas stays isothermal for o \ 10~14 g cm~3.
When c! gets sufficiently close to " at o D 10~14 ggth thec cloud core
cm~3, the radiative
cooling no longer keeps
isothermal against the compressional heating. Thus, thecentral temperature T departs from T . At o D 10~12 gcinit
cm~3 the radiative heating
rate overtakes" c , when thethoptical depth to the thermal radiation reaches the order ofunity. (Note that the condition ““ the radiative heating rate isapproximately the radiative cooling rate ÏÏ corresponds toq Z 1 only when the radiation Ðeld has no external contributions.) Further evolution occurs without efficient coolingand thus enters the adiabatic phase. The central temperature increases linearly in the log T Èlog o plane withcc behavior
the gradient of D2/3 for o Z 10~12 g cm~3.
This
cproves that the isothermal collapse has been replaced by
the adiabatic contraction because in adiabatic evolutions
d log T /d log o equals c [ 1, with c [ 1 \ 2/3 for c \ 5/3.
The adiabatic evolution decelerates the self-gravitational

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

355

solution). Thus, the gas compressional heating rate for selfgravitational collapse becomes! \ C c2 J4nGo ,
(13)
g
1 s
where c is the isothermal sound speed for the initial tems
perature. Our numerical calculations implied that C \ 1.0
1provides a good approximation for ! . Equating " withgth
! in equation (12) using equation (13), one has the centralgdensity when " \ ! as
th
g
C ~2
o
\ 1.4 ] 10~14 g cm~3 1
"/!
1.0
]

A

A B
BA B

i
2 T
6
P
init
.
1.6 ] 10~2 cm2 g~1
10 K

(14)

We regard this analytic expression as the characteristic
density for the violation of isothermality in the opticallythin case. On the other hand, the condition that the
optical depth of the thermal emission becomes unity (q D1) also deÐnes another characteristic density. For ccthe
purpose of calculating the optical depth of the cloud, we
approximate the density distribution as o \ o (r/R )~n,where n \ 0 inside of r \ R and n \ 2 outside. cR isb theb andboundary between the innerb and the outer regions,approximately accords with the Jeans length corresponding
to the current central density. That is,

FIG. 4.ÈThermal evolution at the center of a collapsing cloud core formodel M1a. (a) The heating and cooling rates per unit mass as functionsof the central density. The evolution proceeds from left to right. (b) The
log T Èlog o plane.

collapse and thus diminishes the value of d ln ! /d ln o togcless than 1/2. The abrupt drops of the radiative heating
andcooling rates at o D 10~9 are due to the dust vaporization.
c that our results are consistent with past
It can be shown
works. In Larson (1969) the ““ initial ÏÏ mass and radius of theÐrst core when the central density reaches 170 K are 1031 g
and 6 ] 1031 g cm~3, respectively, which are roughly equal
to our results at a corresponding evolutionary stage (seeFig. 1). Figure 4 shows good agreement with Figure 3 in
Boss (1988), who performed three-dimensional RHD calculations with the three-dimensional Eddington approximation. Moreover, his result for the lower opacity model isalso quite similar to our lower opacity model (model M1c,
see ° 6.4 and Fig. 9). These imply that our results in spherical symmetry would be helpful for investigations in generalsymmetry.
5.2. Characteristic Densities
The characteristic central density when ! reaches "gcan be evaluated as follows. The explicit expressionsof "thth
and ! are
g
d 1
.
(12)
" \ 4i pT 4 , ! \ [p
th
P
g
dt o

AB

The central density evolves with the dynamical timescale
D(4nGo)~1@2 and so that d/dt can be replaced by
C (4nGo3)1@2(d/do), where C is a dimensionless constant of
1 order of unity (C \11.5 for the Larson-Penston
the
1

2nc
s
R \C
(15)
b
2 J4nGocwhere a dimensionless constant C is 0.75 by comparison
2with numerical results. q is then calculated
as follows.
cc
Rcc
Rb
Rcc
r ~2
io dr Dio dr ]
io
dr
1Dq \
c
c R
cc
0
0
Rb
b
R
\ io R 2 [ b .
c bRccNeglecting R /R in the last term, since it is much less thanb ccunity,

P

P

A

P

A B

B

io R \ 1 .
c b 2
Combining equations (15) and (16), we haveo
\ 4.3 ] 10~13 g cm~3
qF1
C ~1
i
~2 T
~1
2
P
init
.
]
0.75
1.6 ] 10~2 cm2 g~1
10 K

A B A

B A B

(16)

(17)

Equations (14) and (17) are roughly consistent with theresults by Inutsuka & Miyama (1997), who derived them fora cylindrical collapse.
Note that oand odepend sensitively both on the"/! andqF1
initial temperatureon the opacity. Therefore thethermal evolution during the Ðrst collapse under other circumstances can be quite di†erent. We shall discuss the e†ectof di†erent initial temperature and opacity in ° 6.
We also note that the condition of q D 1 does not trigger
the adiabatic phase in cases where the compressional
heating rate is still considerably small compared to the radiative cooling rate when o reaches o , because the actual
qF1E, is still larger thancooling rate by radiation, c4ni B [ ciP! even when 4ni B D ci E in suchE cases. We will cong these situations
P in ° 6.3.2,
E
sider
where another criterion is
proposed instead of equations (14) and (17).

356

MASUNAGA, MIYAMA, & INUTSUKA
6. PROPERTIES OF THE FIRST CORE

In this section we deÐne the mass and radius of the Ðrstcore 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 withobservational data. In the present paper we focus our investigation on estimating the size, mass, and accretion luminosity of the Ðrst core and discuss possible interpretationsof observational evidences, taking SEDs as well into consideration.
6.1. DeÐnitions of the Core RadiusFirst we argue how to deÐne the radius of the Ðrst core,
R . Three di†erent (i.e., dynamical, thermodynamical, andfc
geometrical)deÐnitions are examined as follows.
1. R is deÐned where the gas pressure is balanced withfc
the ram pressure of the infalling envelope.
2. R is deÐned at the position of the peak of ds/dt, wheres is the fc
speciÐc entropy (see Fig. 1h).
3. R is considered as the position where d ln o/d ln r has
fc
the maximum
value (see Fig. 1b).
Figure 5a shows R in each deÐnition as a function of o .
fc the Ðrst core is formed (o Z 10~11c
The lines converge aftercg cm~3) and indicate a few AU for R . We consider that it
fc

Vol. 495

is most desirable to adopt deÐnition 1 [R 4 r(P \fcgasP )], since it shows an intermediate value among the threeram
and is in the best accordance with the analytical prediction
estimated in ° 7.1.
In the remainder of this section, the dependences of R
fc
on the mass, initial temperature, and opacity are argued.
The model parameters are tabulated in Table 1. In eachmodel the calculation is performed until the central temperature reaches D2000 K, because at this temperature theÐrst collapse should actually have been succeeded by thesecond collapse, which is not a subject in this paper.
6.2. Dependence on Cloud Mass
We present here the radii of the Ðrst core for the parent
cloud cores of di†erent masses (M ). Models are denotedccby M2a, M3a, and M4a for M \ 0.1 M , 0.3 M , andcc__
3M , respectively, as well as model M1a for M \ 1 M .
_cc_
The radius of the Ðrst core in each model is deÐned according to deÐnition 1 in ° 6.1 and is plotted in Figure 5.
Interestingly, R as a function of o is not sensitive to
fc not depend on thec initial density and
M . That is, R doesccfc
the size of the parent cloud core. Consequently, the mass ofthe Ð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 the parent cloud

FIG. 5.ÈRadius of the Ðrst core, R . (a) Various deÐnitions for R are shown. The deÐnition R(P \ P ) is represented by the solid line, R(d ln o/d ln
fc \ max) by the dotted line. R(Pfcramr \ max) by the dashed line, and R(ds/dt\ P ) is adopted as R ingasthe following
panels : (b) R for di†erent masses of the
gas
ram
fc
fc
parent cloud core and for an initially inhomogeneous density distribution, (c) for di†erent initial temperatures, and (d) for di†erent
opacities. See text for
denotation of models a, b, and c. The analytical estimation for R assuming 10 K for T and 10~13 g cm~3 for o is shown as well (dot-dashed lines).
fc
init
ad

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

core determines the initial density and its radius. In theisothermal phase the density evolution of the collapsing
cloud core ““ forgets ÏÏ its initial condition and develops into
centrally peaked density proÐle (which resembles theLarson-Penston solution) in the sound crossing time overthe entire cloud core (as noted also by Bodenheimer &Sweigart 1968). Therefore the evolution, regardless of mass,converges to the sequence that is determined uniquely by
the initial temperature. The physical processes that change
the isothermal phase into the adiabatic phase are a†ected
only by the initial temperature and opacity, according toequations (14) and (17), and thus the properties of the Ðrstcore are a†ected by T and i but are independent of theinitPmass 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
density of 10~17 g cm~3 as the initial density distribution,
for which the outer radius truncated at 104 AU. The mass ofthe cloud is 0.87 M . The resultant value of R shows no_ from the other cases, as seen
fc in Figure
signiÐcant di†erences
5b, which supports the forgetful nature of collapsing cloud
cores.
6.3. Dependence on Initial Temperature
The Taurus-Auriga molecular cloud has a typical temperature of 10 K, which motivated us to choose the initial
temperature of T \ 10 K as a typical value. However, theinit stars in the neighborhood of starexistence of massive
forming clouds would raise initial temperatures for the collapse of cloud cores such as in the Orion clouds. In contrast,
the interior of cloud cores may well have a lower temperature than the typical value of the entire cloud because
visual or UV photons from neighboring stars would notreach deep inside cloud cores. Indeed, the value of
9.6 ] 10~2 ergs s~1 g~1 for v , which was found to beCRrequired to sustain the initial cloudcore at T \ 10 K, islarger by 2 orders of magnitude than the typical cosmic ray
heating rate.
Hence, model calculations for T \ 5 K (model C1a)
initto be compared to theand 30 K (model H1a) are carried out
typical case of T \ 10 K. The other conditions are theinitsame as model M1a.
6.3.1. Case for T \ 5 K
init

Equation (14) claims that odecreases quite sensitively"/!
with decreasing initial temperature.
The Planck mean
opacity, moreover, has a lower value for a lower initial
temperature, which makes odiminish with T moreinit
sensitively. For T \ 5 K, we"/!
have
init
o
\ 8.8 ] 10~18 g cm~3 ,"/!
noting that i (5 K) is calculated to be 3.2 ] 10~3 cm2 g~1.
P lowering T by a factor of 2 yields a o
It is found thatthat is lower by 3 orders of init
magnitude. In contrast, o "/!
hasa negative dependence upon the initial temperature qF1(see eq.
[17]). For T \ 5 K one Ðndsinito
\ 2.2 ] 10~11 g cm~3 .
qF1
In this case o
exceeds oby 6 orders of magnitude, in
qF1case of T "/!
contrast to the\ 10 K, where they are comparable to each other. This init
contrast causes a large di†erence inthermal history between the two cases.
Figure 6 illustrates the thermal evolution for T \ 5 K.
init early,
One can see that the isothermal phase ceases rather

357

FIG. 6.ÈThermal evolution at the center of a collapsing cloud core formodel C1a (T \ 5 K). See also the legend of Fig. 4.
init

when o reaches o, and that the central temperaturec gradually "/!
increases
during a long transitional phase before
the adiabatic stage is entered when o exceeds o . In thec o , theqF1
transitional phase, where o\o \compres"/!
c
qF1 the radiative
sional heating rate attempts to dominate
cooling rate, but is pursued immediately by it in response to
the raised central temperature. Consequently ! and " areg graduallythalmost balanced with each other as they increase
with o before the optical depth becomes unity.
Thec isothermal approximation, which many theoreticalworks suppose for collapse calculations, is often considered
to be valid when o \ o
D 10~13 g cm~3. However,under the situations cwhereqF1o>o
as in model C1a,qF1
the transitional phase breaks"/!
the isothermal
condition far
before the adiabatic stage begins, which implies that theisothermal approximation is not so valid even for o \co . This aspect is critical especially for cylindrical colqF1
lapses, where even a slight hardening of the isothermal
equation of state can decelerate the gravitational collapse(Inutsuka & Miyama 1997).
6.3.2. Case for T

init\ 30 K

If the cloud core is surrounded by energy sources such as
OB stars, the temperature of the cloud can be larger fromthe beginning of the collapse. Indeed, cloud cores in theOrion clouds should be exposed to hard radiation Ðelds,
and the average temperature of the Orion clouds is considerably higher than that of dark clouds. Therefore thecalculation with higher initial temperature should be
instructive for the understanding the e†ect of di†erent
environments. Thus in model H1a, we have removed v
CR

358

MASUNAGA, MIYAMA, & INUTSUKA

FIG. 7.ÈThermal evolution at the center of a collapsing cloud core formodel H1a (T \ 30 K). This is a case where v is absent and the incident
radiation Ðeldinitis present. Note the di†erence CRin panel (a) from the othermodels. See also the legend of Fig. 4.

and introduced incident blackbody radiation of T instead(i.e., case 2 in ° 4.2) for T \ 30 K. The results areinit
shown ininit
Figure 7.
Substituting T \ 30 K and i (T \ 30 K) \ 2.1 ] 10~1
init (14) and (17),P we have
cm2 g~1 in equationso\ 1.8 ] 10~9 g cm~3 ,"/!
o\ 8.3 ] 10~16 g cm~3 .
qF1
In contrast to the result of T \ 5 K, ois greatly
init cloud experiences
qF1
exceeded by o
. The collapsing
qD1"/! in the evolution far before the radiativeat the early stage
cooling rate is balanced with the gas compressional heatingrate. However, Figure 7b shows isothermal evolution evenwhen o [ o . Hence oand olose their roles as thec between
qF1 the isothermal"/!
qF1adiabatic phases. The
criterion
andadiabatic phase is set up when the compressional heating
rate, ! , becomes comparable to the cooling rate due to theg di†usion, " . If the internal energy density ofradiative
material per unit mass,dife, is larger than the radiation energy
density per unit volume, E, " can be written as (E/o)/t ,dif t 4 i ol2/c denotes thedifotherwise " equals e/t , where
dif
dif
dif
P
radiative di†usion time with l for the length scale of thesystem. Thus the criterion is now expressed by

A B

E 1
,
(18)
! \ " 4 min e,
g
dif
o t
dif
for low-mass star formation, where e ? E/o, min (e, E/o) isalways replaced by E/o. (Otherwise it should be replaced bye, which leads eq. [18] into a simpler criterion t \ t .)dyn
difVol. 495

FIG. 8.ÈThermal evolution at the center of a collapsing cloud core formodel M1b (larger opacity). See also the legend of Fig. 4.

Substituting 2R for l by recalling the argument aboveb has
equation (17), one
pT 4
,
(19)
i o2R2Pbwhere E 4 aT 4 has been used since we consider the situation where the optical depth has already exceeded unity.
Combining equation (18) with equations (13), (15), and (19),
we have"

dif\

A B A B
B A B

C ~2@3 C ~4@32
o \ 1.2 ] 10~13 g cm~3 1
dif
1.0
0.75
]

A

i
~2@3 T
4@3P
init.
1.6 ] 10~2
10 K

(20)

For T \ 30 K and i (30 K) \ 2.1 ] 10~1 cm2 g~1,
init
Po \ 9.2 ] 10~14 g cm~3 .
difFigure 7a clearly shows the turning point between the isothermal and adiabatic phases at o \ o . In great contrast
dif
to the case of T \ 5 K, where thec transitional
phase interinit
rupts the isothermal stage before q D 1, the isothermal
phase lasts even after the optical depth exceeds unity.
6.4. Dependence on Opacity
We consider here the inÑuences caused by the ambiguityof the dust opacity for submillimeter and millimeter wavelengths. As described in ° 4.3, the three models of di†erent
opacities are examined. We provide models M1b and M1c,that correspond to models b and c in ° 4.3.

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

359

them in the numerical results. We will show in this sectionthat they can be derived from semianalytical considerations.
Possible interpretations of the observational YSOs are also
discussed, on the basis of the accretion luminosities and
SEDs in our results.
7.1. Analytical Derivation of the Core Radius
We approximate the Ðrst core as a polytropic sphere in
equilibrium. Although 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

A B

1 d
dh
m2
\ [hn ,
m2 dm
dm

(21)

where we follow the standard notation of Chandrasekhar(1967). In the present calculations the ratio of speciÐc heats,c, is Ðxed to 5/3, i.e., n \ 1.5. The radius of a polytropic
sphere in equilibrium, R , isLE(n ] 1)K
1@2
o1@n~1m ,(22)
R (o ) \
c
1
LE c
4nG

C

D

where m \ 3.65375 for n \ 1.5 (Chandrasekhar 1967) and1 which relates with the entropy that the system
K 4 p/oc,
has. One immediately obtains K as follows :FIG. 9.ÈThermal evolution at the center of a collapsing cloud core formodel M1c (smaller opacity). See also the legend of Fig. 4.

The opacity of model b gives i \ 7.7 ] 10~2 cm2 g~1,
which is larger by a factor of 5 thanP model a. Equations (14)and (17) for this case yield\ 3.2 ] 10~13 g cm~3 ,o"/!
o\ 1.9 ] 10~14 g cm~3 .
qF1
Considering o\o, one might preferqF1"/!
o \ 4.2 ] 10~14 g cm~3
dif
as the central density where the isotherm is connected withthe adiabat. However, under situations where ois closeenough to o
, as in this case, orepresents qF1
the critical"/! the compressionalqF1
density because
heating rate is no longer
negligible when the optical depth reaches unity as seen in
Figure 8.
For model c, where the dust-to-gas ratio is decreased by afactor of 10 from model a, i equals 1.6 ] 10~3 cm2 g~1,
P value of model a. Equations
which is of course a tenth of the(14) and (17) giveK \ c2 o1~c ,
(23)
s0 ad
where c is the sound speed at the initial isothermal phaseand o s0
is deÐned as the intersection of extrapolated lines ofadthe isotherm
and the adiabat in the log oÈlog T plane (asillustrated in Fig. 10). The log oÈlog T planes of numericalsolutions are shown in Figures 4, 6, 7, 8, and 9, which readthe typical value of 10~13 g cm~3 for o as tabulated inad
Table 1. Substituting equation (23) into equation(22), oneobtains the analytical description of R as follows.
fco~1@6 T
1@2cinit
R (o ) \ 5.3 AU
fc c
10~9 g cm~3
10 K

A

B A B

]

A

(24)

Figure 5 delineates R (o ) in equation (24) (dot-dashedfc toc be 10 K and 10~13 g cm~3,line), with T and o set
init
adrespectively. One will Ðnd that the analytical prediction is inexcellent accordance with the numerical results.

o\ 1.4 ] 10~16 g cm~3 ,"/!
o\ 4.3 ] 10~11 g cm~3 .
qF1
Figure 9 shows that the transitional phase appears again
where o\o \o(see ° 6.3.1). It would be worth"/!the transitionalc
qF1 phase also appears for the lower
noting that
opacity model in Figure 3 in Boss (1988).
7. DISCUSSION

In the previous section, we noted the importance of determining the radius and mass of the Ðrst core and deÐned

B

o
~1@3
ad
.
10~13 g cm~3

FIG. 10.ÈDeÐnition of o (see text)
ad

360

MASUNAGA, MIYAMA, & INUTSUKA

Apparently o relates with o, o , and/or o .
ad"/! qF1
dif
Because all of o
, o , and o depend solely on T and
"/!
qF1
dif
init
i as seen in equations (14), (17), and (20), R also depends
Pfc
only on T and i when o is given, which supports theinitPcnumerical results that R is independent of the mass of the
fc
parent cloud core (see ° 6.2).
Interestingly, o seems never to depart from 10~13 gad
cm~3 further than an order of magnitude for any models
examined here, in contrast to oand othemselves,"/!
qF1
which vary quite sensitively to the initial temperature.
Moreover, considering the weak dependences of R on o ,fc
c
T , and o in equation (24), one will Ðnd that the typicalinit
advalue of 5 AU for R would be a reasonable estimate forfc
any cases unless extreme conditions for T
and i areinit of theP Ðrstsupposed. Indeed, Figure 5 shows that the radiuscore converges to several AU in every model.
7.2. Estimate of the Accretion L uminosity
The emergent Ñux of the radiation is the most important
quantity when we compare our results to 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 :
GM M0
fc ,
(25)
R
fc
where M and M0 are the mass of the Ðrst core and the massaccretionfcrate, respectively. For determining L acc, M andfc
fc
M0 are estimated as follows.
L acc \
fc

A

1 dh
n
M \ 4no R3 [
fc
c fc
m dm

B

,

(26)

m/m1

c3
M0 \ a s0 .
G

(27)

We deÐne M as the mass of an equilibrium polytropic
fc radius of R in equation (26), where
sphere of the\ 5.5642 ]fc10~2 for n \ 1.5 (Chan[(1/m)(dh /dm)]
1nm/m
drasekhar 1967). In equation (27) a is a dimensionless coefficient, which is found to be D10 by comparison with thenumerical results for model M1a. For another example, Shu(1977) found a \ 0.975 in the expansion-wave solution. Thisdi†erence arises because the expansion-wave solution
assumes an initially hydrostatic singular sphere in contrast
to our homogeneous density distribution. Note that themass accretion rate would become smaller with time inmuch later stages. This aspect will be discussed in detail in
our subsequent papers.
Substituting equations (26) and (27) into equation (25), we
obtain

A BA
A B A

a
L acc \ 0.07 L
fc
_ 10
]

B
B

o
2@3
c
10~9 g cm~3

T
5@2
o
~2@3
init
ad
.
10 K
10~13 g cm~3(28)

Table 1 tabulates L acc evaluated by equation (28) as well as
the numerical result,fc L . By comparing L acc with L , it is seen
fc
that equation (28) o†ers a reasonable approximation
for theluminosities in most of the models.

Vol. 495

7.3. Case for c \ 7/5
Hydrogen molecules in molecular cloud cores of D10 K
have the ratio of speciÐc heats, c, equal to 5/3 because rotational and vibrational levels of H are not excited at low
2
temperatures. However, the rotational degree of freedom
raises the ratio of speciÐc heats up to 7/5 before the temperature reaches several hundred kelvin (Black & Bodenheimer 1975). The calculations presented in this paper
suppose c \ 5/3 throughout the evolution, for simplicity.
To study this e†ect, another calculation was performed
under the same conditions as model M1a, except that c is
set to 7/5 instead of 5/3. We found that R does not show
fc
apparent di†erences between the two cases, and M isfc
smaller by only a factor of 2 for c \ 7/5. These results are
consistent with past works where the rotational degree of
freedom was taken into account (Larson 1969 ; Tscharnuter& Winkler 1979). Therefore this simpliÐed treatment for cremains useful for this kind of study.
However, we will use a more realistic equation of state in
the subsequent papers, where we consider the dissociation,
ionization, and the rotational and vibrational degrees of
freedom for molecules, as well as the nonideal e†ects of theequation of state such as the pressure ionization and thedegeneracy of electrons.
7.4. Comparison with Observations
Our nongray radiation hydrodynamic calculations
provide information on the frequency-dependent radiativeÑux. Recently a number of SEDs of young stellar objects
became available. In this section we examine possible suggestions and interpretations of the observations from ourresults.
In Figure 3 we presented SEDs for the typical model,
whose shapes were found to be quite insensitive to thethermal evolution in the center throughout the Ðrst collapse
because the envelope is opaque. In Figure 12, 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 the shape of the SEDs, it isdifficult to distinguish cloud cores that have already formedthe Ðrst (or possibly more evolved) cores from starless coresor precollapse cores. However the total luminosity should
experience a considerable increase even at this early stage ofthe collapse that is hidden in the opaque envelope.
Ward-Thompson et al. (1994) presented SEDs and luminosities for four starless cores, on the basis of the results of a
submillimeter continuum survey. Their masses are D1È2
M and are comparable with our models. The SEDs (Fig. 5
in _Ward-Thompson et al. 1994) can be Ðtted by graybody
spectra of T D 10È20 K, which are quite similar to our
results in Figures 3 and 12. Furthermore, the estimatedluminosities of the observed starless cores, 10~1 L [ L [_ that2 L , seem to be rather higher than expected for cores_
have not yet experienced dynamical collapse. Indeed,
Figure 11 depicts the evolutions of the luminosities, whichreach at most D10~1 L even at the end of the Ðrst col_ a relatively luminous (D1 L )lapse. These results show_
star-forming cloud core that might be recognized to havealready 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 that are hidden by opaque envelopes, rather than as

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

361

FIG. 11.ÈSEDs that would be observed with the beam size of 103 AU for di†erent models. The SED for model H1a has been omitted because it supposesa surrounding isotropic radiation Ðeld.

preprotostellar cores as concluded by Ward-Thompson etal. (1994).
Boss & Yorke (1995) calculated SEDs at the Ðrst collapse
stage, which is ““ class [I ÏÏ stage by their three-dimensional
numerical code. No essential discrepancy is recognized
between their SEDs and ours, although the SEDs of thepole-on view by Boss & Yorke (1995) show somewhat more
excess of the component of higher temperatures. (Note thatthe vertical axis in Figs. 1 and 2 in Boss & Yorke 1995
denotes lL by their deÐnition, but F by our deÐnition.)l
Thus the pole-onview of young stellar lobjects would possibly reveal the dramatic evolution at the center of the cloud
cores.
8. SUMMARY

In this paper, we have investigated the properties of theÐrst collapse and the Ðrst core by means of radiation hydrodynamic numerical calculations, supplemented by analytical considerations. Our Ðndings are summarized as follows.
The Ðrst core Ðnally has the mass of D0.05 M and the_radius of D5 AU for T \ 10 K and the typical opacity.
If
init
the initial temperature and opacity are given, these values
for mass and radius are uniquely determined and are independent both of the mass of the parent cloud core and of theinitial density proÐle.
The thermal evolution at the center of a collapsing cloudcore is characterized by the conditions of ! \ " , q D 1,th
and/or ! \ " (see eqs. [14], [17], and [20],g respectively),gdifdepending sensitively on the initial temperature andopacity. These expressions determine how and when theearly isothermal evolution changes to nearly adiabatic evolution.

362

MASUNAGA, MIYAMA, & INUTSUKA

Vol. 495

FIG. 12.ÈEvolutionary curves of the luminosity. L vs. o (upper panel) and L vs. t (middle and lower panels) for models M1a, M2a, M3a, and M4a. Thecluminosity remains at a constant value since the collapse is initiated(i.e., t \ 0) until the luminosity begins to increase.

We have found the analytic expression for the radius ofthe Ðrst core (eq. [24]) and the accretion luminosity (eq.
[28]) as functions of the central density. They are weakly
dependent on the initial temperature and opacity. Those
expressions are useful for understanding the numericalresults.
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, ofcourse, does not exclude the possibility that the Ðrst cores
might be observed on the pole-on view, where the envelopematerial would actually be more depleted than predicted byspherically symmetrical calculations. The luminosity
increases up to D0.1 L at the end of the Ðrst collapse_
phase.
Dense cloud cores with luminosities higher than D0.1 L_

possibly have young protostars inside the opaqueenvelopes, even if observations show no evidence for protostellar activities.
In the subsequent paper, we will discuss the evolution
after the Ðrst collapse, in which the dissociation of hydrogenmolecules causes the second collapse. The main goals will
be to examine the properties of the main accretion phaseand to answer the question : ““ When does a protostar
emerges before us in near-infrared or optical light ? ÏÏ
The authors greatly appreciate valuable comments andhelp of the referee, Alan P. Boss, concerning improvements
to the paper. The research of H. M. is supported in part by
the Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists.

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

363

APPENDIX A
THE COMPUTATIONAL METHOD IN SPHERICAL SYMMETRY
A1. THE TRANSFER EQUATION SOLVER

We solve equation (6) numerically, dropping the time derivative term, in the transfer equation solver to determine thevariable Eddington factors. The transfer equation to be solved is reduced toLjl \ [h ,
l
Lq
l

(A1)

Lh
l \ [j ] S ,
l
l
Lq
l

(A2)

where
j 4 1 [I (]k) ] I ([k)] ,
(A3)
l 2 l
l
h 4 1 [I (]k) [ I ([k)] .
(A4)
l 2 l
l
Note that an independent variable l is noted as a subscript hereafter. Figure 13 illustrates the paths along which equations(A1) and (A2) should be integrated. The speciÐc intensity on each point at which an integral path intersects a radial grid is
integrated over k on each ring to calculate the Ðrst three moments for radiation,

P
P
P

1
j dk ,
(A5)
l
0
1
(A6)
H \ kh dk ,ll
0
1
K \ k2j dk ,
(A7)
l
l
0
where k is the direction cosine, i.e., k \ cos h, where h is the angle with respect to the outward normal (see Fig. 1). The Ðrst two
moments for radiation, J and H , correspond to the monochromatic radiation energy density and radiation Ñux, respectively.
l radiation pressure, which is the rr-component of the radiation stress tensor in spherical
The third, P or (4n/c) Kl , is thellcoordinates. Both the hh and // components equal P [ (3P [ E )/2 (see MM).
l
l
l
J \
l

FIG. 13.È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, k, is deÐned as k \ cos h, where h is the angle between the lines oxll
and or.

364

MASUNAGA, MIYAMA, & INUTSUKA

Vol. 495

The scalar variable Eddington factor is
K (r)f (r) \ l .
l
J (r)
l
When the radiation Ðeld is nearly isotropic, i.e., I(k) \ I ] I k, K is found to equal 1 J using equations (A5) and (A7). On the0
13other hand, the streaming limit condition I(k) \ Id(k [ 1) yields K \ J. The Eddington factor thus varies between 1 and 1,
3
except for some special cases where f goes down below 1 , according to the degree of anisotropy of the radiation Ðeld.
3
We adopt the monochromatic transfer equation solver in our scheme. It has two advantages over the gray transfer equation
solver. First, the energy mean and Ñux mean opacities are i and i for the moment equations are obtained (see ° A3), whereas
EF
the Planck mean or Rosseland mean opacities are substituted in the gray approximation. Another advantage is that one canobtain SEDs directly from the monochromatic radiation Ñuxes. The obtained SEDs will help us to compare theoreticalmodels to observations, as mentioned in ° 2.
A2. THE MOMENT EQUATIONS SOLVER

In spherical symmetrical Lagrangian (comoving) coordinates, equations (1)È(4), (7), and (9) lead to the following (MM) :
dr
\v ,
dt

(A8)

3 LM
r,
4n L(r3)

(A9)

o\

Dv
GM
L(p ] Q) sr [ 4nr2
\[
] FF,
(A10)
Dt
r2
LM
oc
r
E
D 1
v
D 1
L(4nr2F)De]] f[ (3f [ 1)E ] (p ] Q)
\v [
,
(A11)
CR
o
Dt o
or
Dt o
LM
Dt
r
s F
1 DF 4nr2o L( fqE)]\[ F ,
(A12)
c2 Dt
q
LM
c
r
D 1
v
1
L(4nr2F)
D E] f[ (3f [ 1)E \ (4ns B [ cs E) [
,
(A13)
P
E
Dt o
or
o
LM
Dt o
r
where Q is the artiÐcial viscosity in the formulation of von Neumann & Richtmeyer (1950), and the sphericality factor q isdeÐned by (Auer 1971)

A B C ABD

AB C ABD

AB

P

r 3f [ 1
dr@ .
fr@
Rin
These equations are closed when the scalar Eddington factor and the equation of state (5) are given. The internal energy
density e and the gas pressure p are replaced by T and o using the equation of state, which is mentioned in ° 4.4.
The six variables r, o, v, T , E, and F are solved as functions of mass. For spatial di†erencing of the equations, o, T , and E aredeÐned at the center of each grid cell, and r, v, and F are on the interfaces that bound the cell. We thus obtain vectors of thevariables (r , o,v,T,E, F ), which are to be solved in a fully implicit scheme. The system is linearized around thei i`1@2
trial estimates(r*, o*i i`1@2, v*, T *i`1@2, E*i, F*) to obtain an (I ] I) block tridiagonal system of (6 ] 6) matrices, where I isi cells.
i`1@2Thei solutioni`1@2 vectors
i`1@2 (dri , dothe number of grid, dv , dT, dE, dF ) are solved by Gaussian elimination.
i are
i`1@2constrained
i
i`1@2
i
In contrast to implicit method, explicit methods
by thei`1@2
CFL condition,
which assures numerical stability(Courant, Friedrichs, & Levy 1928),
ln q \( o v o ] c)*t ¹ *r ,(A14)where v, c, *t, and *r are the Ñuid velocity, the speed of transport of information, the time step, 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 compared to the evolutionary timescale because c in (A14) is the speed of light in this case. Then, the moment equations ofradiation (A12) and (A13) with the Ñuid energy equation (A11) are often solved on an implicit scheme while the Ñuid motion
can be solved explicitly for some cases. However, the CFL condition that rules the equation of motion (A10)Èi.e., equation(A14) in which c is now replaced by the speed of soundÈstill imposes such a severe restriction that explicit methods are ruledout in numerical calculations for the whole evolution of protostellar collapse (Tscharnuter & Winkler 1979). It is therefore
imperative that the full set of equations (A8)È(A13) is solved with an implicit scheme.
The frequency-averaged absorption coefficients s , s , and s and Eddington factor f are needed in order to solve
P , are given by direct integrationEequations (A8)È(A13). The weighting values E and FE, orFJ and H
of I over k. We explain
l
l
l
l
l
the method for calculating J and H afterward. The frequency-averaged values are obtained using the solved
spectral proÐles
ll
J and H asll/= s H dl
s \ 0 l l ,
(A15)
F
H

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

365/= s J dl
s \ 0 l l ,
E
J

(A16)

/= s B dl
s \ 0 l l ,P
pT 4/n(A17)/= f J dl
f \ 0 l l ,
E
J

(A18)

and

where
J\
H\

P
P

=

0
=

J dl ,l

H dl .
l0
The entire procedure at each time step is carried out as follows.
1. Trial estimates of density and temperature structures [o(r) and T (r)] are assumed.
2. Equations (A1) and (A2) are solved under the given density and temperature distributions to calculate the Ðrst threemoments for radiation equations (A5), (A6), and (A7).
3. The integrations over frequency equations (A15)È(A18) are carried out to obtain f (r), i (r), i (r), and i (r).
E
E and Fthen yieldPimproved density
4. Equations (A8)È(A13) are solved using the given Eddington factor and mean opacities
and temperature distributions, with which the transfer equation is solved again.
5. We iterate steps (1)È(4) until it converges sufficiently.
One thus obtains the radiation hydrodynamic quantities that are consistent with the radiation Ðeld at each time step. Trial
estimates of an iterative procedure are given by the solutions at the previous time step.
A3. THE ENERGY MEAN AND FLUX MEAN OF OPACITY

In the transfer equation solver, the mean opacities (i , i , and i ) are calculated, as well as the Eddington factor. We note
E calculatingF
P the mean opacities in equations (A16)È(A17).
here how to obtain J (r) and H (r), which are required forll
MM suggested that J (r) and H (r) should be obtained by solving the monochromatic moment equations of radiation using
the current estimates of lT (r) and f l(r) (MM p. 484). However, one can obtain J (r) and H (r) directly also by integrating I (r, k)llll
and kI (r, k) over k. The former method,
however, has a serious problem for nongray
calculations
as described below.
l
In a nongray calculation the Planck function, B , increases exponentially with temperature at a frequency in the Wienlregime, in contrast to the gray approximation, for which
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, anlevolved protostar has a central core of T Z 3000 K with a cold envelope of T D 10 K. Such a situation involves
a great rangeof I over 100 orders of magnitude at l D 1014 Hz throughout the computational region. It exceeds the available numericall by far for ordinary computational tools. Thus, J , H , and K in the cold region at high frequency should equal zerorange
l obtainedlnumerically. In such cases the Eddington factor is notl easily
because f 4 J /K is not well-deÐned numerically,ll lvalues of J and K are to zero. If
although the actual Eddington factor should behave regularly however close the individualll
the Eddington factor is not obtained, one cannot solve the moment equations.
In contrast to the method suggested by MM, we calculate the moments of radiation directly by integrating I (r, k) over kl even in theinstead of solving the monochromatic moment equations in the transfer equation solver. Our method works wellcase where J and K are very small.
ll
A4. COORDINATES FOR COMPUTATION

Calculations are performed on the 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. This method was found to successfully
resolve the accretion shock at the surface of the Ðrst core in our calculations. For Eulerian calculations we introduce
second-order accurate advection, supplementing the Lagrangian RHD equations (A8)È(A13).
For discretization we use 200 radial grid cells and 20 frequency bins. The paths along which the transfer equation is
integrated are set tangentially to each radial grid (see Fig. 13). The innermost zone has 30 additional paths, so as to provide a
sufficient number of angle-mesh points for innermost rings of the radial grid. The initial radial grid is set up as * log
r \ constant.
APPENDIX BTEST CALCULATIONS FOR RADIATIONTo check the abilities of the moment equations solver, the transfer equation solver, and the convergence between them, test
calculations for the numerical code are presented brieÑy in this section. Because the moment equations of radiation are366

MASUNAGA, MIYAMA, & INUTSUKA

Vol. 495

difficult to solve analytically when they are coupled with hydrodynamics, e†ects of radiation are investigated separately fromhydrodynamics here (see Ensman 1994). Test calculations for hydrodynamics are not presented in this paper, but we havechecked that numerical results of some hydrodynamic tests show good agreements with analytical predictions.
For radiation tests, material is assumed to be static, with the right-hand side of equation (A10) set to zero. For simplicity,
the density distribution is assumed to be homogeneous and radiation to be gray. Then, suppose small Ñuctuations oftemperature in radiative equilibrium. The Ñuctuation should be smoothed out by radiation quasi-statically if the radiationenergy density is sufficiently small compared to the internal energy density of material. This mode is called the ““ thermal
relaxation mode ÏÏ of radiation (MM p. 510). The relaxation rate n(k) for a perturbation with wavenumber k is (Spiegel 1957)

C

D

s
s
n(k) \ l 1 [ 0 cot~1 0 ,
k
k

(B1)

where
16ps T 3
0 0,
(B2)
oc
v
where c is the heat capacity of material and subscript zero represents unperturbed quantities.
v
Spiegel (1957) derived equation (B1) for a perturbation of plane wave, i.e., eik Õ x. However a spherical wave solution isrequired for the test of our numerical scheme. We have found that the zeroth-order spherical Bessel function, j (r) 4 sin kr/kr,
0 Appendix C
is the eigenfunction in the spherical case, and it yields the exactly same dispersion relation as equation (B1). Seefor derivation.
Figure 14 shows relaxation of a Ñuctuation of the spherical wave j (kr) in the case of the total optical depth of 10. The
0unperturbed temperature distribution is homogeneous at 10 K. A temperature
disturbance (Fig. 14a) damps toward theunperturbed state. The radiation temperature T 4 (E /a)1@4 (Fig. 14b) damps correspondingly to the relaxation of therad moreradmaterial temperature, but slightly more rapidly. For
opaque material, of course, the di†erence between T and T wouldradbe smaller.
The dispersion relations are shown in Figure 15. The ordinate represents n(k)t , where t 4 (cs)~1 is the mean Ñight time ofjj a photon travels over a mean freea photon. The ordinate therefore means how much a temperature disturbance decays
while
path. First, to test the moment equations solver solely, we calculate the thermal relaxation mode in the Eddington approximation, Ðxing f 4 1 instead of solving the transfer equation solver. In this case equation (B1) reduces to the dispersion
E Eddington
3
relation under theapproximation (Unno & Spiegel 1966),l\ln (k) \.
(B3)
E
1 ] 3(s /k)20
The numerical results show good agreement with equation (B3) in Figure 15a. Therefore the moment equations are found tobe solved correctly.
Second, accuracy in the transfer equation solver is tested. In Figure 15b, the numerical result that was obtained using the
Eddington factors calculated by the transfer equation solver is plotted with the analytical dispersion relation (eq. [B1]). They
agree, again, with each other within a relative error of one percent. These results guarantee the accuracy of the method over awide range in optical depth.

APPENDIX CTHERMAL RELAXATION MODE IN SPHERICAL SYMMETRY
The thermal relaxation mode of a radiating Ñuid was investigated Ðrst by Spiegel (1957). He derived the dispersion relation(B1) for a perturbation of a plane wave, which inevitably requires that the temperature distribution must be in slab symmetry.
We will show that a temperature disturbance in spherical symmetry decays in the same rate as given by equation (B1).
We assume that there are no Ñuid motions and no heat sources, and the radiation Ðeld is gray. The Ñuid energy equation (3),then, reduces tooLe\ 4ns(J[B) ,
Lt

(C1)

where J 4 (1/4n) { I d) \ (c/4n)E is the mean intensity and s is a gray opacity per unit volume. The source function in LTE,B, can be replaced by (p/n)T 4.
The speciÐc intensity I is determined by the transfer equation,LI\ [s(I[B) ,Ls
where we have neglected the transit time of radiation.

(C2)

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

367

FIG. 14.ÈResults of test calculations for the thermal relaxation mode. (a) Relaxation of a small disturbance of temperature. The given perturbation (thecurve of the largest amplitude) decays to a homogeneous distribution. (b) Relaxation of radiation temperature. (c) Relaxation of radiative Ñux.

For linearizing the above two equations, we write the temperature as follows, neglecting the terms of higher orderÑuctuation.
T (r, t) \ T (r) ] T (r, t) .
0
1

(C3)

LB
B(r, t) \ B (r) ] 0 T 4 B ] B ,
0
0
1
LT 1

(C4)

Ls
s(r, t) \ s (r) ] 0 T 4 s ] s .
0
0
1
LT 1

(C5)

The other quantities are written as

Using equations (C3), (C4), and (C5), equations (C1) and (C2) are reduced to the following equation, which describes therelaxation of the temperature disturbance (MM p. 509).

368

MASUNAGA, MIYAMA, & INUTSUKA

Vol. 495

FIG. 15.ÈDispersion relation for the thermal relaxation mode of radiation. See text for the physical meanings of the axes. T op : The dispersion relation
under the Eddington approximation. Solid line is the analytical solution, and open circles denote the numerical results with f 4 1 . Bottom : The exact
dispersion relation without the Eddington approximation. Dotted line corresponds to the exact solution, and open circles denote 3the numerical results
obtained by using the variable Eddington factor.

C

Q P

D

=
s ds T (r Ô ns)e~s0s ,
(C6)
0
1
0
where n is the unit direction vector along a light ray and l is deÐned in equation (B2).
The terms including s have vanished because of the assumptions that the unperturbed radiation Ðeld is isotropic and inradiative equilibrium. All1the assumptions introduced above are the same as those by Spiegel (1957).
Equation (C6) admits separable solutions of the following form for spherical symmetric cases,
LT (r)
1
1 \ [l T (r) [
1
Lt
4n

d)

T (r, t) \ /(k, t) j (kr) ,(C7)
10
where j (x) 4 sin x/x is the 0th order spherical Bessel function.
0
Substituting
equations (C7) and r@ 4 r Ô ns into (C6) rewritten in the cylindrical coordinates (r, k), where k is the direction
cosine,C

1L/\ [l 1 [
2
Lt

P P
1

~1

dk

=

D

j (kr@)
dq e~q 0/,
j (kr)
0
0(C8)

where q 4 s s.
0
For the purposeof replacing r@ with q and k, let the path of a light ray s be written down explicitly (see Fig. 16),s \ ^[r@2 [ r2(1 [ k2)]1@2 ] rk .

(C9)

The negative sign in the Ðrst term in right-hand side is required only for r@ \ r. After rearranging equation (C9) and replacing swith q, we havem 4 kr@ \k[q2 [ 2qs rk ] (s r)2]1@2 ,
0
0
s
0

(C10)

No. 1, 1998

AN RHD MODEL FOR PROTOSTELLAR COLLAPSE. I.

369

FIG. 16.ÈIllustration of the path of integration, s. See Appendix C.

ands m
dk \ [ 0 dm ,qrk2kk(q ] s r), m 4 m(k \ 1) \ o q [ s r o .
0
1
0
s
s
0
0
After integrating over m instead of k, the last term in equation (C8) equalsm~1

4 m(k \ [1) \PP= e~q cos m [ cos ms=sin (k/s q)
1~1 \ 00 .
dq
dq e~qqj (kr)kq
0
0
0
The last integral is reduced to cot~1s /k). Hence, equation (C8) becomes0ssL/\ [l 1 [ 0 cot ~1 0 / .
k
k
Lt
s
0
2rk2

C

A BD

(C11)

The thermal relaxation rate n(k) is deÐned asL/\ [n(k)/ .
Lt

(C12)

Comparing equations (C11) and (C12), we have equation (B1).
REFERENCES
Adams, F. C., & Shu, F. H. 1986, ApJ, 308, 836
Mizuno, A., Onishi, T., Hayashi, M., Ohashi, N., Sunada, K., Hasegawa,Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122
T., & Fukui, Y. 1994, Nature, 368, 719
Appenzeller, I., & Tscharnuter, W. 1974, A&A, 30, 423
Narita, S., Nakano, T., & Hayashi, C. 1970, Prog. Theor. Phys., 43, 942
Auer, L. H. 1971, J. Quant. Spectrosc. Radiat. Transfer, 11, 573
Penston, M. V. 1969, MNRAS, 144, 425
Black, D. C., & Bodenheimer, P. 1975, ApJ, 199, 619
Pollack, J. B., Hollenbach, D., Beckwith, S., Simonelli, D. P., Roush, T., &Bodenheimer, P., & Sweigart, A. 1968, ApJ, 152, 515
Fong, W. 1994, ApJ, 421, 615
Boss, A. P. 1984, ApJ, 277, 768
Preibisch, T., Sonnhalter, C., & Yorke, H. W. 1995, A&A, 299, 144
ÈÈÈ. 1988, ApJ, 331, 370
Shu, F. H. 1977, ApJ, 214, 488
Boss, A. P., & Myhill, E. A. 1992, ApJS, 83, 311
ÈÈÈ. 1995, in Molecular Clouds and Star Formation, ed. C. Yuan &Boss, A. P., & Yorke, H. W. 1995, ApJ, 439, L55
J. H. You (Singapore : World ScientiÐc), 97
Buchler, J. R. 1979, J. Quant. Spectrosc. Radiat. Transfer, 22, 293
Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23
Castor, J. I. 1972, ApJ, 178, 779
Spiegel, E. A. 1957, ApJ, 126, 202
Chandrasekhar, S. 1967, An Introduction to the Study of Stellar Structure
Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637(New York : Dover)
Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 (SMN)
Chick, K. M., Pollack, J. B., & Cassen, P. 1996, ApJ, 461, 956
Tscharnuter, W. M., & Winkler, K.-H. 1979, Comp. Phys. Comm., 18, 171
Courant, R., Friedrichs, K. O., & Levy, H. 1928, Math. Ann., 100, 32
Unno, W., & Spiegel, E. A. 1966, PASJ, 18, 85
Efstathiou, A., & Rowan-Robinson, M. 1990, MNRAS, 245, 275
von Neumann, J., & Richtmeyer, R. D. 1950, J. Appl. Phys., 21, 232
ÈÈÈ. 1991, MNRAS, 252, 528
Ward-Thompson, D., Scott, P. F., Hills, R. E., & Andre, P. 1994, MNRAS,
Ensman, L. 1994, ApJ, 424, 275
268, 276
Goldsmith, P. F., & Langer, W. D. 1978, ApJ, 222, 881
Winkler, K.-H. A., & Newman, M. J. 1980a, ApJ, 236, 201 (WNa)
Inutsuka, S., & Miyama, S. M. 1997, ApJ, 480, 681
ÈÈÈ. 1980b, ApJ, 238, 311 (WNb)
Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676
Yorke, H. W. 1979, A&A, 80, 308
Lada, C. J. 1987, in IAU Symp. 115, Star Forming Regions, ed. M. PeimÈÈÈ. 1980, A&A, 85, 215
bert & J. Jugaku (Dordrecht : Reidel), 1
Yorke, H. W., Bodenheimer, P., & Laughlin, G. 1993, ApJ, 411, 274
Larson, R. B. 1969, MNRAS, 145, 271
ÈÈÈ. 1995, ApJ, 443, 199
Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York : Oxford Univ. Press) (MM)