THE ASTROPHYSICAL JOURNAL, 531 : 350È365, 2000 March 1( 2000. The American Astronomical Society. All rights reserved. Printed in U.S.A.

A RADIATION HYDRODYNAMIC MODEL FOR PROTOSTELLAR COLLAPSE. II. THE SECONDCOLLAPSE AND THE BIRTH OF A PROTOSTARHIROHIKO MASUNAGA1 AND SHU-ICHIRO INUTSUKA2
Division of Theoretical Astrophysics, National Astronomical Observatory, Mitaka, Tokyo 181-8588, Japan
Received 1999 April 19 ; accepted 1999 November 15

ABSTRACT
We carry out radiation hydrodynamic calculations to study physical processes in the formation of a 1
M protostar. Following our previous work, calculations pursue the whole evolution from the beginning_of the Ðrst collapse to the end of the main accretion phase. The adiabatic core formed after the initialcollapse (i.e., the Ðrst core) experiences further gravitational collapse triggered by dissociation of molecular hydrogen, which leads to the formation of the second core, i.e., the birth of a protostar. The protostar
grows in mass as accreting the infalling material from the circumstellar envelope, while the protostar
keeps its radius at D4 R during the main accretion phase. These typical features in the evolution are in_good agreement with previous
studies. We consider two di†erent initial conditions for the density distribution : homogeneous and hydrostatic cloud cores with the same central density of 1.415 ] 10~19 g
cm~3 . The homogeneous core has the total mass of 1 M while the hydrostatic core has 3.852 M . For_ rapidly rises to the maximum value of _25 L
the initially homogeneous model, the accretion luminosity_
just after the birth 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 di†erence
arises because the mass accretion rate varies depending on the inward acceleration at the initial stage,
which a†ects the luminosity curve. A less massive hydrostatic core would possess the similar properties
in the luminosity curve to the 3.852 M case, because a hydrostatic cloud core with mass lower than_
3.852 M can be shown to provide a smaller
mass accretion rate after the birth of a protostar and a_more gradual rise in the luminosity curve. Our numerical code is designed to provide the evolution of
the spectral energy distribution (SED) along with the dynamical evolution in our spherically symmetriccalculations. We conÐrm that the SED evolves from a 10 K graybody spectrum to hotter spectra typicalfor class I and Ñat spectrum sources. The SED for the class 0 sources corresponds to the age of 2 ] 104
yr, which is smaller by an order of magnitude than the typical age of class I objects. Considering possible nonspherically symmetric e†ects, we suggest that observed class 0 sources should consist of the““ genuine ÏÏ class 0 objects, which are as young as 104 yr, and more evolved protostars on edge-on view(““ class 0Èlike class I ÏÏ objects). The contamination of edge-on class I objects into class 0 sources are notnegligible because they are more abundant than genuine class 0 objects. Since observations indicate thatthe class 0 sources are typically more luminous than class I sources, the initially hydrostatic model,
where the luminosity increases monotonically with time, does not match the observations. The initiallyhomogeneous model, in contrast, shows a tendency consistent with the observations. Compiling our
results and other theoretical and observational evidence, we illustrate an evolutionary picture of protostar formation. In terms of the evolutionary time and the inclination to an observer, we Ðnd that protostellar objects are clearly categorized.
Subject headings : hydrodynamics È ISM : clouds È methods : numerical È radiative transfer Èstars : formation
1. INTRODUCTION

timescale is relatively short for such young objects, so thatthe expected number of the objects is small. Recently, observations have made remarkable progress in the studies ofprotostar formation. Andre, Ward-Thompson, & Barsony(1993) proposed the ““ class 0 ÏÏ category for the objects of
L /L [ 2 ] 104, where L
is the 1.3 mm luminosity,
bol they
1.3 considered class 0 sources
1.3
and
as candidates for young
protostars under dynamical growth. In addition to theabove deÐnition, Barsony (1994) listed some empirical characteristics for class 0 sources as follows : (1) no detection in
the near- and mid-infrared, (2) a cold spectral energy distribution (T ¹ 30 K), and (3) presence of a molecularoutÑow. While the class 0 category has become familiar, its
evolutionary position is not conÐrmed by detailed theoretical calculations. Past theoretical studies on protostar formation (Winkler & Newman 1980a, 1980b, hereafter WN
collectively ; Stahler, Shu, & Taam 1980, hereafter SST, andreferences therein) were not designed to explain the recent

Investigation of star formation processes has been a signiÐcant subject both for theoretical astrophysics and forobservational astronomy. Nevertheless, many problems are
to date unresolved, especially in the early phase in the evolution in which 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.
Objects at very early stages of the evolution are difficult
to observe for the following two reasons : (1) a growing
protostar is veiled by the infalling envelope, and (2) the

1 Present address : Earth Observation Research Center, National Space
Development Agency of Japan, 1-9-9 Roppongi, Minato-ku, Tokyo 1060032 Japan ; masunaga=ccsr.u-tokyo.ac.jp.
2 inutsuka=th.nao.ac.jp.

350

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.
observations by means of observable quantities such as
spectral energy distributions and molecular line spectra.
Now a theoretical model for protostar formation and its
spectral synthesis are required to account for latest observations. This motivates us to a series of this study, whichstarted with our previous paper (Masunaga, Miyama, &Inutsuka 1998, hereafter Paper I). We have also performedmodel calculations for line spectral synthesis, which are presented in Masunaga & Inutsuka (2000).
We have developed a radiation hydrodynamic (RHD)
numerical code in spherical symmetry, where the angle- andfrequency-dependent radiative transfer equation is solved
along with a full set of the radiation hydrodynamic equations. Therefore our calculations invoke neither the Eddington (and di†usion) approximation nor the grayapproximation. In addition, the monochromatic radiativetransfer equation provides spectral energy distributions(SEDs), and the evolution of a SED is obtained as well as
the dynamical evolution from our numerical results. This
enables us to compare computational results directly withobservations.
In Paper I, we investigated the evolution of protostellar
collapse before the formation of a protostar. Main subjectsof Paper I are the gravitational collapse of molecular cloud
cores (the so-called Ðrst collapse), violation of the isothermal evolution (discussed further in Masunaga & Inutsuka
1999a), and the formation of a central adiabatic core (theÐrst core). Previous versions of our numerical code did notinclude dissociation of molecular hydrogen, and thus theevolution succeeding to the Ðrst collapse was outside thescope of Paper I because the Ðrst core experiences further
gravitational collapse (i.e., the second collapse) triggered bydissociation of H , which leads to the birth of a protostar.
2 of the code takes into account all the
The present version
necessary physical processes (described in ° 2.1) in order topursue the entire evolution from gravitational collapse of amolecular cloud core to the end of mass accretion from theinfalling gas envelope onto a protostar. A main purpose ofthis paper is to construct a theoretical model for protostar
formation incorporating updated data for EOS and opacities that were not available in days of WN and SST.
Another important aim of the present work is to clarify
evolutionary positions of observed protostellar objects(starless cores, class 0 sources, etc.) in terms of SEDs andbolometric temperatures (Myers & Ladd 1993). Some
recent theoretical studies yielded SEDs of protostellar
objects by multidimensional RHD calculations, where thecalculation stopped before a protostar is born (the ““ class I ÏÏ
phase ; Boss & Yorke 1995) or a protostar is put as a central
luminosity source (Yorke, Bodenheimer, & Laughlin 1993,
1995). Their calculations used the Eddington approximation (Boss & Yorke 1995) or the Ñux-limited di†usion
approximation (Yorke et al. 1993, 1995) for simplicity,where the hydrodynamic quantities and radiation Ðelds are
not treated in a self-consistent way. In this work we haveestablished comprehensive theoretical models for the wholeevolution of a protostar, solving self-consistently the physical structures of the protostellar interior and the circumstellar envelope.
In ° 2.1 we explain the revisions that have been madesince our previous work in the numerical scheme and in the
input physical processes. Computational results are shownin ° 3, where the evolutionary process in protostar formation is described in terms of the hydrodynamical quantities,

351

SEDs, and luminosity curves. Section 4 is devoted to discussions on the relations between our results and statisticaltrends in observations. We there illustrate an evolutionary
picture of protostar formation by taking into accountexpected nonspherically symmetric e†ects. Conclusions are
summarized in ° 5. The basic equations that we solve innumerical calculations are shown in Appendix.
2. CALCULATIONS

