Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://ip.rsu.ru/~vorobyov/collapse.pdf Äàòà èçìåíåíèÿ: Thu Sep 27 10:44:24 2007 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:45:07 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: http www.astronomy.ru forum index.php topic 4644.0.html

Mon. Not. R. Astron. Soc. 000, 000­000 (0000)

Printed 10 August 2005

A (MN L TEX style file v2.2)

The effect of a finite mass reservoir on the collapse of spherical isothermal clouds and the evolution of protostellar accretion
E. I. Vorobyov
1 2

1, 2

Department of Physics and Astronomy, University of Western Ontario, London, Ontario, N6A 3K7, Canada Institute of Physics, Stachki 194, Rostov-on-Don, Russia

and Shantanu Basu

1

Submitted November 1, 2004

Motivated by recent observations which detect an outer boundary for starless cores, and evidence for time-dependent mass accretion in the Class 0 and Class I protostellar phases, we reexamine the case of spherical isothermal collapse in the case of a finite mass reservoir. The presence of a core boundary, implemented through a constant volume approximation in our simulation, results in the generation of an inward propagating rarefaction wave. This steepens the gas density profile from r-2 (self-similar value) to r-3 or steeper. After a protostar forms, the mass accretion rate M evolves , which is a nonthrough three distinct phases: (1) an early phase of decline in M self-similar effect due to rapid and spatially nonuniform infall in the prestellar phase; (2) for large cores, an intermediate phase of near-constant M from the infall of the outer part of the self-similar density profile which has low (subsonic) infall speed in the prestellar phase; (3) a late phase of rapid decline in M when accretion occurs from the region affected by the inward propagating rarefaction wave. Our model clouds of small to intermediate size make a direct transition from phase (1) to phase (3) above. Both the first and second phase (if the latter is indeed present) are characterized by a temporally increasing bolometric luminosity Lbol , while Lbol is decreasing in the third (final) phase. We identify the period of temporally increasing Lbol with the Class 0 phase, and the later period of terminal accretion and decreasing Lbol with the Class I phase. The peak in Lbol corresponds to the evolutionary time when 50% ± 10% of the cloud mass has been accreted by the protostar. This is in agreement with the classification scheme proposed by Andr´ et al. (1993); our model adds a physical context to e their interpretation. We show how our results can be used to explain tracks of envelope mass Menv versus Lbol for protostars in Taurus and Ophiuchus. We also develop an analytic formalism which successfully reproduces the protostellar accretion rate from profiles of density and infall speed in the prestellar phase. It shows that the spatial gradient of infall speed that develops in the prestellar phase is a primary cause of the temporal decline in M during the early phase of protostellar accretion. Key words: hydrodynamics ­ ISM: clouds ­ stars: formation.

ABSTRACT

1

INTRODUCTION

Recent submillimeter and mid-infrared observations suggest that prestellar cores within a larger molecular cloud are characterized by a non-uniform radial gas density distribution (Ward-Thompson et al. 1999; Bacmann et al. 2000). Specifically, a flat density profile in the central region of size Rflat is enclosed within a region of approximately r-1

E-mail: vorobyov@astro.uwo.ca (EIV); basu@astro.uwo.ca (SB) c 0000 RAS

column density profile (and by implication an r-2 density profile) of extent Rmid . Beyond this, a region of steeper density ( r-3 or greater) is sometimes detected. Finally, at a distance Redge , the column density N seems to merge into a background, and fluctuate about a mean value that is typical for the ambient molecular cloud. The first two regions, of extent Rflat and Rmid , respectively, are consistent with models of unbounded isothermal equilibria or isothermal self-similar gravitational collapse (e.g. Chandrasekhar 1939; Larson 1969; Penston 1969). In either case, the effect


2

E. I. Vorobyov and S. Basu
the assumption of constant mass and volume of a gravitationally contracting core can affect the mass accretion rate and other observable properties after the formation of the central hydrostatic stellar core. A very important question is: which of the two effects mentioned above - a gradient of infall speed in the prestellar phase, or a finite mass reservoir and associated steep outer density slope - is more relevant to explaining the observations of BATC? The evolutionary tracks of envelope mass Menv versus bolometric luminosity Lbol are another important diagnostic of protostellar evolution (Andr´ et al. 2000). BATC have fit the data using a toy e model in which M decreases with time in exact proportion to the remaining envelope mass Menv , i.e. M = Menv / , where is a characteristic time. We seek to explain the observed YSO evolutionary tracks using a physical (albeit highly simplified) model. We perform high resolution one-dimensional spherical isothermal simulations. The initial peak and decline in the mass accretion rate is modeled through numerical simulations and a simplified semi-analytic approach. A second late-time de cline in M due to a gas rarefaction wave propagating inward from the outer edge of a contracting core, is also studied in detail. Comparisons are made with the observationally inferred decrease of mass accretion rate (BATC), and evolutionary tracks of Menv versus bolometric luminosity Lbol (from Motte & Andr´ 2001). e Numerical simulations of spherical collapse of isothermal cloud cores are described in § 2. The comparison of the model with observations is given in § 3. Our main conclusions are summarized in § 4. An analytical approach for the determination of the mass accretion rate is presented in the Appendix.