2.1. Numerical Scheme and Input Physics
As a computational scheme for handling radiationhydrodynamics, we adopt the variable Eddington factor(VEF) method, where we solve the system of the hydrodynamic equations and the radiation moment equations
which are coupled with the angle-dependent radiativetransfer equation via the variable Eddington factor. The
VEF method enables us to solve the RHD equations with adesirable accuracy and with a reasonable computationalcost. In contrast to WN, who also adopted this method, we
have designed our scheme to yield SEDs as computational
byproducts by considering the monochromatic radiativetransfer. Detail on the computational procedure of the VEF
method is described in Paper I.
Although the outline of our numerical scheme remains
unchanged throughout the revisions, some signiÐcant
improvements have been made in the numerical code since
Paper I, as summarized in the following :1. Input physical processes are more elaborate in thepresent version. Dissociation of molecular hydrogen, ionization of atomic hydrogen and helium, and nonideal equation of state (EOS) have been incorporated. We use thetabulated data provided by Saumon, Chablier, & Van Horn
(1995) for the nonideal EOS. The speciÐc heat of molecularhydrogen, which was given as a constant in Paper I, is now
calculated as a function of temperature by solving the partition functions for the rotational and vibrational energy
levels of H in addition to the translational component.
2 and dust opacities have been also improved.
2. The gas
The gas opacity is taken from Iglesias & Rogers (1996) inthe high-temperature regime [log T (K) [ 3.75] and from
Alexander & Ferguson (1994) for lower temperatures. The
dust opacity was provided by H. W. Yorke (1998, private
communication), where the three components, i.e., amorphous carbon, silicate, and dirty ice-coated silicate, are considered on the basis of the opacity model by Preibisch et al.
(1993). The scattering by dust grains are taken into accountas well as the absorption.
3. The dust temperatures are solved independently of thegas temperature, while they are assumed to be the same in
WN, SST, and Paper I. The dust and gas temperatures aredecoupled with each other when the density is so low thatthe exchange of the thermal energy by collisions between
gas and dust grains are less e†ective. Such a situation isfound to occur where the infalling envelope becomes
tenuous in the late stages of the evolution.
See Masunaga (1999) for more detailed descriptions on
the input physical processes and on the numerical scheme.
2.2. Initial Density Structure
An appropriate initial condition for gravitational collapse of a cloud core is related with mechanism for the
initiation of star formation, which involves a controversial

352

MASUNAGA & INUTSUKA
09Vol. 531

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)

(d)
-1

9

20

9

-2
15

6

-3

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. 1.ÈNumerical results for the initially homogeneous model before the main accretion phase begins are shown. The succeeding evolution is depicted inFig. 4. Distributions of temperature (a), density (b), infall velocity (c), and mass (d) are presented. Thick line denotes the initial condition. Elapsed times aretabulated in Table 1.

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 1995 ;and references therein) 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 initialcondition for the dynamical collapse, while the ““ triggered
star formation ÏÏ would not necessarily. Although somerecent observations resolved spatial structures of molecular
cloud cores, the initial condition is still difficult to determine
from observational evidence.

Larson (1969, hereafter L69) and WN supposed a homogeneous distribution as the initial density structure in theircalculations for protostar formation. A homogeneous cloudcore is not in hydrostatic equilibrium and hence the gravitational collapse is initiated with a considerable inward acceleration. The dynamical evolution is known to be insensitive
to the initial density conÐguration (Bodenheimer & Sweigart 1968 ; Paper I), 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 small di†erence in the mass accretionrate, which possibly hastens or retards the growth of a protostar.
In this paper, two limiting cases are considered for the
initial density distribution : the homogeneous distribution

No. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.

and the hydrostatic distribution. The homogeneous cloudcore has a total mass of 1 M and is bounded by the Ðxedboundary at R \ 104 AU, _which yields the density of
1.415 ] 10~19g cm~3. This condition is equivalent to thoseby L69 and WN. The hydrostatic distribution is obtainedby solving the equation of hydrostatic equilibrium for an
isothermal sphere :

A B

4nGo
1 d r2 do
\[
,
(1)
c2
r2 dr o drswhere c is the isothermal sound speed. We numericallysintegrate equation (1) outwardly from the center with thesame central density as in the homogeneous model (i.e.,o 4 1.415 ] 10~19). For the initially hydrostatic models
c outer boundary is Ðxed at R \ 4 ] 104 AU and thethetotal mass is then found to be 3.852 M . We adopt a larger_
value of the total mass than in the homogeneous model
because the mass should exceed a critical value so that an
isothermal hydrostatic core is gravitationally unstable(Ebert 1955 ; Bonner 1956).
2.3. Initial T emperature Structure
For the thermal energy balance at the initial stage, we
consider a more realistic condition than in Paper I and thanin any other works concerned with protostar formation.
Previous studies assumed isotropic blackbody radiation of
10 K surrounding the cloud core (L69 ; WN) or heat sources
embedded homogeneously with no incident radiation(Paper I). In contrast, the initial temperature proÐle here istaken as the thermal equilibrium state between the thermalcooling of dust grains and the heating by the external radiation and cosmic rays (e.g., 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 6000 K with the intensity integrated over frequencyof 2 ] 10~4 ergs s~1 cm~2 for the initially homogeneousmodel and 1 ] 10~4 ergs 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 under the Ðxed density distribution. The stellar
radiation sustains the cloud core to be nearly 10 K (Hayashi& Nakano 1965) with the cosmic ray heating rate of v \CR
6.4 ] 10~28n(H ) ergs s ~1cm~3 (Goldsmith & Langer21978). The initial cloud core has an optical depth close tounity to optical photons. Incident optical photons, therefore, are attenuated as passing through the cloud core andthe temperature in deep interior is slightly lower than thetemperature near the cloud edge (see Fig. 1).
In order to consider the thermal balance between thethermal cooling and external radiation, a frequencydependent treatment for radiation transport is imperative.
The dust opacity for optical wavelengths is larger by 4
orders of magnitude than the value for the submillimeterand then gray approximations are totally invalid.
3. RESULTS

3.1. Overview
We present computational results of protostellar collapseof an initially homogeneous cloud core. The result of theinitially hydrostatic model is found to be qualitatively thesame as the initially homogeneous one except for the behavior of the luminosity curve. Hence we describe the evolu-

353

tionary properties only for the initially homogeneous modelin this section, and the luminosity curve will be discussed
for each case in ° 3.2.
3.1.1. T he Early Isothermal Phase

Figure 1 shows the evolution before a protostar is born.
Evolution at the early stage, i.e., the Ðrst collapse stage,proceeds along the line described in Paper I, where thedensity distribution is reminiscent of the isothermal selfsimilar solutions such as the so-called Larson-Penston solution (L69 ; Penston 1969, hereafter the LP solution). The
results are consistent with Foster & Chevalier (1993), where
a numerical calculation was performed for the isothermal
collapse of a marginally unstable sphere. A major di†erence
from Paper I arises in the thermal evolution because of themore realistic assumption for the external radiation Ðeld(° 2.3). In contrast to the previous models, where the temperature never descends below the initial value, temperaturein the central region once decreases to D5 K and then rises
again. The reason for the decrease of T is that increasingcolumn density of the collapsing cloud core prevents opticalphotons from penetrating into the interior. The cloud corethen cools down to 5 K, at which the cosmic ray heating
support the temperature against the thermal cooling.
However, the compressional heating rate (! ) immediately
overtakes the thermal cooling rate (4ni B) atg o B 10~17 g
P merged,c and thecm~3, where all the lines except i cE areE
central temperature begins to increase (Fig. 2). This behavior is supported by a result in our previous work (° 3.1 inPaper I), where the isothermal evolution of a cloud core of
T \ 5 K was found to cease very early ato D o
\ 8.8
c
"/!
] 10~18 g cm ~3. Cooling of gas due to collisions
with dust
grains (" in Fig. 2a) is closely coupled with the radiative
gd and with the compressional heating rate forcooling rateo Z 10~17 g cm~3.
c
3.1.2. T he First Core

Although the central temperature begins to rise far before
the optical depth arrives at unity, the evolution becomescompletely adiabatic when q B 1 at o B10~12 g cm~3, atc overtakes the radiwhich the radiative heating rate (i cE)
Eative cooling rate. For adiabatic evolutions the e†ective
adiabatic exponent, c 4 d ln T /d ln o ] 1, equals the ratio
eff 7/5 for T Z 100 K) as one can seenof speciÐc heats (c \effin Figure 2b. The gravitational collapse is decelerated by theincreasing gas pressure, which leads to the formation of theÐrst core. The physical properties of the Ðrst collapse andÐrst core were discussed in detail in Paper I. Typical scalesof radius and mass of the Ðrst core agree with the values
obtained there. We do not repeat further discussions on theÐrst collapse here.
3.1.3. T he Second Collapse

The second collapse takes place (curves 6È8 in Fig. 1)
when the central density and the central temperature reacho D10~7 g cm~3 and T D 2000 K, where dissociation ofccmolecularhydrogen reduces
the e†ective adiabatic exponent to D1.1. Dissociation of H ceases when the disso2 the center at o D 10~2
ciation degree approaches unity inc core isand T D 104 K. This is the instance when the secondcformed, i.e., when a protostar is born. The dynamical timescale of the second collapse is of the same order as thefree-fall time corresponding to o \ 10~7 g cm~3, which is
0.1 yr. The second collapse therefore begins and ends almost
instantaneously compared to the Ðrst collapse, which lasts

354

MASUNAGA & INUTSUKA

Vol. 531

2

20

(a)
1

0

R [AU]

10
-1

-2

0
-3

-15

-10

-5

0

-4
-3

5

(b)

-2.5

-2

-1.5

-1

-.5

0

.5

FIG. 3.ÈTrajectories of mass elements of the cloud core are delineatedjust after the second collapse. The innermost and uppermost mass shells
correspond to M \ 10~3 M and M \ 2 ] 10~2 M , respectively. Difr_r_
ference in mass between the neighboring shells is constant in logarithmic
scale, where * log M \ 5.42 ] 10~2. See legend of Table 1 for the deÐnir
tion of t .
0

4

3

2

1

-15

-10

-5

0

FIG. 2.ÈThermal evolution at the center of the cloud core is depicted
for the initially homogeneous model. The upper panel shows the energy
exchange rate between gas and dust (" ), the radiative heating rate (i cE),
gd the compressional heating rate
E of
the radiative cooling rate (4ni B), andP! ) as functions of the central density.
For o Z 10~12 g cm~3 the radiativegheatingand cooling rates are merged. The clower panel shows the central
temperature as a function of the central density.

for a free-fall time of the initial cloud core, i.e.,
n
t \
ff 2

S

R3
\ 1.77 ] 105 yr .
2GM

(2)

The adiabat recovers its steepness after the second collapse
is halted, but c does not reach 5/3 as expected for monoeff
atomic gas because
the partial degeneracy of electrons
begins to dominates the pressure in the center.
All the components of dust grains evaporate at o D 3] 10~7 g cm~3 and the radiative heating and cooling crates
drop suddenly. The increasing density and temperature,
however, again raise the radiative heating and cooling rates.
The compressional heating rate of gas increases with ountil the collapse is almost halted at o D 1 g cm ~3, wherecc more than 6 orders
the compressional heating rate drops byof magnitude and the quasi-static contraction of the protostar follows instead of the dynamical collapse.
3.1.4. T he Birth of a Protostar

L69 and WN reported rebound of the second core justafter its formation. In our results, indeed, a shock wave
propagates outward as if the second core rebounded (seeFig. 1c). The rebound, however, does not actually occur
because the velocity always remains negative everywhere.

The feature that the infalling gas is decelerated (but notreversed) by an outgoing shock wave is observed moreclearly in Figure 3, which illustrates the trajectories of masselements during the second collapse. Here the o†set of time,t 4 1.7526 ] 105 yr, is taken as the instance when the0secondcollapse begins (Table 1). As expected, t is very
0 propaclose to t deÐned by equation (2). A shock wave
ffgates outwardly as settling the infalling material onto thesecond core.
The very short second collapse phase is followed by themain accretion phase (curves 9È13 in Fig. 4), where thecentral protostar grows in mass by the steady accretionfrom the infalling envelope. The duration time of the mainaccretion phase is characterized by the accretion time, t ,which is deÐned by the time required for a mass elementacc
in
TABLE 1
ELAPSED TIMES FOR THE
INITIALLY HOMOGENEOUS
MODEL

Label

t[t
0
(yr) a1 .......
2 .......
3 .......
4 .......
5 .......
6 .......
7 .......
8 .......
9 .......
10 . . . . . .
11 . . . . . .
12 . . . . . .
13 . . . . . .

[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

a The o†set of time, t 4
1.7526 ] 105 yr, represents 0theinstance when the second collapse begins.

No. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.

6

12

13

13

0

(a)

355

(b)

9
10

5

-5

9
4
log T [K]

9
10

-10
3

9

2

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)
11

-1
200

10
-2
9

150
-3
13

100

-4
11
50
-5

10
9
0

-6
-4

-2

0
log R [AU]

2

4

-4

-2

0
log R [AU]FIG. 4.ÈNumerical results for the initially homogeneous model in the main accretion phase are shown. The preceding evolution is depicted in Fig. 1.
Elapsed times are tabulated in Table 1.

the infalling envelope to be settled on the central star :
R
t ^ ,
(3)
acc v
ff
where R is the radius of the infalling envelope, and v 4ff
J2GM/R is the free-fall velocity. Rearranging equation (3),
we obtain
t

acc

^

S

R3
,
2GM

(4)

which gives essentially the same value as the free-fall time(eq. [2]) except for numerical factors of order of unity.
Figure 4 shows the evolution succeeding to Figure 1. The
second core, which is a ““ seed ÏÏ of a protostar, grows monotonically in mass as accreting the infalling material (Fig. 4d).

The radius of the protostar, which is indicated by the““ wall ÏÏ in the density distribution (Fig. 4b), is settled at
2 ] 10~2 AU or 4 R during the main accretion phase.
_This value of the protostellar
radius is in agreement with
WN and SST.
3.1.5. T he T emperature ProÐle in the Protostar

The internal structure of the protostar shows an o†centered peak in the temperature distribution, as also foundby WN and SST. The o†-centered maximum of temperaturearises because of partially degenerate electrons in thecentral 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 enough and356

MASUNAGA & INUTSUKA

the thermal energy dominates the Fermi energy, the gaspressure is more sensitive to temperature and the o†centered temperature becomes higher than in the center.

f=1Vol. 531µ
r3.1.6. T he Infalling Envelope

The density structure in the infalling envelope is approximately o P r~1.5 in the inner region and o P r~2 in theouter, and the boundary between these regions shifts outwardly with time. This feature resembles the singular isothermal self-similar solutions, e.g., the solutions obtained by
Shu (1977) along with expansion-wave solution and thoseby Hunter (1977) which include the continuous counterpartof the LP solution.
Plateau in the mass distribution (Fig. 4d) indicates thecurrent mass of the central protostar. The protostellar massmonotonically reaches 0.73 M at the end of the calcu_
lation.

1/3 < f < 1µr

3.1.7. T he Radiation Field