of an outer boundary is considered to be infinitely far away (i.e. Rmid in our terminology). In numerical simulations of gravitational collapse in which there is a qualitative change in the physics beyond some radius (e.g. a transition from magnetically supercritical to subcritical mass-toflux ratio: Ciolek & Mouschovias 1993; Basu & Mouschovias 1994), the development of a very steep outer density profile is also seen. Finally, larger scale simulations of core formation in clouds with uniform background column density (Basu & Ciolek 2004) show an eventual merger into a near-uniform background column density, demonstrating the existence of Redge . The implication of an outer density profile steeper than r-3 is that there is a finite reservoir of mass to build the star(s), assuming that the gas beyond Redge is governed by the dynamics and gravity of the parent cloud, and thus does not accrete on to the star(s) formed within the core. An important constraint of the observations are the actual sizes of the cores. For example, in the clustered star formation regions such as Ophiuchi protocluster, Redge 5000 AU, and Redge /Rflat 5, while in the more 20000 AU, extended cores in Taurus, 5000 AU Redge 10 (see Andr´ et al. 1999; Andr´ e e, and 5 Redge /Rflat Ward-Thompson, & Barsony 2000). Clearly, only the latter case may approach self-similar conditions. Once a central hydrostatic stellar core has formed, the mass accretion rate is expected to be constant in isothermal similarity solutions (Shu 1977; Hunter 1977; Whitworth & Summers 1985). For example, for the collapse from rest of a singular isothermal sphere (SIS) with density profile SIS = c2 /(2 Gr 2 ), where cs is the isothermal sound speed, s Shu (1977) has shown that the mass accretion rate ( M ) is constant and equal to 0.975 c3 /G. However, two effects can s work against a constant M in more realistic scenarios of isothermal collapse: (1) inward speeds in the prestellar phase are not spatially uniform as in the similarity solutions, and tend to increase inward, meaning that inner mass shells fall in with a greater accretion rate; (2) the effect of a finite mass reservoir will ultimately reduce accretion. The first effect has been clearly documented in a series of papers (e.g. Hunter 1977; Foster & Chevalier 1993; Tomisaka 1996; Basu 1997; Ciolek & K¨ onigl 1998; Ogino, Tomisaka, & Nakamura 1999). It is always present since the outer boundary condition for collapse is distinct from the inner limit of self-similar supersonic infall found in the Larson-Penston solution. Rather, the outer boundary condition must represent the ambient conditions of a molecular cloud, which do not correspond to large-scale infall (Zuckerman & Evans 1974). Additionally, the finite mass reservoir and steeper than r-3 profile as a source of the declining accretion rate has been studied analytically by Henriksen, Andr´ & Bontemps (1997) and e, Whitworth & Ward-Thompson (2001), although they did not account for the physical origin of such a steep density slope. Indeed, a study of outflow activity from young stellar ob jects (YSO's) by Bontemps et al. (1996; hereafter BATC) suggests that M declines significantly with time during the accretion phase of protostellar evolution. Specifically, BATC have shown that if the CO outflow rate is proportional to M , then Class 0 ob jects (young protostars at the begin ning of the main accretion phase) have an M that is factor of 10 greater (on average) than that of the more evolved Class I ob jects. In this paper, we investigate in detail how

2 2.1

ISOTHERMAL COLLAPSE Mo del Assumptions

We consider the gravitational collapse of spherical isothermal (temperature T = 10 K) clouds composed of molecular hydrogen with a 10% admixture of atomic helium. The models actually represent cloud cores which are embedded within a larger molecular cloud. The evolution is calculated by solving the hydrodynamic equations in spherical coordinates: 1 +2 r2 vr t r r 1 r2 vr vr (vr ) + 2 t r r e 1 2 +2 (r evr ) t r r = = = 0 p GM - 2 , r r p 2 -2 ( r vr ) r r - (1) (2) (3)

where is the density, vr is the radial velocity, M is the enclosed mass, e is the internal energy density and p = ( - 1)e is the gas pressure. The ratio of specific heats is equal to = 1.001 for the gas number density n 1011 cm-3 , which implies isothermality (the value of is not exactly unity in our implementation in order to avoid a division by zero). We define the gas number density n = /m, where m = 2.33 mH is the mean molecular mass. When the gas number density in the collapsing core exceeds 1011 cm-3 , we form the central hydrostatic stellar core by imposing an adiabatic index = 5/3. This simplified treatment of the transition to an opaque
c 0000 RAS, MNRAS 000, 000­000


Finite Mass Reservoir Col lapse
protostar misses the details of the physics on small scales. Specifically, a proper treatment of the accretion shock and radiative transfer effects is required to accurately predict the properties of the stellar core (see Winkler & Newman 1980 for a detailed treatment and review of work in this area). However, our method should be adequate to study the protostellar accretion rate, and has been used successfully by e.g. Foster & Chevalier (1993) and Ogino et al. (1999) for this purpose. We use the method of finite-differences, with the time-explicit, operator split solution procedure used in the ZEUS-1D numerical hydrodynamics code; it is described in detail by Stone & Norman (1992). We have introduced the momentum density correction factor, as advocated by M¨ hmeyer & Muller (1989), to avoid the development of onc ¨ an anomalous density spike at the origin (see Vorobyov & Tarafdar 1999 for details). The numerical grid has 700 points which are initially uniformly spaced, but then move with the gas until the central stellar core is formed. This provides an adequate resolution throughout the simulations. We impose boundary conditions such that the gravitationally bound cloud core has a constant mass and volume. The assumption of a constant mass appears to be observationally justified by the sharp outer density profiles described in § 1. Physically, this assumption may be justified if the core decouples from the rest of a comparatively static, diffuse cloud due to a shorter dynamical timescale in the gravitationally contracting central condensation than in the external region. A specific example of this, due to enhanced magnetic support in the outer envelope, is found in the models of ambipolar-diffusion induced core formation (see, e.g. Basu & Mouschovias 1995). The assumption of a constant volume is mainly an assumption of a constant radius of gravitational influence of a cloud core within a larger parent diffuse cloud. The radial gas density distribution of a self-gravitating cloud with finite central density that is in hydrostatic equilibrium (e.g. Chandrasekhar 1939) can be conveniently approximated by a modified isothermal sphere, with gas density = c 1 + (r/rc )
2

3

Table 1. Model parameters Model I1 I2 n
c

r

c

r

out

c out

M

cl

T 10 10

5 .0 5 .0

0.033 0.033

0.16 0.5

24 324

5 24
-3

All number densities are in units of 104 cm masses in M , and temperatures in K.

, lengths are in pc,

sity perturbation in order to initiate evolution. Use of the modified isothermal sphere simplifies the analysis a little bit since there is a clear transition from flat central region to a power-law outer profile. The choice of central density c and outer radius rout determines the cloud mass. We study many different cloud masses - two models are presented in this section and other models are used to fit observational tracks in § 3. We also add a (small) positive density perturbation of a factor of 1.1 (i.e. the initial gas density distribution is increased by a factor of 1.1) to drive the cloud (especially the inner region which is otherwise near-equilibrium) into gravitational collapse. Table 1 shows the parameters for two model clouds presented in this section. The adopted central number density nc = 5 â 104 cm-3 is roughly an order of magnitude lower than is observed in prestellar cores (Ward-Thompson et al. 1999). Considering that these cores may be already in the process of slow gravitational contraction, our choice of nc is justified for the purpose of describing the basic features of star formation. In both models, the outer radius rout is chosen so as to form gravitationally unstable prestellar cores with central-to-surface density ratio c /out 14 (since our initial states are similar to Bonnor-Ebert spheres). In model I1, c /out 24 and by implication rout /rc 4.7, whereas in model I2 c /out 324 and rout /rc 18. Model I2 thus represents a very extended prestellar core. Models I1 and I2 have masses 5 M and 24 M respectively; the `I' stands for isothermal. 2.2 Numerical Results

(4)

(Binney & Tremaine 1987), where c is the central density is and rc a radial scale length. We choose a value rc = 1.1 cs / Gc , so that the inner profile is close to that of a Bonnor-Ebert sphere, rc is comparable to the Jeans length, and the asymptotic density profile is 2.2 times the equilibrium singular isothermal sphere value SIS = c2 /(2 Gr 2 ). s The latter is justified on the grounds that core formation should occur in a somewhat non-equilibrium manner (an extreme case is the Larson-Penston flow, in which case the asymptotic density profile is as high as 4.4 SIS ), and also by observations of protostellar envelope density profiles that are e, often overdense compared to SIS (Andr´ Motte, & Belloche 2001). For small radii (r rc ), the initial density is very close to the equilibrium solution for an isothermal sphere with a finite central density. However, at large radii it is twice the value of the equilibrium isothermal sphere, which converges to SIS . Hence, our initial conditions resemble those of other workers (Foster & Chevalier 1993; Ogino et al. 1999) who start with Bonnor-Ebert spheres and add a positive denc 0000 RAS, MNRAS 000, 000­000

Fig. 1 shows the temporal evolution of the radial gas density profiles (the upper panel) and velocity profiles (the lower panel) during the runaway collapse phase (before the formation of the central hydrostatic stellar core) in model I1. The density and velocity profiles are numbered according to evolutionary sequence, starting from the initial distributions (profile 1; note that the cloud core is initially at rest) and ending with those obtained when the central number density has almost reached 1010 cm-3 (profile 5). The dashed lines in the upper panel of Fig. 1 show the power-law index d ln /d ln r of the gas distribution for profiles 1, 2, and 3. By the time that a relatively mild center-to-boundary density contrast 150 is established (profile 2), the radial density profile starts resembling those observed in Taurus by Bacmann et al. (2000): it is flat in the central region, then gradually changes to an r-2 profile, and falls off as r-3 or steeper in the envelope at r 0.08 pc. The sharp change in slope of the density profile (e.g. at r 0.08 pc in profile 2 of Fig. 1) is due to an inward ly-propagating gas rarefaction wave caused by a finite reservoir of mass. The self-similar re-


4

E. I. Vorobyov and S. Basu

Figure 1. The radial gas density (the upper panel) and velocity (the lower panel) profiles obtained in model I1 before the central hydrostatic stellar core is formed. The number 1 corresponds to the initial profiles (note that initially the cloud core is at rest) and the number 5 labels profiles when the central gas number density has almost reached 1010 cm-3 . The dashed lines in the upper panel show the power-law index d ln /d ln r of the gas distributions 1, 2, and 3. Profiles 2, 3, 4, and 5 are reached at times 0.309 Myr, 0.378 Myr, 0.392 Myr, and 0.394 Myr, respectively. For reference, the dotted line is the density profile of a singular isothermal sphere.

gion with r-2 density profile is of the Larson-Penston type, with density somewhat greater than the equilibrium singular isothermal sphere value (c2 /2 Gr 2 ). The velocity profiles s in Fig. 1 also show a distinct break at the instantaneous location of the rarefaction wave. Furthermore, the peak infall speed is clearly supersonic (since cs = 0.19 km s-1 ) by the time profile 4 is established, again consistent with LarsonPenston type flow in the inner region. Fig. 2a shows the temporal evolution of the accretion rate at a radial distance of 600 AU from the center in model I1.1 The evolution is characterized by a slow initial gravitational contraction and then a very rapid increase until about 0.4 Myr. Subsequently, a central hydrostatic stellar core forms and the mass accretion rate reaches its maximum value of 2.8 â 10-5 M yr-1 (or 17.4 c3 /G). After stellar s core formation, the evolution of the mass accretion rate has possibly three distinct phases, of which two are on display
1

Figure 2. The temporal evolution of the mass accretion rate at the radial distance r = 600 AU from the center of a cloud core obtained in a) model I1 and b) model I2. The model cloud I1 has mass 5 M and the model cloud I2 has mass 24 M . The solid lines show M during the runaway collapse phase, prior to the formation of a central stellar core. The dashed and dotted lines plot M after stellar core formation; the dashed lines show M before the gas affected by the inwardly propagating rarefaction wave has reached r = 600 AU, whereas the dotted lines show M after this gas has reached r = 600 AU. The numbers in parentheses reflect the percentage of cloud mass remaining in the envelope at the given times.

We note that the accretion rate is not expected to vary significantly in the range 0.1 AU < r < 1000 AU according to Masunaga & Inutsuka (2000).

in Fig. 2a. The early phase, plotted with the dashed line in Fig. 2a, is characterized by accretion of material that has not yet been affected by the rarefaction wave propagating inward from the outer boundary. The accretion rate is declining, even though the density profile near the stellar core was nearly self-similar at the moment of its formation. This decline is due to the gradient of infall velocity in the inner regions, an effect not predicted in the similarity solutions. However, if there is a large outer region with mass shells that are falling in at significantly subsonic speeds when the central stellar core forms (see discussion of Fig. 2b below), the accretion rate will eventually stabilize to a constant value that is consistent with the standard theory of Shu (1977). In that picture, progressively higher shells of gas lose their partial pressure support and start falling from rest on to the central stellar core almost in a free-fall manner. This would be the intermediate phase of accretion. However, the late phase of very rapid decline of the accretion rate starts at roughly 0.46 Myr (before the intermediate
c 0000 RAS, MNRAS 000, 000­000


Finite Mass Reservoir Col lapse
phase can be established in the 5 M cloud), when gas affected by the inwardly propagating rarefaction wave reaches the inner 600 AU. This results in a sharp drop of M as shown in Fig. 2a by the dotted line. The existence of the (in principle) three distinct phases of mass accretion is clearly seen Fig. 2b, where M of the more extended cloud (rout /rc 18) is plotted (hereafter, model I2). The outer boundary is now at rout = 0.5 pc and it takes a time 1 Myr for the influence of the rarefaction wave to reach the inner 600 AU. As a result, the mass accretion rate has time to stabilize at a constant value of 1.34 â 10-5 M yr-1 (the dashed line in Fig. 2b), before it sharply drops at later times (the dotted line in Fig. 2b). According to Shu (1977), the collapse from rest of a power-law profile that has a density equal to twice SIS yields a mass accretion rate 5.58 c3 /G = 8.86 â 10-6 M yr-1 . Our stable s intermediate accretion rate is roughly consistent with this prediction since the density in the power-law tail is actually somewhat greater than twice SIS . It is equal to 2.42 SIS in the initial state, and grows to greater overdensities in the innermost regions. However, the bulk of the matter, which is in the outer tail, has density within 2.5SIS . Further experiments with our numerical simulations show that the intermediate phase of constant accretion rate is observed only in rather extended prestellar cores with rout /rc 15. Foster & Chevalier (1993) found an even stronger criterion 20. Since more extended cores tend to be more rout /rc massive as well, we may expect to observe the intermediate phase more frequently in the collapse of massive cores.

5

Figure 3. The temporal evolution of the mass accretion rate. The solid line shows the results of isothermal numerical simulations (model I1). The dashed line shows M obtained in a pressure-free approximation if such collapse begins from rest with the relatively mildly concentrated density profile 2 in Fig. 1. The dotted line shows the result for pressure-free collapse if a non-zero initial velocity (that of profile 2 in the bottom panel of Fig. 1) is also used. The agreement of the pressure-free model and full numerical simulation are quite good in the latter case.

2.4

Semi-analytic Mo del

2.3

Effect of Boundary Condition

Our standard simulation does not contain an external medium explicitly. In order to explore the effect of such a medium, we ran additional simulations in which the cloud core is surrounded by a spherical shell of diffuse (i.e. nongravitating) gas of constant temperature and density. The outermost layer of the cloud core and the external gas are initially in pressure balance. We found that the value of M in the late accretion phase may depend on the assumed values of the external density and temperature. For instance, if the gravitating core is nested within a larger diffuse nongravitating cloud of T = 10 K and = out , the accretion rate increases slightly as compared to that shown in Fig. 2 by the dotted line. A warmer external non-gravitating environment of T = 200 K and = out /20 shortens the duration of the late accretion phase shown in Fig. 2 by the dotted line. This phase may be virtually absent if the sound speed of the external diffuse medium is considerably higher (by a factor 1000) than that of the gravitationally bound core. This essentially corresponds to a constant outer pressure boundary condition (see Foster & Chevalier 1993). However, such a high sound speed contrast is not expected for star formation taking place in a dense (n 104 cm-3 ) environment like Ophiuchi (Johnstone et al. 2000). We believe that the constant volume boundary condition, and resulting inward propagating rarefaction wave, are best at reproducing the steep outer density profiles and the low (residual) mass accretion rate necessary to explain the Class I phase of protostellar accretion.
c 0000 RAS, MNRAS 000, 000­000

Finally, we compute M of a pressure-free cloud using the analytical approach developed in the Appendix. This ap proach allows for the determination of M for a cloud with given initial radial density (r0 ) and velocity v0 (r0 ) profiles, if the subsequent collapse is pressure-free. We find that the success or failure of the analytical approach to describe the mass accretion rate of the isothermal cloud depends on the adopted (r0 ) and v0 (r0 ) profiles. For instance, if (r0 ) is determined by profile 2 (the upper panel of Fig. 1) and v0 (r0 ) = 0, the pressure-free mass accretion rate shown in Fig. 3 by the dashed line reproduces only very roughly the main features of the isothermal accretion rate (the solid line in Fig. 3). However, if we take into account the non-zero and non-uniform velocity profile v0 (r0 ) plotted in the lower panel of Fig. 1 (profile 2), then the pressure-free M shown by the dotted line in Fig. 3 reproduces that of the isothermal cloud much better. This example demonstrates the importance of the velocity field prior to stellar core formation in determining the accretion rates after its formation. The success of our analytical pressure-free approach also shows that the collapse of the isothermal cloud can be regarded as essentially pressure-free from the time of a relatively mild central concentration c /out 150, when the central number density 2 â 105 cm-3 .

3

ASTROPHYSICAL IMPLICATIONS

Class 0 ob jects represent a very early phase of protostellar evolution (see Andr´ et al. 2000), as evidenced by a relae tively high ratio of submillimeter luminosity to bolometric luminosity: Lsubmm /Lbol > 0.5%. Class 0 ob jects also drive powerful collimated CO outflows. A study of outflow ac-


6

E. I. Vorobyov and S. Basu
Table 2. Model parameters for Ophiuchus M
cl

nc 1 5 2 2 â â â â 107 106 106 106

r

out

c / 15.0 14.1 18.0 18.0

out

r

out

/rc

2.0 2.0 2.4 4.0

0.17 0.23 0.55 0.9

1600 1900 4000 4000

3.7 3.6 4.1 4.1

All number densities are in cm M.

-3

, lengths in AU, and masses in

Table 3. Model parameters for Taurus M Figure 4. CO momentum flux Fco versus envelope mass Menv for 41 sources of Bontemps et al. (1996). The Class 0 and Class I ob jects are plotted with the open and filled circles, respectively. The model Fco - Menv tracks of three prestellar clouds of Mcl = 0.2 M , 0.7 M , and 2 M are shown by the dotted, dashed, and solid lines, respectively.
cl

nc 1 1 1 8 .5 â 10 â 106 â 106 â 105
6

r

out

c /out 18 26 71 73
-3

r

out

/r

c

1.9 1.8 1.5 1.8

0.46 0.65 1.0 1.5

5000 6000 10000 12000

4.1 5.0 8.4 8.5

All number densities are in cm M.

, lengths in AU, and masses in

tivity in low-mass YSO's by BATC suggests that the CO momentum flux Fco declines significantly during protostellar evolution. Specifically, Fco decreases on average by more than an order of magnitude from Class 0 to Class I ob jects. This tendency is illustrated in Fig. 4, where we plot Fco ver sus Menv for 41 sources listed in BATC. We relate Fco to M by F
co

=f

ent

â ( M w / M ) Vw â M ,

(5)

where fent is the entrainment efficiency that relates Fco to the momentum flux Mw Vw of the wind. Based on theoretical models in the literature, BATC suggested that the factor fent and the outflow driving engine efficiency (Mw /M )Vw do not vary significantly during protostellar evolution. This implies that the observed decline of Fco reflects a corre sponding decrease in M from the Class 0 to the Class I stage. Following BATC, we take fent = 1, Mw /M = 0.1, -1 and Vw = 150 km s and use equation (5) to compute Fco from our model's known mass accretion rate M (see Fig. 2). Since the sample of Class 0 and Class I ob jects listed in BATC includes sources from both the Ophiuchus and Taurus star forming regions, we develop model clouds which take into account the seemingly different initial conditions of star formation in these regions. As mentioned in § 1, the two most prominent differences between these two regions are: (1) The cores in Ophiuchus have outer radii 5000 AU) which are smaller than in Taurus, where (rout 20000 AU (Andr´ et al. 1999; Andr´ e e 5000 AU rout et al. 2000); (2) The radial column density profiles of the protostellar envelopes of Class 0 ob jects in Ophiuchus are at least 2-3 times denser than a SIS at T = 10 K, whereas in Taurus the protostellar envelopes are overdense compared to the SIS by a smaller factor 2 (Andr´ et al. 2001). e This implies that radial column density profiles of prestel lar cores in Ophiuchus and Taurus may follow the same tendency. We develop a set of Ophiuchus model cores which have rout 5000 AU, and a set of Taurus model cores which have 5000 AU rout 20000 AU. Furthermore, the factor

(by which our model density profiles are asymptotically overdense compared to SIS ) is taken to be 2.0 for Ophiuchus and < 2.0 for Taurus. Clearly, there is no unique set of model cloud parameters that would be exclusively consistent with the observational data, given the measurement uncertainties. We have chosen a set of core central densities c , radii rout , and overdensity factors so as to reasonably reproduce the observed properties of the cores in the two regions. The parameters of the model density distributions for Ophiuchus and Taurus are listed in Table 2 and Table 3, respectively. We also note that we have ensured that the cores satisfy the gravitational instability criterion rout /rc 3.6, which is similar to that for Bonnor-Ebert spheres. The Ophiuchus model cores are clustered near this limiting value of rout /rc , but the Taurus model cores are allowed to be somewhat more extended, again in keeping with observed properties. We also note that the masses of prestellar cores with the radial density profile given by equation (4) scale as 1/0.5 , c if the ratio rout /rc is fixed. The sample of 41 sources in BATC contains Class 0 and Class I ob jects from both Ophiuchus and Taurus. Hence, in Fig. 4 we take three representative prestellar clouds of Mcl = 0.23 M (Ophiuchus), 0.65 M (Taurus), and 2.0 M (Taurus), for which the Fco - Menv tracks are shown by the dotted, dashed, and solid lines, respectively. Both the data and model tracks show a near-linear correlation be tween Fco ( M ) and Menv . A slightly better fit of the model tracks to the data can be obtained by adjusting one or more of the estimated parameters fent , Mw /M , and Vw by factors of order unity. Based on the near-linear correlation of Fco and Menv , BATC developed a toy model in which M decreases with time in exact proportion to the remaining envelope mass Menv , i.e. M = Menv / , where is a characteristic time. Furthermore, if one assumes that the bolometric luminosity derives entirely from the accretion on to the hydrostatic stellar core, i.e., Lbol = GMc M /Rc , where Mc and Rc are
c 0000 RAS, MNRAS 000, 000­000


Finite Mass Reservoir Col lapse

7

Figure 5. Temporal evolution of the bolometric luminosity Lbol and CO momentum flux Fco . a) The solid and dashed lines show Lbol and Fco obtained in model I1, respectively. b) The solid and dashed lines show Lbol and Fco obtained in model I2, respectively. Note that Lbol is still increasing during the early phase of accretion rate decline but only declines later due to the more severe accretion rate decline caused by the inward propagating rarefaction wave. The vertical dash-dotted line is the temporal dividing line between the Class 0 and Class I phases for each model; the numbers below give the mass of the central stellar core as a percentage of the total cloud mass. Crosses indicate the time when 50% of the initial cloud mass has been accreted by the protostar. Note the use of a logarithmic scale for time, so that the Class I phase is still longer than the Class 0 phase for model I2.

the mass and radius of the stellar core, respectively, then the bolometric luminosity reaches a maximum value when half of the initial prestellar mass has been accreted by the protostar and the other half remains in the envelope. The evolutionary time when Mc = Menv was defined by Andr´ e et al. (1993) as the conceptual border between the Class 0 and Class I evolutionary stages. The solid and dashed lines in Fig. 5a and Fig. 5b show Lbol and Fco obtained in model I1 and model I2, respectively. Since we do not follow the evolution of a protostar to the formation of the second (atomic) hydrostatic core, we take Rc = 3 R and let Lbol = GMc M /Rc . The radius Rc depends on the accretion rate and stellar mass (see Fig. 7 of Stahler 1988) and may vary from Rc 1.5 R for small stellar cores Mc 0.2 M and low accretion rates M 2 â 10-6 M yr-1 to Rc 5.0 R for large stellar cores Mc 1.0 M and high accretion rates M 1 â 10-5 M yr-1 . However, this variation constitutes roughly a factor of 2 change in the adopted average value of Rc = 3 R . Indeed, we performed numerical simc 0000 RAS, MNRAS 000, 000­000

ulations with a varying Rc (assuming a normal deuterium abundance) and found that it has only a minor qualitative influence on our main results. The stellar core mass Mc is computed by summing up the masses of the central hydrostatic spherical layers in our numerical simulations. An obvious difference in the temporal evolution of Fco and Lbol is seen in Fig. 5. The temporal evolution of Fco after the central hydrostatic core formation at t 0.4 Myr goes through the same phases as shown for M in Fig. 2. The temporal evolution of Lbol shows two distinct phases: it increases during the early phase (unlike Fco ) and starts decreasing only when gas affected by the inward propagating rarefaction wave reaches the central hydrostatic core. Thus, in our model, only the rarefaction wave acts to reduce Lbol during the accretion phase of protostellar evolution. This is a physical explanation for the peak in Lbol that also occurs in the toy model of BATC. In that model, the bolometric luminosity reaches a maximum value when exactly half of the initial prestellar mass has been accreted by the protostar. In our simulations, the peak in Lbol corresponds to the evolutionary time when 50% ± 10% of the matter is in the protostar (higher deviations up to +15% are found in very massive and extended prestellar clouds). Finally, in Fig. 6 we show the Menv - Lbol evolutionary tracks. We use eight representative prestellar cloud core masses as listed in Table 2 and Table 3. Fig. 6a shows the overlaid data for YSO's in Taurus, while Fig. 6b has overlaid data for Ophiuchus. The data for both samples are taken from Motte & Andr´ (2001). The open circles repe resent bonafide Class 0 ob jects, the solid circles represent the bonafide Class I ob jects, while the triangles represent the so-called peculiar Class I ob jects observed in Taurus. We note that the envelope masses given in Table 2 of Motte & Andr´ (2001) and plotted in their Fig. 5 and Fig. 6 are dee termined within a 4200 AU radius circle. While this should relatively well describe the total envelope masses in Ophiuchus, a substantial (a factor of 3) portion of the envelope mass may be missing in the Taurus cores, which have sizes as large as 15000-20000 AU. For this reason, we plot in Fig. 6a the total envelope masses given in Table 4 of Motte & Andr´ e (2001) for a set of resolved Taurus cores. The loci of maximum Lbol in the Menv - Lbol tracks roughly separate two phases in the evolution of a protostar: a shorter one characterized by accretion of matter from the envelope not yet affected by the rarefaction wave (i.e. characterized by the r-2 gas density profile or shallower) and a longer one characterized by accretion of matter from the rarefied envelope (i.e. characterized by the r-3 profile or steeper). The turnover also corresponds to the evolutionary time when 50% ± 10% of the matter is in the protostar and a corresponding amount remains in the envelope as shown by the crosses in Fig. 6. This is in agreement with the observational requirements and toy model of BATC. Given that the peak in Lbol is our conceptual dividing line between two distinct phases of accretion, we conclude that in Taurus, most of the so-called Class I ob jects would tend to fall into the Class 0 category in our scheme. They may indeed be more evolved than the already identified Class 0 ob jects, having lower values of M and Menv , but would not be in a qualitatively distinct phase of evolution (see Motte & Andr´ 2001 e for a similar conclusion). In contrast, the so-called peculiar Class I ob jects in Taurus would be proper Class I ob jects in


8

E. I. Vorobyov and S. Basu
One problem should be pointed out here. While our model Lbol - Menv tracks in Fig. 6 explain well the measured bolometric luminosities in Ophiuchus, they seem to overestimate Lbol in Taurus by a factor of 5-10. This is the so-called "luminosity problem" that was first noticed by Kenyon, Calvet, & Hartmann (1993). As a consequence, the position near the turnover in Lbol - Menv tracks for Taurus is scarcely populated. This implies that while spherical collapse models may be appropriate for the determination of Lbol in Ophiuchus, they tend to overestimate Lbol in Taurus. It is possible that a significant magnetic regulation of the early stages of star formation in Taurus, as implied by e.g. polarization maps (Moneti et al. 1984) would yield more flattened envelopes which result in a lower accretion rate on to the central protostar and a smaller bolometric luminosity. Interestingly, Kenyon et al. (1993) also concluded that envelopes in Taurus should be highly flattened in order to explain their spectral energy distribution. Two dimensional simulations are required to address this issue. Finally, it is worth noting that our Taurus model cores are generally more massive than the Ophiuchus model cores. This in agreement with observations, and can be justified theoretically on the basis of a lower mean column density (hence greater Jeans length and Jeans mass of a sheetlike configuration) in Taurus compared to regions of more clustered star formation in e.g. Ophiuchus and Orion. Taken at face value, our models then imply that Taurus protostars should be more massive in general than Ophiuchus protostars. While such a conclusion must be tempered by the fact that we do not model magnetic support or feedback from outflows, there is some evidence that Taurus does have a significantly higher mass peak in its initial mass function than does the Trapezium cluster in Orion (see Luhman 2004 and references within).

Figure 6. Envelope mass Menv versus bolometric luminosity Lbol for 37 protostellar ob jects taken from Motte & Andr´ (2001). Die agrams are shown for a) Taurus, and b) Ophiuchus. The Class 0 and Class I ob jects are plotted with the open and filled circles, respectively. The triangles represent the observed peculiar Class I sources. The model Menv - Lbol tracks of eight prestellar clouds with masses Mcl = 1.5 M , 1.0 M , 0.65 M , and 0.46 M (Taurus) and 0.9 M , 0.55 M , 0.23 M , and 0.17 M (Ophiuchus) are shown by the solid, dashed, dotted, and dotted-dashed lines, respectively. Crosses indicate the time when 50% of the cloud mass has been accreted by the protostar.

4

CONCLUSIONS

our scheme since they are likely in a phase of declining Lbol . In Ophiuchus, the currently identified Class 0 and Class I ob jects do seem to fall on two distinct sides of the peak in Lbol . Fig. 5 indicates that in extended clouds (as of model I2) the phase of increasing Lbol is longer than in compact clouds (as of model I1). This may explain why this phase is more populated in Taurus than in Ophiuchus. In addition, extended clouds have a longer phase of accretion from the envelope not yet affected by the rarefaction wave and, as a consequence, a higher probability of having a quasi-constant accretion phase. This is in agreement with the previous suggestions made by Henriksen et al. (1997) and Andr´ et al. (2000) e that the accretion history in Taurus is closer to the SIS scenario than in Ophiuchus. However, we note that none of our representative prestellar cores listed in Table 3 and used to 15) fit the data for Taurus are extended enough (rout /rc to have a distinct phase of constant accretion.

Our numerical simulations indicate that the assumption of a finite mass reservoir of prestellar cores is required to explain the observed Class 0 to Class I transition. We start our collapse calculations by perturbing a modified isothermal sphere profile (eq. [4]) that is truncated and resembles a bounded isothermal equilibrium state. Specifically, we find that · Starting in the prestellar runaway collapse phase, a shortage of matter developing at the outer edge of a core generates an inward propagating rarefaction wave that steepens the radial gas density profile in the envelope from r-2 to r-3 or even steeper. · After a central hydrostatic stellar core has formed, and the cloud core has entered the accretion phase, the mass ac cretion rate M on to the central protostar can be divided into three possible distinct phases. In the early phase, M decreases due to a gradient of infall speed that developed during the runaway collapse phase (such a gradient is not predicted in isothermal similarity solutions). An intermediate phase of near-constant M follows if the core is large enough to have an extended zone of self-similar density profile with relatively low infall speed during the prestellar phase. Finally, when accretion occurs from the region affected by the
c 0000 RAS, MNRAS 000, 000­000


Finite Mass Reservoir Col lapse
inward propagating rarefaction wave, a terminal and rapid decline of M occurs. · A pressure-free analytic formalism for the mass accretion rate can be used to predict the mass accretion rate after stellar core formation, given the density and velocity profiles in a suitably late part of the runaway collapse phase. Our formulas can estimate M at essential ly any radial distance from the central singularity. This makes it possible to obtain M as a function of radial distance at any given time. We have demonstrated the importance of the velocity field of a collapsing cloud in determining M ; our approach successfully estimates the accretion rate if the velocity field is taken into account. It demonstrates that the initial decline in M is due to the gradient of infall speed in the prestellar phase. · From an observational point of view, we can understand evolutionary Menv - Lbol tracks using core models of relatively small mass and size, so that there is not an extensive self-similar region, in agreement with the profiles observed by e.g. Bacmann et al. (2000). This means that in the accretion phase, M makes a direct transition from the early decline phase to the late decline phase when matter is accreted from the region of steep (r-3 or steeper) density profile that is affected by the inward propagating rarefaction wave. In the first phase (which we identify as the true Class 0 phase), the bolometric luminosity Lbol is increasing with time, even though M and the CO momentum flux Fco are slowly decreasing. In the second phase (which we identify as the Class I phase), both Lbol and Fco decline with time. Hence, our simulations imply that the influence of the rarefaction wave roughly traces the conceptual border between the Class 0 and Class I evolutionary stages. Regions of star formation with more extended cores, like Taurus, should reveal a larger fraction of protostars in the phase of increasing Lbol . Our Fig. 6 reveals that this is indeed the case, if most of the so-called Class I ob jects in Taurus are reclassified as Class 0, according to our definition. The so-called peculiar Class I ob jects in Taurus would be bona-fide Class I ob jects according to our definition (see Motte & Andr´ 2001 for a e similar conclusion on empirical grounds). · Luminosities derived entirely from the accretion on to the hydrostatic stellar core tend to be larger than the measured bolometric luminosities Lbol in Taurus by a factor of 5-10, while they seem to better explain the measured Lbol in Ophiuchus. This implies that physical conditions in Ophiuchus may favour a more spherically symmetric star formation scenario. Our results should be interpreted in the context of models of one-dimensional radial infall. They illuminate phenomena which are not included in standard self-similar models of isothermal spherical collapse, by clarifying the importance of boundary (edge) effects in explaining the observed Fco - Menv and Menv - Lbol tracks. Important theoretical questions remain to be answered, such as the nature of the global dynamics of a cloud which could maintain a finite mass reservoir for a core. A transition to a magnetically subcritical envelope may provide the physical boundary that we approximate in our model. For example, Shu, Li, & Allen (2004) have recently calculated the (declining) accretion rate from a subcritical envelope on to a protostar, under the assumption of flux freezing. An alternate or complementary
c 0000 RAS, MNRAS 000, 000­000

9

mechanism of limiting the available mass reservoir is the effect of protostellar outflows. Our main observational inference is that a finite mass reservoir and the resulting phase of residual accretion is necessary to understand the Class I phase of protostellar evolution. Our calculated mass accretion rates really represent the infall onto an inner circumstellar disk that would be formed due to rotation. Hence, our results are relatable to observations if matter is cycled through a circumstellar disk and on to a protostar rapidly enough so that the protostellar accretion is at least proportional to the mass infall rate on to the disk. This is likely, since disk masses are not observed to be greater than protostellar masses, but needs to be addressed with a more complete model. In future papers, we will investigate the role of non-isothermality (using detailed cooling rates due to gas and dust), rotation, magnetic fields, and non-axisymmetry in determining M and implied observable quantities.

ACKNOWLEDGMENTS We thank Sylvain Bontemps, the referee, for an insightful report which led us to make significant improvements to the paper. We also thank Philippe Andr´ for valuable comments e about the observational interpretation of our results. This work was conducted while EIV was supported by the NATO Science Fellowship Program administered by the Natural Sciences and Engineering Research Council (NSERC) of Canada. EIV also gratefully acknowledges present support from a CITA National Fellowship. SB was supported by a research grant from NSERC.

REFERENCES Andr´ P., Motte, F., Bacmann, A., Belloche, A., 1999, in e, Nakamoto, T., ed., Star Formation 1999. Nobeyama Radio Observatory, Nobeyama, p. 145 Andr´ P., Ward-Thompson, D., Barsony, M., 1993, 406, e, 122 Andr´ P., Ward-Thompson, D., Barsony, M., 2000, in e, Mannings, V., Boss, A. P., Russell, eds., Protostars and Planets IV. Univ. Arizona Press, Tucson, p. 59 Andr´ P., Motte, F., Belloche, A., 2001, in Montmerle, T., e, Andr´ P., eds, ASP Conf. Ser. Vol. 243, From Darkness to e, Light. Astron.Soc.Pac., San Francisco, p. 209 Bacmann, A., Andr´ P., Puget, J. L. et al., 2000, A&A, e, 314, 625 Basu, S., 1997, ApJ, 485, 240 Basu, S., Ciolek, G. E., 2004, ApJ, 607, L39 Basu, S., Mouschovias, T. Ch., 1994, ApJ, 432, 720 Basu, S., Mouschovias, T. Ch., 1995, ApJ, 453, 271 Binney, J., Tremaine, S., 1987, Galactic Dynamics. Princeton Univ. Press, Princeton Bonnor, W. B., 1956, MNRAS, 116, 351 Bontemps, S., Andr´ P., Terebey, S., Cabrit, S. 1996, A&A, e, 311, 858 (BATC) Chandrasekhar, S., 1939, An Introduction to the Study of Stellar Structure. Univ. Chicago Press, Chicago Ciolek, G. E., Mouschovias, T. Ch., 1993, ApJ, 418, 774 Ciolek, G. E., K¨ onigl, A., 1998, ApJ, 504, 257


10

E. I. Vorobyov and S. Basu
t= arccos r/r0 + 0.5 sin(2 arccos 2GM (r0 )/r
3 0

Ebert, R., 1957, Afz, 42, 263 Foster, P. N., Chevalier, R. A., 1993, ApJ, 416, 303 Henriksen, R., Andr´ P., Bontemps, S., 1997, A&A, 323, e, 549 Hunter, C., 1962, ApJ, 136, 594 Hunter, C., 1977, ApJ, 218, 834 Johnstone, D., Wilson, C. D., Moriarty-Schieven, G. et al., 2000, ApJ, 545, 327 Kenyon, S. J., Calvet, N., Hartmann, L. 1993, ApJ, 414, 676 Larson, R. B., 1969, MNRAS, 145, 271 Luhman, K. L., 2004, ApJ, 617, 1216 Masunaga, H., Inutsuka, S., 2000, ApJ, 531, 350 Moneti, A., Pipher, J. L., Helfer, H. L., McMillan, R. S., Perry, M. L., 1984, ApJ, 282, 508 M¨ hmeyer, R., Muller, E., 1989, A&A, 217, 351 onc ¨ Motte, F., Andr´ P., 2001, A&A, 365, 440 e, Ogino, S., Tomisaka, K., Nakamura, F., 1999, PASJ, 51, 637 Penston, M. V., 1969, MNRAS, 144, 425 Shu, F. H., 1977, ApJ, 214, 488 Shu, F. H., Adams, F. C., Lizano, S., 1987, ARA&A, 25, 23 Shu, F. H., Li, Z.-Y., Allen, A., 2004, ApJ, 601, 930 Stahler, S. W., 1988, ApJ, 332, 804 Stone, J. M., Norman, M. L., 1992, ApJS, 80, 753 Tomisaka, K., 1996, PASJ, 48, L97 Vorobyov, E. I., Tarafdar, S. P., 1999, A&ATr, 17, 407 Ward-Thompson, D., Motte, F., Andr´ P., 1999, MNRAS, e, 305, 143 Whitworth, A. P., Summers, D., 1985, MNRAS, 214, 1 Whitworth, A. P., Ward-Thompson, D., 2001, ApJ, 547, 317 Winkler, K.-H. A., Newman, M. J., 1980, ApJ, 236, 201 Zuckerman, B., Evans, N. J., 1974, ApJ, 192, L149

r/r0 )

.

(A3)

The velocity v (r, t) at a given radial distance r and time t can now be obtained from equations (A2) and (A3). The values of r and t are sufficient to determine r0 (a value > r but rout , where rout is the radius of a cloud) from equation (A3). Subsequently, we use the obtained value of r0 in equation (A2) to obtain v (r, t). Provided that the shells do not pass through each other (i.e. the mass of a moving shell is conserved, dM (r, t) = dM (r0 , t0 )), the gas density of a collapsing cloud is
2 0 (r0 ) r0 dr0 , (A4) 2 dr r where 0 (r0 ) is the initial gas density at r0 . The ratio of dr0 /dr determines how the thickness of a given shell evolves with time. The relative thickness dr0 /dr is determined by differentiating r = r0 cos2 with respect to r0 , yielding

(r, t) =

dr r = - r0 sin 2 arccos dr0 r0

r r0

d . dr0

(A5)

Next, d /dr0 is determined from an alternate form of equation (A3): + 0.5 sin 2 = t 2GM (r0 ) . 3 r0 dM (r0 ) 3 M ( r0 ) - dr0 r0 (A6)

Differentiating with respect to r0 yields d = dr0 G 2 M (r0 )r
3 0

.

(A7)

Now that the density (r, t) and velocity v (r, t) distributions of a collapsing pressure-free sphere are explicitly determined, the mass accretion rate at any given radial distance r and time t can be found as M (r, t) = 4 r 2 (r, t)v (r, t). A2 Collapse with non-zero initial velo city

APPENDIX A: PRESSURE-FREE COLLAPSE A1 Collapse from rest

In a general case of non-zero initial radial velocity profile v0 (r0 ), integration of eqaution (A1) yields v= dr =- dt 2GM (r0 ) 1 1 2 - + v0 (r0 ), r r0 (A8)

The equation of motion of a pressure-free, self-gravitating spherically symmetric cloud is GM (r) dv =- , r2 dt where v distance inside a to yield v= is the velocity of r from the center sphere of radius r. the expression for 2GM (r0 ) (A1) a thin spherical shell at a radial of a cloud, and M (r) is the mass equation (A1) can be integrated velocity v (A2)

where v0 (r0 ) is the initial velocity of a shell at r0 . Equation (A8) can be reduced to an integrable one by means of a substitution r = r0 cos2 sin(2 ) d tan + a
2

= dt

2GM (r0 ) , 3 r0

(A9)

dr =- dt

1 1 , - r r0

where r0 is the initial position of a mass shell at t = 0, and M (r0 ) is the mass inside r0 . Here, it is assumed that all shells are initially at rest: v0 (r0 ) = 0. Equation (A2) can be integrated by means of the substitution r = r0 cos2 (see Hunter 1962) to determine the time it takes for a shell located initially at r0 to move to a smaller radial distance r due to the gravitational pull of the mass M (r0 ). The answer is

where a = r0 v0 (r0 )/2GM (r0 ). Another substitution sin2 = x and integration over x from x = 0 to x = sin2 finally gives the time t it would take for a shell located initially at r0 and having a non-zero initial velocity v0 to move to a smaller radial distance r: t = (1 - a)2GM (r0 ) 3 r0
-1 -1 2

1-

- - (1 + ) tan +(1 + ) tan
-1



r + r0

r r0

, (A10)

1 - r/r0 + r/r
0

c 0000 RAS, MNRAS 000, 000­000


Finite Mass Reservoir Col lapse

11

Figure A1. The mass accretion rate as a function of time and radial distance from the center of a pressure-free cloud that has initial radial gas density distribution = c /[1+(r/rc )2 ], in which c = 7.5 â 104 cm-3 and rc = 0.033 pc.

Figure A2. The same as Fig. A1 but for the initial gas density distribution = c /[1 + (r/rc )3 ].

A3.2

A steeper profile

where = a/(1 - a). In the case of a non-zero initial velocity profile, it is more complicated to obtain a simple analogue to equations (A5)-(A7) and explicitly determine a density distribution (r, t), as done in the previous example. Instead, we obtain the mass accretion rate by computing the mass that passes the sphere of radius r during time t, i.e. M ( r0 + M (r, t) = r 0 ) - M (r0 ) , t (A11)

where M (r0 ) is the mass inside a sphere of radius r0 . A time interval t is the time that it takes for two adjacent shells of radius r0 and r0 + r0 to move to the radial distance r. The value of t can be found by solving equation (A10) for fixed values of r0 , r0 + r0 , and r. A3 Applications

The submillimeter and mid-infrared observations of WardThompson et al. (1999) and Bacmann et al. (2000) suggest that the gas density in the envelope of a starless core falls off steeper than r-2 . As a second example, we consider a pressure-free cloud with the initial gas density profile = c /[1 + (r/rc )3 ] and plot the corresponding mass accre tion rate M (r, t) in Fig. A1. The values of c and rc are retained from the previous example. As is seen, the tempo ral evolution of M strongly depends on the radial distance 4 r. At r 10 AU, the mass accretion rate has a well-defined maximum at t 0.21 Myr, when the central gas density has exceeded 1010 cm-3 (the central stellar core formation). Af ter stellar core formation, M drops as t-1 . At r > 104 AU, the temporal evolution of M does not show a well-defined maximum.

As two examples, we consider two different initial gas density profiles and determine the pressure-free mass accretion rate M (r, t) as a function of radial distance r and time t. A3.1 Modified isothermal sphere

First, we consider the radial gas density profile of a modified isothermal sphere: = c /[1 + (r/rc )2 ] (Binney & Tremaine 1987), where c is the gas density in the center of a cloud and rc is the radial scale length. Figure A1 shows M (r, t) of a pressure-free cloud with c = 5.5 â 104 cm-3 and rc = 0.033 pc. The mass accretion rate increases with time and appears to approach a constant value at later times t > 0.7 Myr. Note that the temporal evolution of the mass accretion rate depends on the radial distance r: M approaches faster a constant value at smaller r. This behav ior of M (r, t) is independent of the adopted values of c and rc .
c 0000 RAS, MNRAS 000, 000­000