The upper panel in Figure 5 depicts spatial distributionsof the Eddington factor in the main accretion phase. The
Eddington factor is kept at 1 in the protostellar interior3(R [ 10~2AU), where the opticaldepth is quite large(q D 1015). At the accretion shock front on the protostellar
surface the Eddington factor sharply drops below 1 , which
3 layer,
occurs because the optically thin radiative relaxation
where the shocked gas cools down by radiating away its
thermal energy, is Ðlled with radiation hotter than in thepreshock gas and the postrelaxation layer (Fig. 6 ; see legend

f = 1/3µr1

Eddington factor.8.6

0 < f < 1/3

13

µ

11 12

r

.4
10

-4

-2

0
log R [AU]

2

4

f=0

-4

µ
r

10
11

-5

12

13
-6

-4

-2

0

2

4

FIG. 5.ÈUpper panel depicts the Eddington factor as a function ofradius in the main accretion phase. Elapsed times are tabulated inTable 1. The lower panel shows the mass accretion rate, i.e., M0 (r) 44nr2o(r)[ [ v(r)].

FIG. 6.ÈSchematic picture for the Eddington factor, f 4/1 Ik2dk//1 Idk, is illustrated. Arrows represent the angular distribu~1 of the speciÐc
~1
tion
intensity, I. A special case where the Eddington factor
drops below 1 as well is shown. For the isotropic radiation Ðeld3
(I \ constant), f \ 1 (middle). As the radiation Ðeld becomes outwardly
3
peaking (e.g., the central
luminous source with a di†use and extended
radiation Ðeld) f increases toward unity, which is the extreme case where
only radial component exists [I P d(k [ 1)] (top). In contrast, f descends
below 1 where the radial component of the radiation Ðeld is dominated by3
the orthogonalcomponent, and f falls onto zero in another extreme casewhere I P d(k) (bottom). This case occurs near the strong radiation shock
front, where the narrow radiation relaxation layer Ðlled with hot radiation
is sandwiched between the adjacent cool preshock and postrelaxation
layers.

No. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.

for more details ; see also Balluch 1988). The Eddington
factor rapidly rises toward 1 since the Planck-mean opticaldepth is smaller than unity between the protostellar surface
and the dust sublimation front, at which the dust grains
vapor away, but it decreases again onto 1 near the dust
3
sublimation front because the opaque dusty envelope reprocesses the outwardly peaking radiation from the center intoa 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 and theoutwardly peaking radiation dominates the isotropic Ðeld.
At the Ðnal stage (curve 13) the Eddington factor departs
from 1 even near the dust sublimation front. This is because3
the entire
dusty envelope turns to be optically thin due tothe depletion of the infalling material.
The lower panel in Figure 5 shows the distribution of themass accretion rate, M0 . On the protostellar surface M0 (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~2 AU) is D10~5
M yr~1 independently of radius except near the outer_ in the early main accretion phase (curves 10È12). Thisedge
value of M0 is consistent with SST, who adopted M0 4 10~5
M yr~1 D c3 /G as a free parameter, where c is the e†ec_ sound velocitys0
tive
in the initial cloud core, s0
However, M0
decreases to D10~6 M yr~1 (curve 13) at later stages in_ Our calculation supposes no incithe main accretion phase.
dent mass Ñux through the outer boundary and therefore
the mass accretion rate should decline after a signiÐcant
fraction of the total mass is accreted by the protostar.
The mass accretion rate shows a somewhat di†erent
behavior in the initially hydrostatic model, as seen in ° 3.2
below.
3.2. L uminosity Curves
While the luminosity remains [0.01 L during the Ðrstcollapse as shown in Paper I, it rapidly_ rises when thesecond core, or a protostar, is born. The luminosity reaches
the maximum value of B25 L and then decreases grad_ model (Fig. 7). The lumiually in the initially homogeneous
nosity of a protostar surrounded by an infalling envelope isdominated by the accretion luminosity :P

GM (t)M0 (t) GM0 (t) t*** \M0 (t@) dt@ ,
(5)
*
acc
R (t)
R (t)
t0
*
*
where M is the protostellar mass, R is the protostellar* M0 is the mass accretion rate
* onto the protoradius, and*star. As the protostellar radius is settled at 4 R during themain accretion phase (°3.1), the shape of the_ luminosity
curve reÑects the time evolution of M and M0 , which* rise in the* lumiclosely relate with each other. The initial
nosity 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 7. The timeaveraged mass accretion rate reads 5.7 ] 10~6 M yr~1
_ M -t(0 \ t [ t \ 1.4 ] 105 yr) from the gradient of the0curve, while M0
just after the birth of a protostar * is*
1.3 ] 10~5 M yr~1.
_
A homogeneouscloud core begins the gravitational collapse with inward acceleration because it has no pressure
gradient and thus the force balance is not satisÐed initiallyagainst the self-gravity. On the other hand, collapse of a
L

(t) \

357

1.5

1

.5

0

-.5

-1
0

50000

0

50000

1

.8

.6

.4

.2

0

FIG. 7.ÈLuminosity curve (upper panel) and the protostellar mass
(lower panel) are shown as functions of time after a protostar is formed. t is0deÐned in Table 1. Solid and dotted lines denote the initially homogeneousmodel and the initially hydrostatic model, respectively.

hydrostatic cloud core is initiated without signiÐcant acceleration, which yields a lower mass accretion rate than in theinitially homogeneous case. Although the dynamical evolution ““ forgets ÏÏ the initial density proÐle in the soundcrossing time (Paper I), small surviving inÑuence on themass accretion rate causes a qualitative di†erence in the
behavior of the luminosity curve as seen in Figure 7.
In contrast to the initially homogeneous model, the initially hydrostatic model produces neither a rapid rise fort Z t nor a decrease after the maximum in the luminosity
0 and instead shows a gradual and monotonic increase
curve,of the luminosity. The M -t diagram presents a linear*
growth of the protostellar mass
for this model, which indicates that 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 massaccretion rate is5.3 ] 10~6 M yr~1, which is very close to
the time-averaged value in _the initially homogeneousmodel. The self-similar collapse solutions yield the massaccretion rate of a factor times c3/G independently of theinitial density conÐguration and ofs the total mass of a cloud
core, which explains the reason why the time-averaged massaccretion rates in the two models accord with each other.
Nevertheless, the mass accretion rate as a function of time
exhibits di†erent behaviors depending on the initial condition and the luminosity curve sensitively reÑects the di†erence in M0 between the two models.
The di†erence in the total mass also inÑuences the behavior of the luminosity curve. A less massive hydrostatic core

358

MASUNAGA & INUTSUKA

should be more centrally condensed so as to be gravitationally unstable, and hence the free-fall time (P1/JGo) is
shorter. In contrast, the condition when a protostar is born
depends on neither the initial conÐguration nor the total
mass of a cloud core because the condition of the violationof isothermality depends only on the initial temperatureand opacity (Masunaga & Inutsuka 1999a). This meansthat a less massive cloud core is not allowed to spend sufficient time before the birth of a protostar to raise the massaccretion rate as large as that for a more massive cloud core.
In other words, the Ðrst collapse phase of a less massive coreis ended at an earlier stage when the convergence to the LP
solution is less complete. As a consequence, any hydrostatic
cloud cores would yield only a gradual rise in the luminosity curve after the birth of a protostar unless the total
mass is substantially larger than 3.852 M . A less massive_hydrostatic core, say a 1 M core, should therefore possess_
the similar properties in the luminosity curve to the 3.852
M case, and the di†erence in the total mass is not_
responsible
for the di†erent behavior in the luminosity
curves between our initially homogeneous and hydrostatic
models.
SST found that L increases monotonically with time. Evidently this increase of L shown by SST originates their““ steady state approximation ÏÏ in which M0 is Ðxed at 10~5
* L
M yr~1 throughout the calculation, where
inevitably_
rises as the central protostar grows in mass.accThe steadystate approximation by SST is valid for our initially hydrostatic model but is invalid for the initially homogeneousmodel. WN adopted a homogeneous cloud core as theinitial condition, and thus they found the same tendency as
in our initially homogeneous model on the contrary to SST.
Even for the initially homogeneous model where theluminosity declines with time, the Ðnal value of the luminosity (D10 L ) is still too high for the typical value_
for class I sources(D1 L , Kenyon et al. 1990). This““ luminosity problem ÏÏ for _class I protostars, which waspointed out by Kenyon et al. (1990), indeed arose in ourspherically symmetric calculation. The luminosity problem,
however, can be resolved by considering a possible di†erence in the mass accretion rates between the circumstellar
disk and the infalling envelope (Kenyon, Calvet, & Hartmann 1993), which is not taken into account in our spherically 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 which halts themass accretion onto a protostar would include complicatedprocesses for actual YSOs. Nakano, Hasegawa, & Norman(1995) showed that the mass outÑows from YSOs Ðnallyhalt further mass accretion by blowing o† the residual
envelope matter until the envelope becomes no longer
gravitationally bound. This e†ect would suppress M0 at* tolate stages of the evolution, which may also contributeresolve the luminosity problem.
Which of the decreasing or increasing luminosity accounts for observations ? To discuss observational counterparts of our results, we present results of SEDs in the next
section.
3.3. Spectral Energy Distributions
The computational result for the emergent SED in theinitially homogeneous model is shown in Figure 8. The
SEDs are calculated from the monochromatic outgoing

Vol. 531

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. 8.ÈEvolution of the spectral energy distribution for the initiallyhomogeneous model is shown. Elapsed time of each curve is listed in Table
1. The upper panel delineates the emergent Ñux, F , in the unit of Jy, whilel
the lower panel depicts lF in the cgs unit. An observeris assumed to be ata distance of 150 pc with a ltelescope of the resolution of 10A.

radiation through the cloud edge, convolved by the Gaussian beam pattern with a given FWHM. In Figure 8 an
imaginary telescope is assumed to observe the object at a
distance of 150 pc with the resolution of 10A (FWHM). Theinterstellar extinction is not considered.
At the initial stage the SED reÑects optically thin thermal
emissions by dust grains of nearly 10 K (the curve denotedby ““ 1 ÏÏ). The SED once becomes colder (curves 2 and 3),
due to a slight decrease of temperature mentioned in ° 3.1,
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 theouter envelope as found also in Paper I.
A deep absorption feature at 1013.5 Hz (or 10 km) corresponds to a striking feature in silicate. The absorption isextremely 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 km 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 theprotostar.
The object would be detectable in near-infrared wavelengths at the last stages (curves 12 and 13). The SED
labeled ““ 12 ÏÏ shows typical features for class I protostars : a

No. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.

large excess of the far-infrared component and undetectability in optical light. The SED labeled ““ 13 ÏÏ corresponds to a more evolved YSO such as a class II object, forwhich the enhancement of luminosity in the far-infrared has
disappeared and the peak in lF has been shifted to thelnear-infrared. SED 13 shows a Ñat spectrum in the infraredand hence this object might correspond to a Ñat spectrum T
Tauri star if the optical component appears in the SED
owing to nonspherically-symmetric e†ects such as bipolar
cavities (° 4.2). We stopped the calculation at this stage
because the further evolution of SED should involve a nonspherically symmetric treatment, where emissions from thecircumstellar accretion disk and photons that escapethrough the bipolar cavities a†ect the SED in the nearinfrared and the optical.
Figure 9 shows the SED in the initially hydrostatic
model. One Ðnds no essential di†erence from the initiallyhomogeneous model (Fig. 8) except that the Ñux in the
millimeter and submillimeter (l \ 12.5 Hz) never decreases
even at the Ðnal stage of the evolution. This is because ourinitially hydrostatic model contains a larger mass of thecloud core (° 2.2) than the initially homogeneous model,
which keeps providing the thermal emission of the infallingenvelope in millimeter and submillimeter wavelengths.
However, the spectral appearance in shorter wavelengths isquite insensitive to the initial density proÐle and, more
remarkably, to the total mass (or the radius) of the initial

10
8
5,6

0

3
-2

1

-4
11

11.5

12

12.5

13

13.5

TABLE 2
ELAPSED TIMES FOR THE
INITIALLY HYDROSTATIC
MODEL

Label

t[t
0,e
(yr)a

1 .......
2 .......
3 .......
4 .......
5 .......
6 .......
8 .......
9 .......
10 . . . . . .

[2.4674 ] 105[3.6097 ] 104[7.6469 ] 102[8.7500 ] 10~1
[1.7188 ] 101
1.7656 ] 100
2.1568 ] 104
6.5776 ] 104
1.7614 ] 105

a The o†set of time, t 42.4674 ] 105 yr, represents0,etheinstance when the second collapse begins.

cloud core. The reason for this is that the optical depth,q \ / io dr, is insensitive to the outer boundary when thedensity decreases with radius more rapidly than o P r~1,
which is the case in protostellar envelopes.
4. DISCUSSION

In this section we discuss on possible interpretations ofstatistical trends in observations on the basis of our computational results. In ° 4.1 we examine the bolometric temperature, which is a quantitative expression of SEDs, torelate our results with observed objects. Considering expected nonspherically symmetric e†ects, we illustrate an
evolutionary picture of protostar formation in ° 4.2.

4

2

359

14

14.5

15

4.1. Bolometric T emperature
As a quantitative measure of evolutionary status of
YSOs, Myers & Ladd (1993) proposed the bolometric temperature, T , deÐned bybol/= lB (T )dl0l bol \ l6 ,(6)/= B (T )dl
0 l bol
where-8/= lF dll
l6 4 0(7)/= F dl
0 lon the basis of an observed SED, F . Combining equationsl(6) and (7), one obtains-10f(4) hl6
T 4\ 1.25 ] 10~11l6 Hz ~1K ,bol 4f(5) k

-6

(8)

where f(n) is the Riemann zeta function of argument n :
-12

-14
11

11.5

12

12.5

13

13.5

14

14.5

15

FIG. 9.ÈEvolution of the spectral energy distribution for the initiallyhydrostatic model is shown. Elapsed time of each curve is listed in Table 2.
The upper panel delineates the emergent Ñux, F , in units of Jy, while thellower depicts lF in cgs units. An observer is assumed
to be at a distance ofl
150 pc with a telescopeof the resolution of 10A.

P

= tn~1dt .
(9)et [ 10
Myers & Ladd (1993) showed that T
successfully disbol from protostars
tinguishes evolutionary positions of YSOs
to weak-line T Tauri stars. Gregersen et al. (1997) foundthat the mean value of T is 44 K for 16 class 0 sources.
bol
To understand the statisticaltendency in T of observedbol
protostars, we calculate the bolometric temperature
from
our SED (Fig. 10). The bolometric temperature monotonf(n) 4

1
!(n)

360

MASUNAGA & INUTSUKA

100

10

1
-50000

0

50000

100

Vol. 531

core is (marginally) Jeans unstable. A cloud core that is notin hydrostatic equilibrium obtains a relatively large M0*
because of the initial inward acceleration, which leads tolarge luminosity just after the birth of a protostar (° 3.2). In
contrast, the initially hydrostatic model, where a lower M0*admits neither the rapid growth of a protostar nor theimmediate exhaustion of the infalling material, fails to
reproduce the qualitative tendency in the luminosity ofobserved YSOs and instead the luminosity increases gradually and monotonically. This implies that the observed
evolution of luminosity does not favor the spontaneous
mechanism of star formation (° 2.2).
Nevertheless, we do not exclude the dissipation of turbulence (Nakano 1998) as a possible mechanism for initiatingstar formation in case that the ““ macroturbulence,ÏÏ where
the scale length of turbulent Ñuctuations is comparable withthe scale of the system and hence turbulent motions can nolonger be treated as an e†ective sonic velocity, dominates
the turbulent Ðeld because dissipation of macroturbulence
would not necessarily produce a nearly hydrostatic cloudcore when the dynamical collapse is initiated.

10

4.1.2. Comments on the Generality of the Initial Condition

1-50000050000

FIG. 10.ÈEvolution of the bolometric temperature, T , is shown. Thebol for l6 given bysolid line delineates T estimated by the original deÐnition
boldenote T estimated using l6 with a upper cuto†eq. (7). The other lines
bol where f is 3 ] 1012 Hz for the dottedfrequency, f , in the integral (see text),
c Hz for short-dashed line, andc 1 ] 1013 Hz for the longline, 5 ] 1012
dashed line. The upper panel is for the initially homogeneous model, andthe lower panel isfor the initially hydrostatic model.

ically increases with time. While T remains at 10 K before
a protostar is formed (t [ t \ 0),bol
it abruptly rises up to 40
K after t [ t exceeds 0 and0 then gradually increases to 410
K. Chen et 0al. (1995) found that observed class 0 sources
show T less than 70 K, which corresponds to t [ t ¹ 2bolin our results. This is smaller by an order of 0mag] 104 yr
nitude than the typical age of class I sources. We call thesevery young protostars as ““ genuine ÏÏ class 0 objects. We Ðndno essential di†erences in the evolution of T between the
bol hydrostatic
initially homogeneous model and the initiallymodel as expected from the discussion in ° 3.3, in contrast to
the luminosity curve (° 3.2).
4.1.1. Observational Implication to the Initial Condition

The luminosity curve and T
in the computational
results implies that genuine class bol
0 sources should be moreluminous than class I sources for the initially homogeneousmodel, and less luminous for the initially hydrostatic model.
Observations have revealed that the typical luminosities are
D1 L for class I sources (Kenyon et al. 1990) and 5.7 L_ 0 sources (Gregersen et al. 1997). Obviously protos_
for classtars 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, andany other initial density structures that are not in hydrostatic equilibrium would trace a dynamical evolution
similar to the initially homogeneous model rather than to
the initially hydrostatic model, as long as the initial cloud

An initially hydrostatic cloud core of 1 M , instead of_ luminosity
3.852 M of our choice, would yield a declining_
curve due to the exhaustion of the infalling material as in
the initially homogeneous case. This model, however, does
not match the observations because it cannot produce sufficient luminosity at 2 ] 104 yr, i.e., the genuine class 0 stage,to account for observed class 0 sources. The mean luminosity of observed class 0 sources is 5.7 L , but the genuine
class 0 objects should be signiÐcantly_ more luminous
because the contamination of heavily obscured (or class0Èlike class I) sources should have lowered the mean luminosity of observed class 0. Furthermore, we can present asfollows another evidence for our point that initially hydrostatic cloud cores do not reproduce the observational tendency. The decline of the luminosity curve is observed when
the protostellar mass arrives at D0.5 M in the initiallyhomogeneous model of the total mass of 1_M (Fig. 7). On_
the other hand, it takes 0.8 ] 105 yr for a protostar
to grow
in mass to 0.5 M in the initially hydrostatic model, while it_ the initially homogeneous model. Thisis much earlier for
implies that the luminosity would turn to decrease att [ t \ 0.8 ] 105yr for an initially hydrostatic core of
1 M 0 , or actually more later because a less massive hydrostatic_ core would provide a smaller mass accretion rate after
a protostar is born as discussed in ° 3.2. On the other hand,the class I stage has already begun before that time (see Fig.
11). Therefore the infalling material is exhausted too late to
explain the observational trend of the decreasing luminosity
from class 0 to class I for an initially hydrostatic cloud coreeven if the total mass is limited to 1 M .
It may be helpful to consider also_ a hydrostatic initialcondition with a highly centrally concentrated proÐle,which corresponds to spontaneous star formation with the
ambipolar di†usion of magnetic Ðelds. In this case, the mainaccretion phase begins as the outer envelope remains inhydrostatic equilibrium since the free-fall time in the central
region is quite short (° 3.2) compared to the accretion timeof the whole envelope. The obtained mass accretion rate,
therefore, would be still much smaller and thus the luminosity curve would monotonically increase more graduallythan in the marginally unstable hydrostatic case. ThisNo. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.

t - t 0[105 yr]

361

Inclinationedge-onpole-on-2-1

M[M * ]
0.25

Starless Cores

0
0.2

Genuine Class 0Class 0-like

0.5

Molecular
Outflows(?)

Class I
0.7

1

0.8

Class I

Flat Spectrum
T Tauri Stars
2

FIG. 11.ÈSchematic picture of the evolution of protostar formation is illustrated. The vertical axis represents the evolutionary time in the unit of 105 yr,
where the o†set value, t , indicates the instance of the birth of a protostar (see Table 1). The protostellar mass, M , is shown based on Fig. 7. The epoch in* YSO to an observer. The area occupied
which outÑows associate0 with YSOs is also indicated. The horizontal axis corresponds to the inclination angle of theby each type of YSOs roughly corresponds to the number of the objects expected to be observed.

picture mimics the expansion-wave solution by Shu (1977),
which has the minimum mass accretion rate, 0.975c3/G atsthe origin, among the family of the inside-out solutionsofthe isothermal self-similar collapse.
4.1.3. T he E†ect of the Circumstellar Disk

Our calculations suppose spherical symmetry. Actually,each YSO can 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 high column density, so thatthe object would be seen as if it were less evolved. To
examine how much this e†ect modiÐes T , we also show in
bol would be estiFigure 10 the bolometric temperatures whichmated if emissions were deÐcient in high frequencies. We
calculate the mean frequency, l6 , imposing an upper cuto†
frequency, f , that is,c/fc lF dll ,(10)
l6 4 0/fc F dl
0 lwhich is substituted for equation (7). Figure 10 indicatesthat Tis sensitive to f for 3 ] 1012 Hz \ f \ 1
c] 1013bolHz and T
can bec signiÐcantly underestimatedbol
even after t D 105 when near infrared components are
hidden by locally enhanced column density 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 compared togenuine class 0 sources, which are less abundant by an orderof magnitude than class I sources. Hence the contaminationof edge-on class I objects into class 0 sources is not negligible.
A clue for discriminating genuine class 0 objects fromedge-on class I objects lies in the luminosity of individual
sources since our initially homogeneous model shows thatthe luminosity reaches the maximum value of 25 L at_
2 ] 104 yr, which corresponds to the typical age of genuine
class 0, and the luminosity declines as the evolution proceeds. YSOs with cold T that have luminosities considerbol be genuine class 0 objects. On
ably larger than 1 L would_
the other hand, cold YSOs with comparable luminosities
with class I sources are probably more evolved protostars
than genuine class 0 objects. While the luminosity problem
for class I objects can be resolved by a possible di†erence inthe mass accretion rates between the circumstellar disk andthe infalling envelope (° 3.2), the luminosity problem would
not arise for the protostars substantially younger than class
I, where nonspherically symmetric e†ects are less importantbecause a signiÐcant fraction of the angular momentum of
the rotating cloud core still remains far away from thecentral protostar. Hence the large luminosity for genuine
class 0 objects in our results, i.e., D25 L , would not seriously su†er from the luminosity problem._
Unfortunately, the bolometric luminosity estimated fromobservations su†ers from a large uncertainty for such young
protostars. A major fraction of the luminosity is emitted in
the submillimeter and only an upper-limit value is obtained
for the infrared Ñux. Furthermore, the luminosity of a YSO

362

MASUNAGA & INUTSUKA

on an edge-on view would be underestimated because onlya small fraction of the total emission leaks out through thedisk plane. Direct comparison of the observed luminosities
with our numerical results require a great care.
4.2. An Evolutionary Picture of Protostar Formation
In the preceding section we established a theoreticalmodel for protostar formation which accounts for recent
observations. The theoretical treatment in this model isfully consistent in the framework of 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†ects due
to rotation and magnetic Ðelds. In this subsection we
discuss possible contributions of these nonsphericallysymmetric e†ects to our results.
We pointed out in ° 3 that edge-on class I objects would
be observed as class 0 sources (““ class 0Èlike class I ÏÏ), andthey would contaminate the ““ genuine ÏÏ class 0 objects
which have the age of 2 ] 104 yr and are substantially
younger than class I sources. Similarly, evolved class I
objects would be seen as Ñat spectrum T Tauri stars
depending on the inclination angle to an observer.
Theoretical and observational evidence supports the ideathat Ñat spectrum T Tauri stars are associated with circumstellar envelopes. Natta (1993) reproduced the observed Ñat
spectra by considering the back-warming of the disk by
scattered stellar photons in the envelope. Calvet et al. (1994)
showed that the thermal emission from the extended
envelope with a polar cavity forms the Ñat spectra in addition to the emission from the central star and disk. Kikuchi,
Nakamoto, & Ogochi (2000) performed two-dimensionalradiative transfer calculations and successfully produced
the Ñat spectra by the back-warming of the disk by thescattering and thermal reprocessing of stellar photons in
the envelope. Observations have revealed the presence ofextended envelopes surrounding HL Tau (Hayashi, Ohashi,& Miyama 1993) and T Tauri (Momose et al. 1996), which
are typical Ñat-spectrum T Tauri stars. T Tauri stars withÑat spectra tend to have larger extinction than typical T
Tauri stars (Osterloh & Beckwith 1995). A recent highresolution 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.
Compiling comprehensively our Ðndings and other theoretical and observational evidence, we illustrate an evolutionary picture of protostar formation in Figure 11.
Protostellar evolution proceeds from top to bottom. Beforet reaches t , where t D 2 ] 105 yr indicates the instance of
0 a protostar
0 (see Table 1 in ° 3.1), a cloud core
the birth ofexperiences gravitational collapse without any signiÐcant
changes in SEDs (° 3.3). Such cloud cores would beobserved as starless cores even after the formation of theÐrst core (Paper I). After a protostar is born (t [ t ), the0appearance of a YSO depends not only upon the evolutionary time but on the inclination angle to an observer, whichis indicated 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 view su†ers fromenhanced column density in the disk plane, which causes

Vol. 531

signiÐcant underestimation of the bolometric temperature(° 3.3). A YSO at a certain evolutionary position is therefore
observed as if it were more evolved on a pole-on view andless evolved on an edge-on view. This tendency is schematically illustrated in Figure 11. Observed class 0 sources
consist of genuine class 0 objects and edge-on class I, orclass 0Èlike class I, objects. A Ñat spectrum T Tauri star
would be evolved to a typical classical T Tauri star as thecircumstellar envelope is dispersed or exhausted.
The area occupied by each type of YSOs in Figure 11
roughly corresponds to the number of the objects expected
to be observed. On the basis of this fact, one can give a
constraint on the epoch in which molecular outÑows associate with YSOs by considering a statistical trend in observations. Mardones et al. (1997) listed 23 class 0 sources and
24 class I, where 23 sources of class 0 (100%) and 18 sourcesof class I (75%) show evidence of outÑows. Therefore, theoutÑow epoch should cover the whole region of genuine
class 0 and class 0Èlike class I, and three-fourths of class I in
Figure 11. If an outÑow begins when the main accretion
phase starts (i.e., t \ t ) as predicted by Tomisaka (1998),
0 outÑows should continue until
Figure 11 shows thatt [ t D 1.3 ] 105 yr to naturally reproduce the obser0 This estimation is, of course, very crude because it
vations.
assumes that outÑows are driven universally and steadily inthe same epoch regardless of the individuality of eachobject. More elaborate estimation requires a further progress in both observational and theoretical studies in thedriving 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 class 0Èlike class I
objects should have large inclination angles.
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(25 L , see ° 3.2) and have a small protostellar mass
(0.1 M_ ). Momose et al. (1998) derived the protostellar
mass as_D0.1 M for L1551 IRS5, whose luminosity is D30
L (Cohen et _al. 1984). Hence L1551 IRS5 may be a_genuineclass 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 it appears as a class I object
despite its extreme youth.
A more quantitative picture than Figure 11 should
involve numerical calculations based on multidimensionalradiation hydrodynamics. Such a calculation requires verylarge computational costs and is not yet available atpresent, but it will be a challenging task in future.
5. SUMMARY

Our Ðndings in this paper are summarized as follow.
1. RHD numerical calculations are carried out to study
physical processes in the formation of a 1 M protostar.
_
The whole evolution are pursued from the beginningof theÐrst collapse to the end of the main accretion phase. All thephysical phenomena that appear in protostellar formation,(e.g., the dynamical collapse of a cloud core, the formationof the pressure-supported core, the accretion shock front on
the protostellar surface, and the supersonically infalling gas)
are incorporated self-consistently in the calculations.

No. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.

2. A protostar is born after the initial isothermal collapseof a cloud core (the Ðrst collapse) and the subsequent collapse triggered by dissociation of molecular hydrogen (thesecond collapse). The protostar grows in mass as accreting
the infalling material from the circumstellar envelope while
the protostar keeps its radius at D4 R during the main_accretion phase. These typical features in the evolution are
in good agreement with previous studies.
3. For the initially homogeneous model, the accretion
luminosity rapidly rises to the maximum value of 25 L 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 initiallyhydrostatic model. This discrepancy arises because the massaccretion rate varies depending inward acceleration at theinitial stage, which a†ects the luminosity curve. Any hydrostatic cloud cores would yield only a gradual rise in theluminosity curve after the birth of a protostar unless thetotal mass is substantially larger than 3.852 M .
4. We conÐrm that the SED evolves from _a 10 K graybody spectrum to hotter spectra typical for class I and Ñat
spectrum sources. The bolometric temperature derived by
the obtained SED shows that cold sources with T lessthan 70 K are younger than 2 ] 104 yr (““ genuine ÏÏ bolclass 0objects). This indicates that class 0 objects have younger
ages by an order of magnitude than class I objects.
However, T can be signiÐcantly lower even after t D 105
bol
yr if the infallingmaterial with angular momentum falls
onto a circumstellar disk and near-infrared components are
hidden by locally enhanced column density in the disk
plane. Observed class 0 sources should be the compound of
the genuine class 0, which is as young as 104 yr, and more
evolved protostars on edge-on view (““ class 0Èlike class I ÏÏ).
The contamination of edge-on class I objects are not negligible because they are intrinsically abundant than genuine
class 0 objects.
5. Since observations indicate that class 0 sources are
typically more luminous than class I sources, the initiallyhydrostatic model, in which the luminosity increases monotonically with time, does not explain the observation. The
initially homogeneous model, in contrast, shows the ten-363

dencies consistent with the observations. This implies thatthe observed evolution of luminosity does not favor the
spontaneous mechanism of star formation. The decreasing
luminosity with time implies that signiÐcantly luminousand cold YSOs (L ?1 L and T \ 70 K) should be_bolgenuine class 0 objects. Less luminous YSOs are more
evolved protostars on edge-on view even if they have coldT .
bol
6. Finally, we propose an evolutionary picture of protostar formation compiling our Ðndings in this paper and othertheoretical and observational evidence. In terms of the evolutionary time and the inclination to an observer, protostellar objects are found to be clearly categorized.
In our numerical code, the convective energy transport
and the nuclear energy production are not taken into
account. The deuterium burning, which is ignited at
T D 106K, should occur in the protostar at the Ðnal stage ofour 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 considered both the convection anddeuterium burning and they found the protostellar radius of3 D 4 R , which accords with our result quite well. Fur_ a slight di†erence in the protostellar structure
thermore,would give no inÑuence on the structure of the infallingenvelope. Hence the SEDs would remain unchanged byincorporating the convection and the nuclear reactions.
The authors are grateful to Shoken M. Miyama, Takenori Nakano, and Munetake Momose for their helpful comments and discussions. H. M. thanks Harold W. Yorke forhis kindly sending me tabulated data of the dust opacities.
Helpful comments and questions of the referee, Dean
McLaughlin, with which some ambiguous points in theoriginal version have become more clear and convincing.
The research of H. M. is supported by the Research Fellowships of the Japan Society for the Promotion of Science for
Young Scientists.

APPENDIX A
BASIC EQUATIONS
In spherical symmetrical Lagrangian (comoving) coordinates, the basic equations of radiation hydrodynamics are written
as follows :
dr
\v ,
dt

(A1)

3 LM
r,
4n L(r3)

(A2)

o\

Lp
igas
v
Dv
GM
r [ 4nr2] F F[ Q ,
\[
o
r2
LM
c
Dt
r
E
D 1
v
L(4nr2F)De]] ( fE ] p)
[ (3f [ 1)E\v ]v [" [
,
CR
Q
gd
o
Dt o
or
LM
Dt
r
1 DF 4nr2o L( fqE)
s F]\[ F ,
c2 Dt
q
LM
c
r

A B

AB

(A3)
(A4)
(A5)

364

MASUNAGA & INUTSUKA

Vol. 531

AB

De
D 1
\ (cigasE [ 4nigasB) ] v ] v [ " .
]p
E
P
CR
Q
gd
Dt
Dt o

(A6)

Here v and v are taken from the tensor artiÐcial viscosity in the formulation by Tscharnuter & Winkler (1979),
Q
Q
L(r2v) Lv L(r2v)
9 L
r3l2o
[
,
v \
Q r Lr3
Lr3 Lr
Lr3

G

C

C

DH

D

L(r2v) Lv L(r2v) 2[,
v \ [4.5l2
Q
Lr3
Lr3 Lrwhere l, a coefficient with dimension of length, is of the order of the local width of grid to smear out the shock front artiÐciallyso that the discretized hydrodynamic equations can handle shock waves. Both v and v are switched to zero where +Év [ 0,QQthat is, the artiÐcial viscosity works only when the gas is compressed.
The Eddington factor, f, and the frequency-averaged quantities are calculated based on solutions of the angle- andfrequency-dependent radiative transfer equation according to the following equations :K
f \ l,l Jl/= i H dl
i \ 0 l l ,
F
H

(A7)
(A8)

/= i J dl
i \ 0 l l ,
E
J

(A9)

/= i B dl
i \ 0 l l ,P
pT 4/n(A10)/= f J dl
f\ 0 l l ,
J

(A11)

and

where J ,H , and K are the Ðrst three moments of the monochromatic speciÐc intensity, and those without a subscript l are
l integratedll over frequency. The sphericality factor q is deÐned by (Auer 1971)
the valuesP

r 3f [ 1
dr@ .
fr@
Rin
Heating or cooling by collisions between gas and dust is taken into account through the energy exchange rate, " :gd" \ ; "(n) ,
gd
gd
(n)
where the superscript (n) denotes the dust components. We consider three components of dust grains, i.e., amorphous carbon,
silicate, and dirty ice-coated silicate, on the basis of the dust model by Preibisch et al. (1993). The energy exchange rate foreach dust component is expressed as (Hayashi & Nakano 1965 ; Leung 1975)
ln q \A B

3 8k 1@2 n
H2 n(n)[na(n)2]T 1@2k[T [ T (n)] ,(A12)"(n) \
gr
gd 2 nmo gr
H2
where n(n) and T (n) are the number density and temperature of dust grains of component (n). The dust temperature is
gr by (Hayashi
grdetermined& Nakano 1965),"(n) ] i(n),abs cE [ 4i(n),abs pT (n) 4 \ 0 .
(A13)
gd
E
P
gr
The nuclear energy productions and the convective energy transport are not considered.
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 is achieved to a desirable extent. One thus obtains the radiationhydrodynamic quantities which are consistent with the radiation Ðeld at each time step. While this procedure works wellduring calculations for the Ðrst collapse (Paper I), convergence is found to be much worse after the second collapse due tocomplicated behaviors of the Eddington factors near the accretion shock front. Therefore convergence between the transfer
equation solver and the 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 the variations of the Eddington factors should not be
too large.
In Eulerian coordinates, the basic equations (eqs. [A1]È[A6]) are rewritten intoLMLMr ] 3vr2
r\0 ,
Lt
Lr3

(A14)

No. 1, 2000

RHD MODEL FOR PROTOSTELLAR COLLAPSE. II.
o\

3 LM
r,
4n Lr3

365
(A15)

Lv
Lv
GM
Lp igasv
r[ ] F F[ Q ,
]v \[
Lr
c
o
Lt
Lr
r2

(A16)

E
L
E
3 L(r2v)v
L(r2F)Le]]ve]] (p ] fE)[ (3f [ 1)E \ v ] v [ " [ 3
,
CR
Q
gd
o
Lr
o
o Lr3
or
Lr3
Lt

(A17)

LF
1 L( fqE)
s F
1 LF
]v
]
\[ F ,
Lr
q Lr
c
c2 Lt

(A18)

Le
Le
3 L(r2v)]v ]p\ (cigas E [ 4nigas B) ] v ] v [ " .
E
P
CR
Q
gd
Lt
Lr
o Lr3

(A19)

A B

A B

A

B

and

Equations (A12) and (A13) remain unchanged.
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 not useful when a sharp structure such as a shock
wave is generated in di†use areas where grid points are sparsely distributed. By this reason, calculations are performed on
Lagrangian coordinates for the initial isothermal collapse, and the Eulerian takes the place after the adiabatic core isformedÈi.e., after the accretion shock is generated. This method was found to successfully resolve the accretion shock at thesurface of the Ðrst and second cores in our calculations. In early versions of our code, second-order accurate advection (i.e., thevan Leer advection) was adopted for Eulerian calculations, supplementing the Lagrangian RHD equations (A1)È(A6) (Paper
I). This method works well for 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 time steps become
unacceptably small compared to the evolutionary timescale at late stages of protostellar collapse. Therefore the discretized
equations in the Eulerian formulation (eqs. [A14]-[A19]) have to be directly incorporated into the linearized matrix, insteadof employing the external advection terms solver. Since upstream di†erentiation of Ðrst-order accurate is adopted, this ““ direct
Eulerian scheme ÏÏ is inferior in accuracy to the ““ Lagrangian with the Eulerian remapping scheme ÏÏ using the van Leer advection. Nevertheless, it has a great advantage that large CFL numbers are permitted and computational time is reduced to areasonable extent.
REFERENCES
Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122
Mardones, D., Myers, P. C., Tafalla, M., Wilner, D. J., Bachiller, R., &Alexander, D. R., & Ferguson, J. W. 1994, ApJ, 437, 879
Garay, G. 1997, ApJ, 489, 719
Auer, L. H. 1971, J. Quant. Spectrosc. Radiat. Transfer, 11, 573
Masunaga, H. 1999, Ph.D. thesis, Univ. Tokyo
Balluch, M. 1988, A&A, 200, 58
Masunaga, H. & Inutsuka, S. 1999a, ApJ, 510, 822
Barsony, M. 1994, in Clouds, Cores, and Low Mass Stars, ed. D. P. ClemÈÈÈ. 2000, ApJ, 535, in press
ents & R. Barvainis (Provo : ASP), 197
Masunaga, H., Miyama, S. M., & Inutsuka, S. 1998, ApJ, 495, 346 (Paper I)
Bodenheimer, P., & Sweigart, A. 1968, ApJ, 152, 515
Momose, M., Ohashi, N., Kawabe, R., Hayashi, M., & Nakano, T. 1996,
Bonner, W. B. 1956, MNRAS, 116, 351
ApJ, 470, 1001
Boss, A. P., & Yorke, H. W. 1995, ApJ, 439, L55
Momose, M., Ohashi, N., Kawabe, R., Nakano, T., & Hayashi, M. 1998,
Calvet, N., Hartmann, L., Kenyon, S. J., & Whitney, B. A. 1994, ApJ, 434,
ApJ, 504, 314
330
Myers, P. C., & Ladd, E. F. 1993, ApJ, 413, L47
Chen, H., Myers, P. C., Ladd, E. F., Wood, D. O. S. 1995, ApJ, 445, 377
Nakano, T. 1979, PASJ, 31, 697
Close, L. M., Roddier, F., Northcott, M. J., Roddier, C., & Graves, J. E.
ÈÈÈ. 1998, ApJ, 494, 587
1997, ApJ, 478, 766
Nakano, T., Hasegawa, T., & Norman, C. 1995, ApJ, 450, 183
Cohen, M., Harvey, P. M., Schwartz, R. D., & Wilking, B. A. 1984, ApJ,
Natta, A. 1993, ApJ, 412, 761
278, 671
Osterloh, M. & Beckwith, S. V. W. 1995, ApJ, 439, 288
Ebert, R. 1955, Z. Astrophys., 37, 217
Penston, M. V. 1969, MNRAS, 144, 425
Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303
Preibisch, Th., Ossenkopf, V., Yorke, H. W., & Henning, Th. 1993, A&A,
Goldsmith, P. F., & Langer, W.D. 1978, ApJ, 222, 881
279, 577
Gregersen, E. M., Evans, N. J., II, Zhou, S., & Choi, M. 1997, ApJ, 484, 256
Saumon, D., Chablier, G., & Van Horn, H. M. 1995, ApJS, 99, 713
Hayashi, C., & Nakano, T. 1965, Prog. Theor. Phys., 34, 754
Shu, F. H. 1977, ApJ, 214, 488
Hayashi, M., Ohashi, N., & Miyama, S. M. 1993, ApJ, 418, L71
ÈÈÈ. 1995, in Molecular Clouds and Star Formation, ed. C. Yuan &Hirano, N., Kameya, O, Mikami, H., Saito, S., Umemoto, T., & YamaJ. H. You (Singapore :World ScientiÐc), 97
moto, S. 1997, ApJ, 478, 631
Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637 (SST)
Hunter, C. 1977, ApJ, 218, 834
Tomisaka, K. 1998, ApJ, 502, L163
Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943
Tscharnuter, W. M., & Winkler, K.-H. 1979, Comput. Phys. Commun., 18,
Kenyon, S. J., Hartmann, L., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869
171
Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676
Winkler, K.-H. A., & Newman, M. J. 1980a, ApJ, 236, 201 (WN)
Kikuchi, N., Nakamoto, T., & Ogochi, K. 2000, PASJ, submitted
ÈÈÈ. 1980b, ApJ, 238, 311 (WN)
Larson, R. B. 1969, MNRAS, 145, 271 (L69)
Yorke, H. W., Bodenheimer, P., & Laughlin, G. 1993, ApJ, 411, 274
Leung, C. M. 1975, ApJ, 199, 340
ÈÈÈ. 1995, ApJ, 443, 199