A&A 398, 1081-1090 (2003)
DOI: 10.1051/0004-6361:20021707

From clouds to stars

Protostellar collapse and the evolution to the pre-main sequence I. Equations and evolution in the Hertzsprung-Russell diagram

G. Wuchterl^1,2 - W. M. Tscharnuter^2,3

1 - Max-Planck-Institut für extraterrestrische Physik, Giessenbachstraße, 85748 Garching, Germany
2 - Institut für Theoretische Astrophysik der Universität Heidelberg, Tiergartenstraße 15, 69121 Heidelberg, Germany
3 - Interdisziplinäres Zentrum für Wissenschaftliches Rechnen der Universität Heidelberg, Im Neuenheimer Feld 368, 69120 Heidelberg, Germany

Received 16 June 1999 / Accepted 8 November 2002

Abstract
We present the first study of early stellar evolution with ``cloud'' initial conditions utilizing a system of equations that comprises a solar model solution. All previous studies of protostellar collapse either make numerous assumptions specifically tailored for different parts of the flow and different parts of the evolution or they do not reach the pre-main sequence phase. We calculate the pre-main sequence properties of marginally gravitationally unstable, isothermal, equilibrium "Bonnor-Ebert'' spheres with an initial temperature of $10 ~ \rm K$ and masses of 0.05 to 10 $M_\odot$. The mass accretion rate is determined by the solution of the flow equations rather than being prescribed or neglected. In our study we determine the protostar's radii and the thermal structure together with the mass and mass accretion rate, luminosity and effective temperature during its evolution to a stellar pre-main sequence object. We calculate the time needed to accrete the final stellar masses, the corresponding mean mass accretion rates and median luminosities, and the corresponding evolutionary tracks in the theoretical Hertzsprung-Russell diagram. We derive these quantities from the gas flow resulting from cloud collapse. We do not assume a value for an "initial'' stellar radius and an "initial'' stellar thermal structure at the "top of the track'', the Hayashi-line or any other instant of the evolution. Instead we solve the flow equations for a cloud fragment with spherical symmetry. The system of equations we use contains the equations of stellar structure and evolution as a limiting case and has been tested by a standard solar model and by classical stellar pulsations (Wuchterl & Feuchtinger 1998; Feuchtinger 1999; Dorfi & Feuchtinger 1999). When dynamical accretion effects have become sufficiently small so that a comparison to existing hydrostatic stellar evolution calculations for corresponding masses can be made, young stars of $2~ M_\odot$ appear close to the location of the Henyey part of the respective classical evolutionary track and at substantially larger ages for given luminosities than those inferred from previous calculations. $1~{M_\odot}$ stars appear at lower luminosities, to the left of the corresponding Hayashi-tracks and are about $1 ~ \rm Myr$ older than an a-priori hydrostatic stellar evolution model at the same luminosity. They burn most of their deuterium during the main accretion phase before mass accretion halts and they become visible. They do not become fully convective during the entire evolution calculated, i.e., up to 1.5 Myr. Altogether the structure of our solar mass star at 1 Myr, with its raditive core and convective envelope, resembles the present Sun rather then a fully convective object. Very low mass stars and proto brown dwarfs close to the substellar limit appear with luminosities close to those at the "top of the tracks'', giving ages roughly in accordance with classical values, tentatively at 0.05 to $0.09~\rm dex$ higher effective temperatures.

Key words: stars: formation - stars: pre-main sequence - evolution - hydrodynamics - convection

1 Introduction

One important aim of star formation theory, which is still an unsolved problem, is to provide the initial conditions for stellar evolution, i.e., the masses, radii and the internal structure of young stars as soon as they can be considered to be in hydrostatic equilibrium for the first time.

Figure 1: Early stellar evolution in the Hertzsprung-Russell diagram: Comparison of results for the collapse of a $1~{M_\odot}$, "Bonnor-Ebert'' sphere (thick line) and a classical pre-main sequence track ( thin line). Both are for mixing length convection, with $\alpha _{\rm ML}=1.4$. The classical PMS-track (thin line on the left) is the $1~{M_\odot}$, "MLT Alexander'', case after D'Antona & Mazzitelli (1994). The thick line for the cloud collapse gives the effective temperature as defined by $L = 4 \pi r_\tau ^2 \sigma T_{\rm eff}^4$. $r_\tau$ is the radius at which the Rosseland mean optical depth, $\tau _{\rm Ross}$ as integrated inward from the outer boundary of the cloud equals 2/3. Alternative temperatures are given for the cloud collapse results: the dashed line is the temperature defined by $L = 4 \pi r^2 \sigma T^4$, at the first point encountered inward starting from the outer boundary of the cloud, i.e., based on the temperature radius as discussed by Baschek et al. (1991). The dotted line gives the temperature at $\tau _{\rm Ross} = 2/3$, an estimate for the photospheric temperature. Note that the photospheres obtained in the collapse calculations are convective at this point and $T=T_{\rm eff}$ is not guaranteed. The two diamonds indicate zero age, defined here as the moment when $\tau _{\rm Ross} = 2/3$ for the first time. Temperatures before zero age (on the dash dotted curve) are central temperatures of the still transparent cloud. The alternative temperatures are given here to illustrate the uncertainties involved in putting a collapsing cloud fragment into the HRD and to emphasise different physical aspects.

Open with DEXTER

Following the evolution of a collapsing cloud fragment requires self-gravity, a description of the (supersonic, compressible) motion of the cloud gas, radiative transfer in extended, non-planeparallel media of low, moderate and very high optical depth, and, sooner or later in the non-isothermal phase, a time-dependent treatment of convection.

We account for all of the mentioned processes that are of importance in clouds and stars. However, we assume spherical symmetry and the grey Eddington approximation for radiative transfer to keep the problem as simple as possible, since this is the first look at the problem of cloud collapse, star formation and early stellar evolution that is based on the solution of one set of equations for all evolutionary stages. In the next section we assemble this system of equations that contains all the physics described above.

To discuss the results and illustrate the differences we focus on the comparison with classical, i.e., fully hydrostatic results that are obtained with model equations (stellar evolution equations) being as close as possible to the "stellar structure limit'' of our set of equations. Comparison of our results with studies that include more physical processes (e.g., disc accretion or frequency dependent photospheric radiative transfer) can then be made by using existing intercomparisons of different hydrostatic studies to our hydrostatic reference study.

To relate our results to quasi-hydrostatic stellar evolution calculations on the pre-main sequence we chose the calculations by D'Antona & Mazzitelli (1994), DM94, as our reference for the following reasons: (1), the photospheric radiative transfer in DM94 is an almost accurate limit of our grey radiative transfer in the Eddington approximation in its "stellar structure limit''. We thereby exclude differences between our calculations and the quasi-hydrostatic calculations which might be due to effects that are specific to the details of the non-grey treatment of the photospheres. (2), the low temperature, molecular opacities for DM94, "Alexander'' sets of PMS-tracks are the same as used here. (3), standard mixing length convection, our "stellar structure limit'' for convection, is one of the convection theories used by DM94.

2 A model for protostellar collapse and the pre-main sequence

2.1 Flow equations

We use the equations of radiation hydrodynamics (RHD) in the grey Eddington approximation and with spherical symmetry (Castor 1972; Mihalas & Mihalas 1984) in their integral form (Winkler & Norman 1986). Convective energy transfer and mixing is included by using a new, time-dependent convection model derived from the model of Kuhfuß (1987). It closely approximates standard mixing length theory in a local, static limit and has been successfully tested in calculations of the Sun and RR Lyrae pulsations (Wuchterl & Feuchtinger 1998; Feuchtinger 1999a,b). The convection model is used with the standard parameters given in Wuchterl & Feuchtinger (1998). The few cases where the mixing length parameter $\alpha _{\rm ML}$ is changed are explicitly emphasised. In all other cases the standard value $\alpha _{\rm ML} = 1.5$ is used, which results in a solar model that is in accord with standard solar models and the Sun for our "stellar structure limit''. The equations and boundary conditions are summarized in Table .1, in the Appendix. Detailed equations of state (Wuchterl 1989, 1990) and opacities for a (proto) solar composition are used (X = 0.77, Y = 0.213, Z = 0.017, $X_{\rm D} = 3.85\times 10^{-5}$). The equations of state have been compared to the MHD equations of state (Mihalas et al. 1988) and agree to better than the table-interpolation errors (cf. Götz 1993). Opacities include the contribution of dust particles with "interstellar'' properties (Yorke 1979, 1980a), Alexander & Ferguson (1994) values in the molecular range and Weiss et al. (1990) Los Alamos high temperature opacities (compilation by Götz 1993). Convective mixing of deuterium and deuterium burning (with reaction rates from Caughlan & Fowler 1988) are included as outlined in the next section.

2.2 Deuterium burning

During early stellar evolution the first nuclear reaction that contributes significantly to the energy budget of a young star is the D burning reaction $\rm {}^2H(p,\gamma){}^3He$ (e.g., Kippenhahn & Weigert 1990). The energy release of $5.5 \rm ~ MeV$ per reaction contributes as a source $\epsilon_{\rm nuc}^{\rm D}$ to our internal energy equation (cf. Table .1). The proton partial density, $\varrho _{\rm p}$, is approximated by the hydrogen partial density, $X \rm\varrho$. Once D becomes significantly depleted, it has to be budgeted by the D balance equation that accounts for D destruction and transfer from regions of different D abundance due to advection and the D flux due to convective mixing, $j_{\rm D}$, (see, Table .1, Eq. (3)). Rate coefficients, $\left<\sigma v\right>$ for $\rm {}^2H(p,\gamma){}^3He$ are taken from Caughlan & Fowler (1988). D burning via $\rm {}^2H(p,\gamma){}^3He$ is the only nuclear reaction that we take into account. As a consequence, our study is limited to evolutionary phases before the onset of hydrogen burning and does not address the evolution of abundances of trace species like $\rm {}^7Li$ and $\rm Be$.

It is assured that our extensions to the system of stellar structure equations do have a fixed relation to standard stellar structure and stellar hydrodynamics. After requiring a correct solar convection zone in the "stellar structure limit'', the equations we use do not contain any additional free parameters and the formation and evolution of a star is determined by the properties of the initial cloud fragment alone. The star formation theory is fully specified with the specification of the initial cloud state and the boundary conditions. See Appendix A for the equations.

3 Collapse of Bonnor-Ebert spheres

Figure 2: The relation between age and luminosity for the collapse of Bonnor-Ebert spheres ( full lines) and classical, fully hydrostatic models of pre-main sequence evolution. These relations are relevant for the determination of PMS-ages. For the comparison we use the quasi-hydrostatic results by D'Antona & Mazzitelli (1994) ("Alexander + MLT'' case) for 0.1, 0.5, 1 and $2~ M_\odot$ ( dotted lines from left to right). Classical PMS-models do not accrete and evolve at constant mass. The new PMS-ages calculated from the collapse of cloud fragments are larger than the classical values by up to a Myr, depending on luminosity. This is due to an age offset because the formation process is modeled and due to a different PMS-contraction behaviour caused by the difference between the assumed initial thermal structure in classical tracks and the thermal structure calculated from the collapse and accretion process. Note that collapse luminosities show an initial rise, close to zero age, as the protostars are "switched'' on after core formation. Ages inferred using collapse PMS-tracks will be larger than those inferred from classical tracks by up to a Myr at ages below 2 Myr. The dashed and dash-dotted line, close to the $1~{M_\odot}$-relation, are obtained by changing the mixing length parameter $\alpha _{\rm ML}$ to 1.2 and 1.4, respectively, for the same fragment mass, see text.

Open with DEXTER

In this section we discuss the results for the collapse of critical Bonnor-Ebert spheres with a temperature of $10 ~ \rm K$ and masses of 0.05 to $10~{M_\odot}$ within a constant volume, as obtained by the methods described above. For a description of the method of solution we refer to Appendix B.

3.1 Structure of results and definitions

The results described below are quantities derived from the flows which are obtained as the solution of equations given in Table .1, for the Bonnor-Ebert collapses. Since there is no a priori hydrostatic part and no a priori photospheric model in our model (they result as special properties of the solutions of our flow equations) we describe how the stellar parameters age, effective temperature, radii and masses are obtained from the radiation hydrodynamical flows. We restrict our discussion mostly to these basic parameters to give an overview over the collapse solutions. We start our discussion with a proposal for the definition of zero stellar age, because as it turns out that most our definitions of stellar parameters will only be valid after this instant of time.

3.2 Stellar zero age

The modeling of the star formation process together with the PMS-contraction opens up the possibility to define a stellar zero age. We require the instant of age zero to be (1), well defined, (2), a unique clock for an individual star; it should be reset by a destruction of the star and stellar mergers, and (3), the gravitationally weakly bound, tenuous states of the cloud should not contribute to stellar age - a long-lasting quasi-static cloud phase would otherwise induce almost arbitrary age additions without any consequences for the subsequent evolution. To compare the luminosities of objects of various masses at a given age a general definition of zero age is necessary. As the thermally controlled Kelvin-Helmholtz contraction is the dominating evolutionary process in young gas spheres, whether they are giant planets, brown dwarfs or stars, we propose to use as age zero the instant of time, when the interior of the gas spheres is thermally enclosed for the first time. The enclosure means that photons cannot escape from the entire object directly but are radiated from a photosphere. Energy transfer to the photospheric bottleneck then delays the cooling and determines how the thermal reservoir in the interior is emptied. The cooling history can then be used to define an age since the heat reservoir was formed. As a practical definition for age zero we propose to use the instant of time when the Rosseland mean optical depth of a gaseous object equals 2/3, i.e.,

$$\mbox{Age} = 0 \!:\ \tau_{\rm Ross} (t) = \int^R_{r^\prime} \kappa_{\rm Ross} \varrho ~ {\rm d}r = \frac{2}{3} \cdot$$ (1)

We use 2/3 here because of the advantage for, e.g., the $T_{\rm eff}$ definition in the comparison to D'Antona & Mazzitelli (1994) but recommend to use $\mbox{Age} = 0$ at $\tau_{\rm Ross} (t) = 1$, in general. $\tau=2/3$ results from an averaging of optical depth over the stellar disk that is specific to the grey Eddington approximation.

For a protostar this instant of time is practically at the end of the isothermal phase of the collapse. Furthermore, with definition (1) the luminosities of the protoplanets can be compared to stars and brown dwarfs of similar age. The proposed definition of zero age has the following properties: (1), it preceeds the short phase of the "final hydrostatic core formation'' for protostars that was used in earlier work (Appenzeller & Tscharnuter 1975; Winkler & Newman 1980a) by only $\sim$ $2000 ~{\rm yr}$. It corresponds to the sharp initial luminosity rise (within 1000 yr for 0.05 to 10 $M_\odot$) that occurs when the first (molecular hydrogen) protostellar cores form. Because this happens during this rapid dynamical density enhancement due to early cloud collapse, the proposed definition is very sharp, in the sense that the properties of the cloud vary rapidly in the vicinity of zero age; $\mbox{Age} = 0$ is reached during a phase of rapid changes in the observables. As a consequence it will be observationally well defined; (2), it is "close'' to the formation time of the oldest "rocks'' in the solar system. They presumably form during the process of planetesimal formation that occurs according to present understanding within $\sim$ $10^4 ~{\rm yr}$ after pre-planetary nebula formation. The pre-planetary nebula forms roughly at the time when the final hydrostatic protostellar core settles into hydrostatic equilibrium (Tscharnuter 1987). Hence, together with the $\sim$ $2000 ~{\rm yr}$ for final core formation, the ages of the oldest meteorites and meteoritic components that are determined by absolute radioactive dating should be only $\sim$ $10^4 ~{\rm yr}$ smaller than the ages defined here. This results in a conceptual consistency of our ages and the meteoritic ages to within probably < $10^5 ~{\rm yr}$. It is a useful property because solar system ages of $4.566\pm 0.002$ Gyr (Allègre et al. 1995; for Allende CAIs) are used to calibrate stellar structure theory by assuming that the age of the Sun is equal to the age of the meteorites and a comparison of stellar-evolution models of the Sun at theoretical ages with properties of the present Sun; (3), relative ages of different objects in a star formation region are physically well defined; (4), the thermal clock which is based on the dominating dynamical process in early stellar evolution, namely gravitational contraction and cooling is (re)set at the proposed zero age; (5), the age spread in young clusters is conceptually accessible without reference to the ZAMS, and knowledge about it. Concepts involving the ZAMS rely on the assumption that it will be reached. This is not guaranteed as it does not exist for brown dwarfs and planets. Hence such concepts would be inapplicable for those objects. Our proposed zero age is a well defined concept for all masses at least from an earth mass up to the largest stellar masses; (6), questions about coevolution of very young binaries can be addressed because individual component ages are well defined; (7), with this zero age we use the earliest instant when an effective temperature can be defined, based on a radius defined by using optical depth. With age definition (1) we can study the luminosity emitted by our cloud volume as a function of age. In the beginning the luminosity comes from the still almost isothermal cloud fragment, in the end from the pre-main sequence star that has formed in the centre. We will use this age to discuss the effect of the cloud collapse on pre-main sequence ages. We can expect two main effects relative to hydrostatic models of early stellar evolution: (i) an age offset that is due to the fact that the collapse phase which is not modeled in fully hydrostatic PMS-studies adds contributions of the order of a free-fall time to the ages usually given (see Tables 1, 2), and (ii) an effect due to the different PMS-contraction behaviour that is due to the difference between the assumed initial thermal structure of classical PMS-calculations and the thermal structure as calculated from the collapse and accretion process in this work.

In our discussion, we use the relation between luminosity and age for two reasons: (1) because of the uncertainties (both theoretical and observational) in the determination of the effective temperature of objects in the vicinity of the Hayashi-tracks, and (2) because these tracks are very steep in the HRD during the early PMS evolution - the luminosity decreases at almost constant effective temperature. These relations for the collapse of Bonnor-Ebert spheres of 0.1, 0.5, 1 and $2~ M_\odot$ are depicted in Fig. 2. The relations determined from collapse (full lines) originate from zero luminosity, i.e., before the cloud fragment starts to radiate significantly. The fully hydrostatic relations (dotted lines) begin at an assumed initial age and at a luminosity that is the consequence of the assumed (usually fully convective, or n=3/2 polytrope) initial thermal structure.

Notice the age offset of about $10^5 ~{\rm yr}$, about a free-fall time, between the beginning of the classical and collapse curves in all cases. This is due to the assumed initial ages of fully hydrostatic calculations. Next we look at the possibility of attributing all the differences to an age offset. This could be done by shifting the classical curves upward to higher ages until the luminosities match the new curves. Such a procedure could be considered successful if the rates of change in luminosity of classical and collapse L(t)-relations match, too. This would be indicated by the luminosity difference between corresponding old and new relations of the same mass to be constant. For the $0.1~{M_\odot}$ case this is fairly successful and the collapse may be essentially viewed as producing an age offset of about $300~ {\rm kyr}$. For the other masses this procedure does not hold. Even for the $0.5 ~{M_\odot}$ relations an addition of $400~ {\rm kyr}$ to the classical relation might just bring the luminosities in accord, but the relations still diverge because of the different contraction behaviour.

The consequences for the inferred ages at given luminosity are already considerable for $0.5 ~{M_\odot}$: at the end of the collapse calculation, a PMS-object with solar luminosity has a new age of $1.5 ~{\rm Myr}$ while following down the $1~ L_\odot$-ordinate to the classical tracks gives $0.5~\rm Myr$. The situation is very similar for the 1 and $2~ M_\odot$ clouds, but the differences in luminosity evolution to the classical relations are even more pronounced. Changing the mixing length, a procedure that probably produces the largest changes in classical models does not significantly change this conclusion. This can be seen by comparing the dashed and dash-dotted lines, for $\alpha_{\rm ML}=1.2$ and 1.4, respectively, that are very close to the standard 1 $M_\odot$ collapse relation for $\alpha _{\rm ML} = 1.5$.

Next the adopted age definition allows us to define the "stellar'' properties, i.e., the properties of the optically thick, ultimately hydrostatic parts of the flow for any instant after zero age. This is described in the next section.

3.3 Radii

Here we focus on the observables of the collapsing cloud fragments and the resulting young stars, so we define our stellar radii to be the optical depth radii, $R_\tau$, obtained by determining the radius r' at which

$$\tau = \int^R_{r^\prime} \kappa_{\rm Ross} \varrho ~ {\rm d}r = 2/3 .$$ (2)

Because of our definition of zero age such a radius exists. This is not necessarily close to the density radius, $R_\varrho$, i.e., the location of the maximum absolute value of the density gradient, that usually defines the stellar radius (density gradients in shocks are ignored for this procedure cf. Winkler & Newman 1980a,b; Baschek et al. 1991). Note that Balluch (1991a-c) used the radius of the accretion shock as an estimate for $R_\varrho$ and the stellar radius, which is a good estimate due to the small stand-off distance of the shock.

The optical depth radius gives a measure of the radius from which the photons escape and which parts of the flow are thermally shielded from the outside. It is the length scale a (bolometric) observer would see. It is very close to the density radius once the protostellar cocoon has been sufficiently diluted by accretion to make it transparent down to the stellar photosphere. The two radii converge when $R_\tau$ "jumps'' from the dust photosphere to the stellar photosphere and intercepts $R_\varrho$.

3.4 Luminosity and effective temperature

With the luminosity at the outer boundary and the optical depth radius we determine an effective temperature by

$$L = 4 \pi R_\tau^2 \sigma T_{\rm eff}^4 .$$ (3)

The so defined $T_{\rm eff}$-values are used in our Hertzsprung-Russell diagrams. To give an estimate for the uncertainty involved when placing a protostar into the Hertzsprung-Russell diagram we additionally plot two alternative temperature definitions in Fig. 1: (1) a temperature T(R_T)based on the temperature radius, R_T, that is governed by

$$r=R_T: L(r) = 4 \pi r^2 \sigma T(r)^4 ,$$ (4)

cf. Baschek et al. (1991). This is the temperature at a distinct point in the flow where the luminosity corresponds to the value that would be obtained by a blackbody radiating at the local flow matter-temperature. (2) a photospheric temperature defined by

$$T_{\rm phot} = T(R_\tau) ,$$ (5)

i.e., the local flow temperature at the optical depth radius. It is a measure of the temperature at the location from which the photons typically escape towards an observer.

In a static, compact, radiative photosphere in Eddington approximation, the definitions result in the same temperatures (cf. Baschek et al. 1991).

3.5 Collapse of a solar mass

The resulting theoretical Hertzsprung-Russell diagram for the collapse of a $1~{M_\odot}$ Bonnor-Ebert sphere is shown in Fig. 1. In this case the mixing length parameter $\alpha _{\rm ML}$ is chosen to be 1.4 to facilitate the comparison with the fully hydrostatic pre-main sequence calculation. The thick line shows the evolutionary track in the Hertzsprung-Russell diagram by using the effective temperature based on the optical depth radius $R_\tau$ (as used in earlier collapse studies by Appenzeller & Tscharnuter 1975 and Winkler & Newman 1980a). The track starts at the instant of zero age, indicated by the diamond symbol. As above noted the adopted definition of zero age coincides with the first instant at which the effective temperature concept can be used.

The initial decrease in $T_{\rm eff}$ is due to the fact that by definition $R_\tau$ starts out very small but rapidly increases during the first core formation. This occurs because the core's density and optical depth increase rapidly during this settling into hydrostatic equilibrium. $R_\tau$ increases at approximately constant luminosity and hence $T_{\rm eff}$ drops. Temperatures given before zero age (dash-dotted part) are central temperatures of the still transparent, almost isothermal cloud. With the onset of accretion onto the core the luminosity rises sharply with a moderate increase in $T_{\rm eff}$. The "mouse-head'' is the hall mark of the formation of the second (inner) atomic hydrogen core that forms after the (second) collapse of the dissociating molecular hydrogen core. At the cusp after the first luminosity maximum the final accretion phase is entered. As the luminosity rises towards the second maximum, the first core is accreted onto the second one and the Mach numbers of the flow increase above ten. The small luminosity decrease after the second maximum accompanies the adjustment of the accretion flow and the mass accretion rate onto the final core to the substantial luminosity now radiated from the accretion shock.

When the luminosity starts rising again the quasi-steady main-accretion phase commences. Shortly after the flat top part of the track is reached, the luminosity, due to D burning, becomes significant and reaches a maximum of 60% at the global luminosity maximum. When the volume is significantly depleted by the accretion process the luminosity starts decreasing and $T_{\rm eff}$ rises as parts of the flow closer to the young star become visible. When the dust cocoon becomes transparent the optical depth radius "jumps'' into the "stellar'' photosphere the effective temperature sharply increases, crosses the classical PMS-track, the Hayashi line and finally reaches the "corner'' of the collapse track, at maximum $T_{\rm eff}$. From this point on $T_{\rm eff}$ coincides with the temperature definition based on the temperature radius (thick and dashed line overlap).

Hence the radiation hydrodynamical results for temperature based on R_T and $R_\tau$ converge to the equality valid in classical, hydrostatic photospheres. The photospheric temperature (dotted) at $\tau=2/3$ passes the shock at this instant and produces a sharp spike to the left as $\tau=2/3$ is at the peak temperature just downstream from the accretion shock, and still upstream from the radiative relaxation layer. As the accretion luminosity becomes negligibly small and after "D burnout'' the collapse track turns almost parallel to the classical track, slightly pointing towards the ZAMS location of the Sun on the classical track. The photospheric temperature remains different from the other two temperatures, and slightly cooler than the $T_{\rm eff}$ value of the fully hydrostatic classical track. Note that the photosphere of the collapse calculation is convective all along the Hayashi-part of the track and $T_{\rm phot}=T_{\rm eff}$, valid for a compact, hydrostatic, radiative photosphere in Eddington approximation is not guaranteed under these conditions.

3.6 Collapse PMS-tracks and isochrones

Figure 3: Evolutionary tracks in the Hertzsprung-Russell diagram for the collapse and early pre-main sequence evolution of 0.05, 0.1, 0.5, 1, 2 and $10~{M_\odot}$ cloud fragments ( full lines). The dashed lines indicate isochrones for the collapse tracks, labeled with the respective ages. Zero age is defined here as the moment when the respective cloud fragment becomes optically thick and the interior is thermally locked as the first photosphere forms (Rosseland mean optical depth reaches 2/3). Quasi-hydrostatic tracks for 0.1, 0.5, 1 and $2~ M_\odot$ from D'Antona & Mazzitelli (1994), "Alexander MLT''-case, $\alpha _{\rm ML}=1.4$ ( dotted lines) are plotted for comparison.

Open with DEXTER

Keeping in mind the above discussion concerning temperatures we now proceed to a discussion of evolutionary tracks for cloud fragments of various masses and the respective isochrones in the Hertzsprung-Russell diagram, based on the effective temperature with the optical depth-radius definition (see Baschek et al. 1991 for a discussion, including the protostellar application). It is important to note that the isochrones are based on the chosen zero age and hence the assumption that the first photospheres of all stars have formed at the same time. Different collapse- and accretion times of fragments of different masses are included. This is in contrast to hydrostatic studies that count their ages from an arbitrary zero point related to there initial condition.

In Fig. 3 an overview over the collapse and PMS evolution of Bonnor-Ebert-sphere cloud fragments with masses of 0.05, 0.1, 0.5, 1, 2 and $10~{M_\odot}$ is given. They are compared to classical, non accreting, fully hydrostatic tracks for 0.1, 0.5, 1, and $2~ M_\odot$ after D'Antona & Mazzitelli (1994). For the comparison an electronic version of their "Alexander, MLT''-tracks for $\alpha _{\rm ML}=1.4$ was used.

The early evolution dominated by the initial cloud collapse is very similar for all masses. At an age of 1000 yr the first hydrostatic cores form, with an age spread of a few hundred years. All fragments enter the main accretion phase at an age slightly above 8 kyr. All tracks remain below effective temperatures of $2000~\rm K$ up to $30~\rm kyr$. The protostars still remain in their dust cocoons and radiate from dust photospheres. The luminosity decline has already started at this age for the lowest mass fragment. At $60~\rm kyr$ the first object has arrived on the pre-main sequence, close to the Hayashi track. It has formed from the lowest mass cloud fragment of $0.05~ M_{\odot}$ and is a young brown dwarf. The higher mass fragments follow until the $2~ M_\odot$ fragment arrives on its PMS-track at $0.5~\rm Myr$. The $10~{M_\odot}$ fragment is still accreting at $0.7~{\rm Myr}$, with a hydrostatic central stellar embryo of $4.5 ~ M_\odot$.

Figure 3 illustrates that star formation is a gravitational heating event that increases temperature by about two orders of magnitude. Due to the negative, so-called gravothermal specific heat of common, thermally supported self-gravitating systems the heating slows down when energy losses slow-down because matter becomes opaque due to the cross-sections of atoms. Ultimately nuclear energy is added which acts as a long-term thermostatic effect.

The early PMS evolution of the collapsed cloud fragments shows a relation to the classical tracks that changes with mass. The $2~ M_\odot$ fragment arrives close to the radiative, Henyey part of the corresponding classical track. The $1~{M_\odot}$ and the $0.5 ~{M_\odot}$ case arrive on the PMS with a luminosity corresponding to about half way down on the convective, Hayashi parts of the classical tracks. The $0.1~{M_\odot}$ fragment arrives with a luminosity corresponding to the top of the classical track. The $0.05~ M_{\odot}$ fragment (a proto brown dwarf) is shown to give an indication of the further mass dependence of the relation to the $0.1~{M_\odot}$ classical track. For the two lowest masses the age difference can be attributed to an offset in zero age. Furthermore, D burning has not begun for the $0.1~{M_\odot}$ case-as might be expected so close to the brown dwarf limit. It follows that, by using the appropriate initial thermal structure, this mass probably will show a quasi-classical PMS evolution starting near the top of a Hayashi track and a somewhat higher effective temperature, but, as found in previous studies, with the slow-down of contraction due to D burning.

On the other hand the results of the collapse calculations indicate substantial corrections in PMS stellar properties for masses of $0.5 ~{M_\odot}$ and higher. PMS objects of these masses already are half-way in the direction of the Henyey-track when they become visible.

3.7 Lines of constant mass-isopleths

Due to our radiation hydrodynamical approach only the masses of the cloud fragments are prescribed. The stellar mass and the corresponding mass accretion rate are calculated from the flow. Hence, quantities like stellar mass and mass accretion rate are determined during the evolution of the flow.

For our discussion it is useful to use the optically thick mass, $M_\tau$, for further considerations. It follows naturally by determining the mass inside the optical depth radius, $R_\tau$ that we have already used in the definition of $T_{\rm eff}$. This corresponds to the hydrostatic mass, defined, e.g., via the density radius (cf. Baschek et al. 1991) for most of the time but has to be distinguished during the very early accretion when calculating the accretion-luminosity from $\dot{M_*} GM_*/R_*$. We use $M_\tau$ for the remainder.

In order to show the stellar properties as accretion proceeds we display lines of constant optically thick mass, $M_\tau$ - which we may call isopleths - in the Hertzsprung-Russell diagram in Fig. 4. It is the actual stellar mass that determines the stellar properties and not the fragment mass that includes the entire circumstellar environment under consideration. The $10~{M_\odot}$ fragment, for example, starts as the other fragments with $M_\tau=0.04~{M_\odot}$ and produces the corresponding luminosities and effective temperatures despite of its $10~{M_\odot}$protostellar envelope.

Figure 4: Optically thick mass, $M_\tau$, in the Hertzsprung-Russell diagram for the collapse and early pre-main sequence evolution of 0.05, 0.1, 0.5, 1, 2 and $10~{M_\odot}$ cloud fragments. $M_\tau$ is the mass interior to the optical depth radius, $R_\tau$. Dash-dotted lines are locations of constant $M_\tau$, isopleths, with values labelled along the $10~{M_\odot}$-track. The locations where $M_\tau$ corresponds to 80, 90 and 99% of the final mass are connected by the thin lines. The thin line for 99% is used as an estimate for the end of significant mass accretion. The lines of constant mass are essentially tangential to the tracks beyond this point. Quasi-hydrostatic tracks for 0.1, 0.5, 1 and $2~ M_\odot$. (D'Antona & Mazzitelli 1994, "Alexander MLT''-case, $\alpha _{\rm ML}=1.4$) are shown for comparison.

Open with DEXTER

The lowest isopleth at $0.01~ {M_\odot}$ illustrates the fact that hydrostatic structures of that mass become optically thick (opacity limit of fragmentation) and, as a consequence, all initial protostellar cores typically form with a few percent of a solar mass. At $M_\tau=0.04~{M_\odot}$ the lowest mass fragment has almost terminated its main accretion phase while the higher mass fragments still grow their cores and accelerate their accretion flows. As $M_\tau$ increases to the final mass the isopleths become tangential to the evolutionary tracks and at constant mass the two curves are identical. Note the decrease in luminosity when the track of a fragment approaches an isopleth close to its final mass. The beginning exhaustion of the mass reservoir becomes then noticeable as a drop in accretion luminosity. To mark the end of significant mass accretion, the point at which 99% of the final mass has been accreted is indicated by a thin line. For the two smallest masses this happens well before they reach the Hayashi phase. Accretion, as a consequence, does not play a role during their PMS-contraction as calculated here. However accretion does play a role in defining their thermal structure at the "top of the track''. For masses from $0.5 ~{M_\odot}$ upward the end of significant mass accretion is terminated after the Hayashi line has been crossed.

4 Global results

To facilitate statistical comparison with observations of star formation and to have values for the factors of order unity in scaling arguments we give some quantities averaged over the star formation time for each fragment in Tables 1 and 2. Ages at $M_\tau(t)=0.99~ M_{\rm final}$, $t_{\rm0.99M}$ are given in units of Myr and the initial free-fall time; mean mass-accretion rates ${<\dot{M}>} = M_\tau(t) / {\rm Age}(t)$ are given in units of $M_\odot /\rm yr$ and $M_{\rm final} / t_{\rm ff}$; M₀ is the fragment mass. For the luminosities the median values <L> are given in solar units.

 

 

Table 1: Duration of accretion and averages of accretion rates and luminosities for the accretion time lasting from zero age to the end of significant mass accretion. For the end of significant mass accretion we use the moment when the young star has reached an optically thick mass $M_\tau$ of 99% of its final mass. The different $1~{M_\odot}$ cases are for different mixing length parameters.

Accretion time duration and averages
M0	$\rm t_{0.99 M}$		$<\dot{M}>$		<L>	Remark
$M_\odot$	$\rm Myr$	$t_{\rm ff}$	$M_\odot /\rm yr$	$M/t_{\rm ff}$	$L_\odot$
1	0.55	3.23	$1.79 \times 10^{-6}$	0.30	13.1	RC1M
1	0.60	3.54	$1.65 \times 10^{-6}$	0.28	14.1	RCML15
1	0.55	3.26	$1.79 \times 10^{-6}$	0.30	12.8	RCML14 $\alpha _{\rm ML}=1.4$
1	0.54	3.18	$1.84 \times 10^{-6}$	0.31	11.9	RCML12 $\alpha_{\rm ML}=1.2$
0.05	0.038	4.51	$1.30 \times 10^{-6}$	0.22	1.02	RC005M
0.1	0.074	4.34	$1.35 \times 10^{-6}$	0.23	2.21	RC01M
0.5	0.32	3.81	$1.56 \times 10^{-6}$	0.26	7.58	RC05M
1	0.55	3.23	$1.79 \times 10^{-6}$	0.30	13.0	RC1M
2	1.28	3.64	$1.54 \times 10^{-6}$	0.27	25.1	RC2M

Quantities are averaged over a time starting at zero age and ending at the end of significant mass-accretion. The end of significant mass accretion is defined as the instant at which 99% of the final mass has been accreted, more precisely when $M_\tau = 0.99~ M_{\rm final}$, cf. the respective thin line in Fig. 4. To indicate the sensitivity on the definition of the end of accretion we give a few averages obtained for the $1~{M_\odot}$ case for periods ending with $M_\tau$ closer to the final mass in Table 2.

 

 

Table 2: Duration of accretion and averages of accretion rates and luminosities for the accretion time lasting from zero age to the end of significant mass accretion. For the end of significant mass accretion we use different fractions $M_\tau /M_{\rm final}$ to show how averages change for the $1~{M_\odot}$ case. The mean mass accretion rate for the $10~{M_\odot}$ fragment corresponds to $0.89 M/t_{\rm ff}$ in the first 0.7 Myr ($\sim$ $0.4 ~ t_{\rm ff}$).

Accretion time and averages for different end of accretion
and a $10~{M_\odot}$ fragment.
$M_\tau$	$\rm Age$	$<\dot{M}>$	M0	$M_{\rm final}$	Remark
$M_\odot$	$\rm Myr$	$M_\odot$/yr	$M_\odot$	$M_\odot$
0.99	0.60	$1.67 \times 10^{-6}$	1	1	RCML15
0.99	0.55	$1.80 \times 10^{-6}$	1	1	RC1M
0.999	0.82	$1.22 \times 10^{-6}$	1	1	RC1M
0.9999	1.08	$9.26 \times 10^{-7}$	1	1	RC1M
0.99999	1.29	$7.78 \times 10^{-7}$	1	1	RC1M
	0.70	$7.11 \times 10^{-7}$	10	4.46	RC10M

To indicate the effects due to the time-dependence of the accretion flow and, in particular, due to the temporal variation of $\dot{M}$, we split the star formation time into 4 phases, based on $M_\tau /M_{\rm final}$. By this procedure we obtain 4 different populations corresponding to different stages of the accretion process. For each of the 4 phases we give averaged quantities, for fragment masses, M₀, of 0.05, 0.1, 0.5, 1 and $2~ M_\odot$, in Table 3.

5 Conclusions

A study with spherical symmetry can certainly not address all aspects of the star formation process, as magnetic fields and angular momentum are ignored, but it has the advantage that the process can be studied by keeping close to basic physical principles-at least as close as stellar evolution theory, with the notorious problem of the convection model.

Based on the concept of "maximum deduction at minimum assumptions'', spherical symmetry provides a useful limiting and reference case. With the symmetry assumption we are able to keep a close and explicit relation to the fundamental conservation laws and other basic physical principles. The essential modeling assumption that goes into this work then is the ever present problem of the convection model. But we have tried not to add further ambiguity that is not already present in the modeling of stellar evolution. As a consequence of our convection model the system of equations we use represents the structure of the Sun correctly. The number of free parameters in any study of star formation cannot be made smaller than in ours unless a parameter-free description of convection becomes available.

The energetics of the star formation process can be addressed rather completely, most notably because all relevant energy transfer processes, i.e., radiation, convection, advection and conduction can be described accurately. Very likely star formation is most efficient with zero angular momentum and hence a lower limit for star formation time and upper limits for luminositites and accretion-rates are obtained.

 

 

Table 3: Duration and averages of accretion rates and luminosities for different accretion phases of cloud fragments of 0.05, 0.1, 0.5, 1 and $2~ M_\odot$. The phases are defined by ranges in optical thick mass $M_\tau$ in units of the final mass $M_{\rm final}$, in the second column. For each phase and fragment mass, M₀, duration, $\Delta t = \min({\rm Age})-\max({\rm Age})$, mean mass-accretion rates, ${<\dot{M}>} = (\min(M_\tau)-\max(M_\tau))/\Delta t$, and the median value for the luminosities, <L>, is given.

Averages for 4 accretion phases separated by $M_\tau/M_{\rm final}=0,0.8,0.9,0.99$ and 0.999 (cf. Fig. 4).

M0

0.05 $M_\odot$

0.1 $M_\odot$

0.5 $M_\odot$

1 $M_\odot$

2 $M_\odot$

$M_\tau$

$\Delta t$

$<\dot{M}>$

<L>

$\Delta t$

$<\dot{M}>$

<L>

$\Delta t$

$<\dot{M}>$

<L>

$\Delta t$

$<\dot{M}>$

<L>

$\Delta t$

$<\dot{M}>$

<L>

$[M_{\rm final}]$

$\rm kyr$

$M_\odot$/yr

$L_\odot$

$\rm kyr$

$M_\odot$/yr

$L_\odot$

$\rm kyr$

$M_\odot$/yr

$L_\odot$

$\rm kyr$

$M_\odot$/yr

$L_\odot$

$\rm kyr$

$M_\odot$/yr

$L_\odot$

$\rm0$

$\le$0.8

16

$2.5 \times 10^{-6}$

1.28

30

$2.7\times 10^{-6}$

2.92

126

$3.2\times 10^{-6}$

12.4

189

$4.2 \times 10^{-6}$

19.5

496

$3.2\times 10^{-6}$

34.1

$\rm I$

(0.8,0.9]

6

$8.1 \times 10^{-7}$

1.19

11

$8.8\times 10^{-7}$

2.63

52

$9.8\times 10^{-7}$

8.21

85

$1.2 \times 10^{-6}$

15.3

200

$1.0\times 10^{-6}$

27.3

$\rm II$

(0.9,0.99]

16

$2.7 \times 10^{-7}$

0.51

32

$2.7 \times 10^{-7}$

0.97

144

$3.0\times 10^{-7}$

4.60

280

$3.2 \times 10^{-7}$

8.22

592

$3.0\times 10^{-7}$

13.9

$\rm III$

(0.99,0.999]

15

$2.9 \times 10^{-8}$

0.21

28

$3.0\times 10^{-8}$

0.33

139

$3.2\times 10^{-8}$

3.03

264

$3.4 \times 10^{-8}$

4.69

554

$3.3\times 10^{-8}$

8.42

In particular spherical collapse provides quantitative predictions for "stellar'' properties for the entire accretion process starting with the formation of a photosphere-when a radiating outer layer forms around an optically thick central cloud region-to the quasi-hydrostatic pre-main sequence contraction. We have shown that properties of stars with "initial thermal structures'' which result from the collapse of cloud fragments differ from those obtained with simplifying assumptions like isentropic or polytropic radial structure. The latter are commonly used as a starting point for the calculation of early stellar evolution.

With our assumptions we can describe the pre-main sequence evolution as a consequence of the star formation process, more precisely as a result of the collapse of Bonnor-Ebert spheres. All "stellar'' masses appear at locations in the HRD that are different from the ones predicted by classical, fully hydrostatic calculations. In particular, they do not start at the top of a classical, initially fully convective track or at the birthline. Those differences cannot be attributed to differences in energy transfer treatment or differences in microphysics since we have constructed our model equations in a way that assures a "stellar structure'' limit, that is approached during pre-main sequence contraction and that coincides with the assumptions made in our classical fully hydrostatic reference-study.

The ages obtained with collapse models differ significantly, i.e., up to 1 Myr, from the ages derived from classical calculations. This is only in part due to the fact that the time needed for collapse and accretion are not accounted for in classical PMS calculations. The primary reason is that the contraction behaviour of the young stars resulting from the collapse differs from the one obtained for the usually assumed fully convective, isentropic, or polytropic initial structures. The fact that deuterium is almost completely burnt during the main accretion phase, i.e., before the PMS is reached for all "truly'' stellar masses calculated here (the $0.1~{M_\odot}$ case is the exception) also contributes to the difference. Hence calculations of early stellar evolution depend on the star formation process. Initial conditions are not forgotten in general. Another way to summarize this is that the star formation time, in our study about 3-4 free-fall times in all cases, is not sufficiently different from the Kelvin-Helmholtz time-scale of young objects to separate initial energy deposition into the star during the accretion phase from energy losses that control the contraction to the main sequence. As a consequence, protostellar collapse determines the PMS evolution.

This effect is enhanced by the off-centre ingnition of D burning and its interaction with the thermal structure. Unlike during most of stellar evolution it is not only controlled by the existing thermal structure but helps also to conserve the global thermal structure outside thermal equilibrium. This differs from the usual behaviour close to complete equilibrium, e.g., of main sequence stars, where nuclear reactions are adjusted to the thermal needs of the star and re-establish global thermal equilibrium.

Our study re-emphasises the fact that star formation is an intrinsically time dependent process, with variations in all flow quantities. We have tried to tentatively group the young population into four different evolutionary phases based on the fraction of the residual circumstellar mass around them. We gave characteristic quantities for those phases showing that mass accretion rates vary significantly over the main accretion phase. Accretion rates are typically a few $10^{-6}~ M_\odot/$yr in the earliest phase, a few $10^{-7}~ M_\odot/$yr for the last 10% and a few $10^{-8}~ M_\odot/$yr for the accretion of the last percent of the final mass. Those values are not too far from what is currently inferred from observations (e.g., Brown & Chandler 1999). Also our star formation times of 0.07 to $1.3~\rm Myr$ for 0.1 to $2~ M_\odot$ are consistent with estimates from a comparison of T Tauri star- and embedded source counts (Kenyon et al. 1990). Our median luminosities between about 0.3 to 10 $L_\odot$, depending on accretion phases, for cloud fragment masses from 0.01 to $0.5 ~{M_\odot}$, make it much harder to construct a "Class I luminosity problem'' than with the assumed constant, canonical mass accretion rate of $10^{-5}~M_\odot/\rm yr$, that is in conflict with observations (Brown & Chandler 1999). Let us finally summarize our main points: (1) the star formation processes, i.e., the protostellar collapse translates initial cloud conditions into initial stellar thermal structure, (2) pre-main sequence evolution depends qualitatively and quantitatively on the collapse flows; (3) D burning ignites and mostly occurs during the accretion phase for cloud fragments of 0.5, 1, 2, and $10~{M_\odot}$; (4) we confirm that protostars and PMS objects have an outer convective shell and a central radiative zone; (5) the radiative core lasts at least to a (new) age of 1.5 Myr for cloud fragments of 0.5, 1, and $2~ M_\odot$ and we expect off-centre ignition of hydrogen burning; (6) collapse pre-main sequence tracks are initially cooler than classical tracks. Later they cross the Hayashi line rapidly and start their quasi-classical decent along a mixed radiative/convective contraction track. After the final stellar photosphere becomes visible they appear to the left of the Hayashi lines. Typically at 0.05 dex hotter effective temperatures than the corresponding classical, fully hydrostatic pre-main sequence model; (7) already during the early pre-main sequence phase the solar mass case is roughly homologous to the Sun with its radiative core rather than to a fully convective star on the Hayashi track; (8) collapse ages below 2 Myr are up to a Myr older than classical ages for the same luminosity.

Acknowledgements

This work was supported by the Fonds zur Förderung der wissenschaftlichen Forschung (FWF) under project numbers S7305-AST and S7307-AST and the Deutsche Forschungsgemeinschaft (DFG), SFB 359 (collaborative research center of the German national science foundation on "Reactive flows, diffusion and transport''). GW thanks M. Feuchtinger for the collaboration on the convection model and sharing the version of the VIP-code and the computational environment that he developed, and H. Herndl for providing a subroutine for the calculation of D reaction rates. We thank E. Dorfi for developing and sharing the original version of the VIP-code and providing substantial amounts of computational resources on his VAX-cluster in Vienna where the computations were perfomed, and for comments and discussion. We thank Christian Straka for reading the manuscript. We are indebted to Prof. H. Petersmann$\dag$ from the Seminar für Klassische Philologie, University of Heidelberg, for a most elucidated advice on the semantics and etymology of the term isopleths.

References

 
• Alexander, D. R., & Ferguson, J. W. 1994, AJ, 437, 879 In the text  
• Allègre, C. J., Manhès, G., & Göppel, C. 1995, Geochim. Cosmochim. Acta, 59, 1445 In the text NASA ADS  
• Appenzeller, I., & Tscharnuter, W. 1975, A&A, 40, 397 In the text NASA ADS  
• Balluch, M. 1991a, A&A, 243, 168 In the text NASA ADS  
• Balluch, M. 1991b, A&A, 243, 187 NASA ADS  
• Balluch, M. 1991c, A&A, 243, 205 NASA ADS  
• Baschek, B., Scholz, M., & Wehrse R. 1991, A&A, 246, 374 In the text NASA ADS  
• Brown, D. W., & Chandler, C. J. 1999, MNRAS, 303, 855 In the text NASA ADS  
• Castor, J. I. 1972, ApJ, 178, 779 In the text NASA ADS  
• Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039 NASA ADS  
• Caughlan, G. R., & Fowler, W. A. 1988, Atom. Data Nucl. Data Tab., 40, 283 In the text NASA ADS  
• D'Antona, F., & Mazzitelli, I. 1994, ApJS, 90, 457 In the text  
• Dorfi, E. A. 1998, in Computational Methods for Astrophysical Fluid Flow, ed. O. Steiner, & A. Gautschy (Springer)  
• Dorfi, E. A., & Drury, L. O'C. 1987, J. Comp. Phys., 69, 175 NASA ADS  
• Dorfi, E. A., & Feuchtinger, M. U. 1995, Comp. Phys. Comm., 89, 69  
• Dorfi, E. A., & Feuchtinger, M. U. 1999, A&A, 348, 815 NASA ADS  
• Feuchtinger, M. U. 1998a, A&A, 337, L31  
• Feuchtinger, M. U. 1999a, A&AS, 136, 217 In the text NASA ADS  
• Feuchtinger, M. U. 1999b, A&A, 351, 103 NASA ADS  
• Feuchtinger, M. U., & Dorfi, E. A. 1994, A&A, 291, 225  
• Foster, N. P., & Chevalier, R. A. 1993, ApJ, 416, 303 NASA ADS  
• Grossmann, A. S., & Graboske, H. C. 1971, ApJ, 164, 475 NASA ADS  
• Gehmeyr, M., & Mihalas, D. 1993, American Astronomical Society Meeting, 183, #48.06 (abstract) http://141.142.3.131/Pubs/TechReports/TechReportsSeries.html  
• Götz, M. 1993, Dissertation, Univ. Heidelberg In the text  
• Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 896 In the text  
• Kippenhahn, R., & Weigert, A. 1990, Stellar Struct. Evol. (Springer, Berlin), 42 In the text  
• Kuhfuß, R. 1986, A&A, 160, 116 NASA ADS  
• Kuhfuß, R. 1987, Ein Modell für zeitabhängige Konvektion in Sternen, Ph.D.-Thesis, TU München In the text  
• Kürschner, R. 1994, A&A, 285, 897 NASA ADS  
• LeVeque, R. 1998, in Computational Methods for Astrophysical Fluid Flow, ed. O. Steiner, & A. Gautschy (Springer)  
• Mihalas, D., Däppen, W., & Hummer, D. G. 1988, ApJ, 331, 815 In the text NASA ADS  
• Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Oxford University Press, New York) In the text  
• Mihalas, D. 1998, in Computational Methods for Astrophysical Fluid Flow, ed. O. Steiner, & A. Gautschy (Springer)  
• Morfill, G., Tscharnuter, W. M., & Völk, H. 1985, in Protostars and Planets II, ed. D. C. Black, & M. S. Matthews (Tucson: Univ. of Arizona Press), 493  
• Tscharnuter, W. M. 1987, A&A, 188, 55 In the text NASA ADS  
• Tscharnuter, W. M. 1989, Low Mass Star Formation and Pre-main Sequence Objects, ed. Bo Reipurth (ESO, Garching bei München)  
• Tscharnuter, W. M. 1991, The Physics of Star Formation and Early Stellar Evolution, NATO ASI, Series C, 342, ed. C. D. Lada, & N. D. Kylafis (Kluwer Dordrecht), 411  
• Tscharnuter, W. M., & Winkler, K.-H. A. 1979, Comp. Phys. Comm., 80, 171  
• Weiss, A., Keady, J. J., & Magee, N. H., Jr. 1990, Atom. Data Nucl. Data Tab., 45, 209 In the text NASA ADS  
• Winkler, K.-H. A., & Newman, M. J. 1980a, ApJ, 236, 201 In the text NASA ADS  
• Winkler, K.-H. A., & Newman, M. J. 1980b, ApJ, 238, 311 NASA ADS  
• Winkler, K.-H. A., & Norman, M. L. 1986, NATO Advanced Research Workshop on Astrophysical Radiation Hydrodynamics, held in Garching bei München, Germany, August 2-13, 1982, ed. K.-H. A. Winkler, & M. L. Norman (Dordrecht: D. Reidel) In the text  
• Wuchterl, G. 1989, Die Entstehung von Gasplaneten, Ph.D. Thesis, Univ. Vienna In the text  
• Wuchterl, G. 1990, A&A, 238, 831992 In the text  
• Wuchterl, G. 1995a, Comp. Phys. Comm., 89, 119  
• Wuchterl, G. 1995b, Earth, Moon Planets, 67, 51 NASA ADS  
• Wuchterl, G., & Feuchtinger, M. U. 1998, A&A, 340, 419 In the text NASA ADS  
• Yorke, H. W. 1979, A&A, 80, 308 In the text NASA ADS  
• Yorke, H. W. 1980a, A&A, 85, 215 In the text NASA ADS  
• Yorke, H. W. 1980b, A&A, 86, 286 NASA ADS

 

Online Material

Appendix A: Equations for protostellar collapse and pre-main sequence evolution

Boundary conditions

At the outer boundary the moment equations of the RTE are closed by

$$J = B(T_{\rm ext}) + \frac{1}{\bar{\mu}} H$$ (A.1)

with $T_{\rm ext}=10~{\rm K}$ and $\bar{\mu}=1/2$ (Yorke 1980b),

       

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} \varrho ~ {\rm d}\tau... ...}S ) = 0 , \hspace*{7cm} \Delta M_{\rm r} = \int_{V(t)} \varrho ~ {\rm d}\tau ,$$ (A.2)

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} \varrho_{\rm D} ~ {\r... ... \frac{ A_{\rm D} }{ N_{\rm L} Q_{\rm D} } \varrho \epsilon_{\rm nuc}^{\rm D} ,$$ (A.3)

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} \varrho u ~ {\rm d}\t... ... ,\hspace*{3.4cm} C_{\rm M} = \int_V \kappa \varrho \frac{F}{c} ~ {\rm d}\tau ,$$ (A.4)

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} \varrho ( e + \omega ... ... = - C_{\rm E} + \int_{V(t)} \varrho \epsilon_{\rm nuc}^{\rm D} ~ {\rm d}\tau ,$$ (A.5)

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} E ~ {\rm d}\tau \righ... ...ace*{4.15cm} C_{\rm E} = \int_V \kappa \varrho ( 4 \pi S - c E ) {\rm d} \tau ,$$ (A.6)

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} \frac{F}{c^2} ~ {\rm ... ...ial r} \right) ~{\rm d}\tau = - C_{\rm M} , \hspace*{3.4cm} P = \frac{1}{3} E ,$$ (A.7)

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} \varrho \omega ~ {\rm... ...tial r} \Pi , \quad \tilde{S}_\omega = \frac{c_{\rm D}}{\Lambda} \omega^{3/2} ,$$ (A.8)

$$j_{\rm w} = \varrho T \Pi , \quad \Pi = \frac{w}{T} ~ u_{\rm c} F... ...= \frac{c_{\rm p}~\kappa~\rho^2~\Lambda^2}{4 \sigma ~ T^3 ~ \gamma_{\rm R}^2} ,$$ (A.9)

$$\epsilon_{\rm nuc}^{\rm D} = \frac{ Q_{\rm D} }{ \varrho } \tilde... ...\rm M} \Lambda \omega^{1/2} \varrho \frac{\partial c_{\rm D}}{\partial r} \cdot$$ (A.10)

  

Table A.1: The equations are applied to spherical volumes. $S_\omega$, contains the Schwarzschild-Ledoux criterion via $-\partial s / \partial r = c_{\rm p}/H_{\rm p} ( \nabla - \nabla _{\rm s} )$. Convectively unstable stratifications occur in the present model when pressure and temperature gradients have the same sign and produce a positive value of $S_\omega$ that then contributes a source of turbulent kinetic energy, $\omega = 3/2 u_{\rm c}^2$ to the balance equation of turbulent kinetic energy. Note that $S_\omega$ can also act as a sink, corresponding to damping of convection due to damping bouyancy effects - a property not usually assured in convection models. For the flux-limiting function $F_{\rm L}[x] = 2/\pi \arctan[\pi/2 ~x]$ has been used, that satisfies $F_{\rm L}[x] \approx x$ for $x \ll 1$ and $F_{\rm L}[x] \approx 1$ for $x \gg 1$. The enthalpy, $w = e + P/\varrho$, is explicitely needed only for the specification of the upper limit of $\Pi$ and hence $j_{\rm w}$. See Wuchterl & Feuchtinger (1998) for a complete formulation of the convection model and definitions but note that a minus sign is missing in Eq. (5), for $S_\omega$. $j_{\rm D}$ is the D flux due to convective mixing, for D-concentration, $c_{\rm D} = \varrho _{\rm D} / \varrho$, mixing length $\Lambda$, and convective velocity $u_{\rm c}$. The proton partial density, $\varrho _{\rm p}$, is approximated by the hydrogen partial density, $\rm X \varrho$. Note that $A_{\rm D}/N_{\rm L} \tilde{r}_{\rm{}^2H(p,\gamma){}^3He} = A_{\rm D} / (N_{\rm L} Q_{\rm D}) \varrho \epsilon_{\rm nuc}^{\rm D}$, $\varrho _{\rm p} N_{\rm L} / A_{\rm p} \varrho _{\rm D} N_{\rm L} / A_{\rm D} = n_{\rm p} n_{\rm D}$. Convection parameters: $\alpha _{\rm ML} = 1.5$, $\alpha _{\rm S} = \alpha _{\rm M}= 1/2 \sqrt {2/3}$, $c_{\rm D} = 8/3 \sqrt {2/3}$, $\beta _r = 0.1$, $\gamma _{\rm rad}=2 \sqrt {3}$.

----------

that is, radiation is streaming freely into an external reservoir emitting blackbody radiation of $10 ~ \rm K$. All other surface terms vanish both at the outer cloud boundary (e.g., no convective flux leaves the fragment) and at the centre, where the interior mass $M_{\rm r}$ is zero. The total volume of the computed domain is assumed to be constant, equal to the volume of the critical Bonnor-Ebert sphere.

Appendix B: Numerical techniques and improvements

The equations are solved with a "protostellar'' version of the VIP-code (Dorfi & Feuchtinger 1995), where convection is added as outlined by Wuchterl (1995a). The equations are discretized following the principles described by Winkler & Norman (1986) and Tscharnuter (1987, 1989) with a tensor pseudo viscosity (Tscharnuter & Winkler 1979) and a self-adaptive grid determined by an equation of the Dorfi & Drury (1987) class. For a recent comparative discussion of discretisation techniques see LeVeque (1998) and for an overview over the self-adaptive grid method see Tscharnuter (1991) and Dorfi (1998). The VIP code contains improvements with respect to other self-adaptive grid codes used to calculate the protostellar collapse. The most important one being a higher order (van Leer) advection scheme (see Dorfi & Feuchtinger 1995). The key advantage in the present context is that a higher order advection scheme allows the grid to move faster, for a given advection error. Hence grid points can be concentrated more rapidly and higher spatial resolution can be obtained faster. The advection errors caused by local grid refinement are an important "side effect'' of using an adaptive grid when very high resolution is obtained on timescales that are much smaller than the physical timescales in other parts of the flow (Kürschner 1994). This side effect is only present when the flow is non steady and, as a consequence, the motion of the gridpoints contributes to the advection errors.

For our calculations of the protostellar collapse, modifications of the VIP-code have been made. They can be understood as a consequence of the quantitative resolution requirements and accuracy demands of a protostellar collapse calculation. Many of these requirements have been discussed previously (Tscharnuter & Winkler 1979; Tscharnuter 1987; Morfill et al. 1985; Balluch 1991a-c; Kürschner 1994). We will briefly describe the modifications as motivated by a view based on the properties of self-gravitating equilibrium gas spheres (with radiation) as they are imposed onto those systems by their gravity. This emphasises the accuracy requirements that become increasingly important toward the end of the star formation process. Increasing fractions of the cloud fragment mass actually settle into hydrostatic equilibrium as the pre-main sequence is approached. Earlier high resolution studies focused on the accurate representation of the flow discontinuities (e.g., Balluch 1991a-c) because the hydrostatic parts were relatively easy to resolve with a self-adaptive grid. Nevertheless, they deserve a careful consideration because they contain a large fraction of the total mass. After all, that fraction should approach unity at the end of a star formation calculation.

B.1 Spatial resolution

A hydrostatic sphere of gas with given temperature, be it a cloud fragment, a star, a brown dwarf or a giant planet, has a typical "pressure'' scale height, $H_{\rm P} = - P / (g \varrho) = c_{\rm T}^2 / g \sim \frac{e}{g}$, where pressure, P, gravitational acceleration, g, and density, $\varrho$ usually vary both in space and time (see Baschek et al. 1991 for a general formulation). $c_{\rm T}$ and e are the isothermal sound speed and the internal energy per unit mass. This mechanical structure is imposed by gravity directly via the force balance of the hydrostatic equilibrium, independent of the details that lead to the specific pressure and density or sound speed in the gas sphere. If the pressure is high the influence of gravity is relatively weak. High pressure systems with small gravity almost behave non-gravitating, with no radial gradients in thermodynamic equilibrium. In our context of star formation the typical example for such a weakly gravitating equilibrium system is the marginally stable isothermal cloud fragment. Length scales below a Jeans wavelength, $\lambda_{\rm Jeans} = \pi c_{\rm T}^2 / (G\varrho)$, are stabilized by gas pressure. A marginally stable cloud with radius $\lambda_{\rm Jeans}$ has a minimum pressure scale height of $H_{\rm P}^{\rm Jeans} = c_{\rm T} ( 16/9 \pi^3 G \varrho )^{-1/2} \sim c_{\rm T} \tau_{\rm ff}$. If we compare the Jeans scale height to the Jeans length we obtain a ratio of the length scale over which the pressure changes relative to the overall cloud size, i.e., $H_{\rm P}^{\rm Jeans} / \lambda_{\rm Jeans} = 3/(4 \pi^2) \sim 0.1$. In other words, a Jeans critical cloud is quasi-homogeneous, with pressure changes within an order of magnitude.

At the end of a star formation calculation gravity becomes much more important for the overall structure, especially in the outermost layers of the stars, where temperature and pressures are smallest. After all, the dominating effect of gravity is the key reason of why spherical symmetry is a good approximation through most parts of stellar evolution. For the Sun we obtain a ratio of the photospheric pressure scale height to the solar radius, of $\left. H_{\rm P}/R \right\vert _\odot \sim 2\times 10^{-4}$.

For a star formation calculation this means that, as the cloud collapses, gravity is becoming more and more dominating and the scale heights are becoming correspondingly smaller. In general the respective dimensionless number to specify the relative length scale over which flow variables vary in a gravitating system close to hydrostatic equilibrium is given by the compactness,

$$\frac{r}{H_{\rm P}} = \frac{ r_{\rm acc} }{r} = \frac{ GM }... ...rac{ R_\odot }{r} \right) \left( \frac{M}{ M_\odot } \right)$$ (B.1)

the meaning becoming more obvious with the $r_{\rm acc}$ expression. $c^2_{\rm T,\odot}$ is the isothermal sound speed for solar effective temperature. In any numerical calculation the local pressure scale height has to be resolved by a few "zones''. This translates into a numerical resolution requirement estimating the length scales which are introduced by gravity alone, i.e., when all dynamical effects are neglected (e.g., no ram pressure, no shocks, etc.). Traditionally this has been met by using the pressure as a coordinate and discretizing in $\log P$, but this is only possible with the strictly monotonic pressure P(r) produced by hydrostatic equilibrium. That approach cannot be used in hydrodynamical studies because pressure cannot be used as a global coordinate.

When hydrostatic parts and a circumstellar gas flow are calculated together, as is the case in our study, the compactness directly translates into a lower limit for the required local radial resolution (Wuchterl 1995a):

$$\frac{r}{\Delta r} > \frac{r}{H_{\rm P}} \cdot$$ (B.2)

Any calculation of the forming Sun has to resolve the outer stellar layers with at least the respective local refinement of $\sim$10⁴. To make progress on an evolutionary timescale the timestep should be at least on the global dynamical timescale $R/c_{\rm s}$ with $\Delta r / r = 10^{-4}$ and hence the numerical scheme must be able to cope with Courant numbers of 10⁴ to at least study global dynamical phenomena. To follow the evolution on the Kelvin-Helmholtz time-scale,

$$\tau_{\rm KH,\odot} = \frac{G {M^2_\odot}}{R_\odot L_\odot} \sim 30 ~ {\rm Myr} ,$$ (B.3)

of thermal contraction-which is typically 10 orders of magnitude larger than the dynamical time-scale for the Sun-the numerical scheme has to progress at Courant numbers of the order of 10¹⁴. Hence any scheme that can provide the necessary resolution has to be implicit.

A hydrodynamic calculation that does not satisfy the resolution requirement (.2) cannot distinguish between a dynamic pressure gradient caused, e.g., by a sound wave and the hydrostatic "background'' structure! In practice there is, of course, much more flow structure to be resolved on a pressure scale height. Radiative transfer, e.g., requires a resolution in optical depth that sometimes is much more stringent as it requires $\Delta \tau < 0.01$ (cf. Balluch 1991a). The condition (.2) above and the resulting Courant numbers therefore must be viewed as a demanding but still relatively weak lower limit. In practice resolutions of $r/\Delta r \sim 10^{10}$ are needed towards the end of the calculations presented here.

The compactness is a dimensionless number that can be used to quantify both the effects of gravity and the respective resolution requirements in astrophysical radiation hydrodynamics. It can also be used to distinguish between systems with different importance of gravity effects (Wuchterl 1995a). The analogy to the square of the Mach number describing the importance of inertial effects is displayed below:

$$\mbox{Mach}^2 = \frac{u^2}{c^2} \sim \frac{ e_{\rm kin} }{ e_{\rm therm} } \frac{ r_{\rm acc} }{r} = \frac{GM}{ r c^2_{\rm T} } \sim \frac{ e_{\rm grav} }{ e_{\rm therm} } \cdot$$ (B.4)

The collapse flow leading to the formation of a $1~{M_\odot}$ star, for instance, involves Mach numbers up to 400 and a compactness of 10⁴. For comparison: Accretion flows on proto giant planets are much easier to calculate because Mach number and compactness involved are both around 10 (cf. Wuchterl 1995b). RR Lyrae and $\delta$ Cephei stars have compactness between a few 100 and 1000 and the Mach numbers for the pulsation-flows are around 10. White dwarfs have compactness ranging from $5\times 10^{4}$ to $2.5\times 10^{5}$; neutron stars finally approach maximum compactness, $c^2/2c_{\rm T}^2$ with radii of a few Schwarzschild radii, $R_{\rm S} = 2GM/c^2$.

B.2 Advection

For a calculation of protostellar collapse the condition (.2) means that, from the equilibrium structure of gravitating systems alone, we have to resolve at least structures starting from 0.1 to 10^-4 of the local radius. The latter is about 10^-10 of the extension of the computational domain, the initial volume of the cloud fragment. In practice this is needed to describe the "background'' structure on top of which the radiation hydrodynamics of the collapse proceeds. In principle, the required resolutions are no problem for a self-adaptive grid method and we could solve the Navier-Stokes equations with the "molecular'' viscosity and omitting any pseudo viscosity as demonstrated by Balluch (1991c). But this pushes present computer floating point operations to its limits (quadruple precision is then needed and slows down the computations considerably), and the "realistic'' molecular viscosity might cause large advection errors during time-dependent flow phases. To meet the resolution requirements of gravitating equilibrium systems while keeping the advection errors as small as possible we use a pseudo viscosity length scale, $l_{\rm visc}$, defined by

$$\frac{1}{ l_{\rm visc} } = \frac{1}{\alpha_{\rm visc}} \left( \frac{1}{r} + \frac{1}{H_{\rm P}} \right) \cdot$$ (B.5)

The harmonic mean with r assures that the length-scale for pseudo viscosity is specified in units of the local radius in weak gravity parts of the flow. Calculations presented here are made for $\alpha_{\rm visc} = 0.01$, comparison calculations to check the non-dependence of the results on this computational parameter have been made for values down to 10^-4. The pseudo viscosity is a numerical tool that is used to calculate viscous dissipation in shocks on length scales larger than those of the molecular mean free path, where it would occur naturally. The resolution requirements for the numerical solving method are relaxed by the corresponding increase in pseudo viscosity relative to molecular viscosity. Of course, it would be satisfactory to adopt the "realistic'' values for molecular viscosity and start out with the Navier-Stokes equations directly, but the resolution requirements are very high $r/\Delta r \sim 10^{11}$ (Tscharnuter 1991), even very early in the collapse after first (outer) core formation. The consequences of such high resolution are either very high grid point numbers or advection errors during time-dependent phases, when a self-adaptive grid is used. Both the high resolution and the high advection speeds cause an additional computational effort. To avoid this we use the above $l_{\rm visc}$ to obtain a compromise between high resolution, large grid point numbers, large advection speeds and large floating point accuracy. With our choice of computational parameters it is assured that shock dissipation takes place on length scales that are two orders of magnitude smaller than the natural length scales of gravitating equilibria. If shocks appear our flow model starts to become viscous on that small relative scale. This avoids the "mixing of scales'' and has the advantage that calculations can be carried out using only double precision (14-16 digit) floating point operations and makes it possible to calculate the collapse with only 1000 grid points. As the compactness of the forming star increases our maximum resolution gradually increases and advection errors are kept to the necessary minimum.

The previous considerations guarantee that the spatial resolution is sufficiently high but low enough to avoid unnecessary advection errors when high resolution is required fast. In other words, grid induced advection is avoided by considering the sources for grid motion during time-dependent evolutionary phases. We also improved the advection scheme by demanding volume conservation, i.e.,

$$\frac{{\rm d}}{{\rm d}t} \left[ \int_{V(t)} {\rm d}\tau \right] + \int_{\partial V} u_{\rm grid} \cdot {\rm d}S = 0$$ (B.6)

to be satisfied algebraically by the advection scheme. This demand results in an expression for $u_{\rm rel}$, the velocity of the flow relative to the grid

$$u_{\rm rel} = u - u_{\rm grid} = u - \frac{1}{3} \frac{\delt... ...} \frac{1}{\delta t} \left( r - \frac{(r^*)^3}{r^2} \right)$$ (B.7)

to be used in the discretisation of the advection terms (see Dorfi 1998 for a geometric discussion). We use advection fluxes calculated from the various densities and velocities in all equations rather than using the temporal change in interior mass (as advocated by Winkler & Norman 1986 and used in the standard version of the VIP-code, cf. Dorfi & Feuchtinger 1995). The calculation of the advection terms from the velocities results in a slightly modified discretisation of the equation of motion (no spatial averaging in the momentum densities). As a consequence of this modification of the advection scheme the advection errors reported by Kürschner (1994), which produce step-like features when grid speeds are very high become negligible even for grids that are similar to the grids he used.

B.3 Hydrostatic parts and total energy

To improve total energy conservation we use the freedom that is available in formulating the grid equations after Dorfi & Drury (1987). First we limit the grid speed by a time-constant for the grid

$$\tau_{\rm grid} = g_\tau \left[ t_{\rm ff}\left(\frac{M}{4\pi r^3}\right) + 1000 \right] [\rm s] ,$$ (B.8)

where $t_{\rm ff}(\varrho)$ is the free-fall time for the enclosed mass, i.e., we specify the time constant for the grid in units of the free-fall time corresponding to the local mean density of the gravitating mass. The weighting factor $g_\tau$ is set 100 to assure that the relative Mach numbers of the grid remain below $\sim$100 everywhere. This reduces the total energy errors due to advection that have been identified (Kürschner 1994) to significantly contribute to the oscillation amplitudes found for first protostellar cores.

Secondly, to enhance the resolution in hydrostatic parts of the flow that carry large fractions of the total energy (and total energy errors, cf. Kürschner 1994) we weight the arc-length in our desired resolution function (sometimes also called the monitor function) for the grid equation motivated by a scaled estimate for the total potential energy content of the respective mass shell,

$$g_{\rm l} = \tanh \left( \frac{ r^2 M_{\rm r} }{ {R^*}^2 M^* } - 1 \right) + 1 ,$$ (B.9)

$$R_{\rm grid} = \Phi_{\rm grid} = \sqrt{ 1 + \sum_f g_{\rm l} g_{\rm f} \left( \frac{Df}{Dr} \right)^2 } ,$$ (B.10)

where Df,Dr: $Dx = \Delta x / \bar{x}$ are the usual, harmonically scaled differences of the flow variables ( $f \in \{ \kappa, \rho, s, T, \nabla_{\rm s} \}$) in the VIP code that enter the desired resolution. Gradients in places with large $E_{\rm grav}$ are strongly weighted in the arc-length for the desired resolution function. For the weights in the resolution function we use $g_{\rm f}=0.3$ for the results shown, and used $g_{\rm f}$-values down to 0.01 in comparison calculations (cf. RC1M vs. RC1ML14 in Table 1).

Our general method to distinguish between different candidate prescriptions for the "arbitrary'' or "personally biased'' part of the Dorfi & Drury (1987) grid equation was to use the total energy conservation properties as a criterion. A better grid is considered to be one that conserves total energy better. By this approach the freedom in specifying the desired resolution function has been somewhat reduced compared to other studies.

B.4 Discretization remarks

The discretisation of the nuclear energy generation source term to the internal energy equation is done analogously to other source terms. The D balance equation is discretized in analogy to the continuity equation. We note that we use D concentration as our additional computational variable and that we use a discretization for $j_{\rm D}$ which is symmetric in the grid index.

B.5 Constructing the initial conditions

The critical, i.e., marginally gravitationally unstable Bonnor-Ebert sphere was determined numerically in the following way: We started with a homogeneous, isothermal sphere of one tenth of the Jeans-critical density. After dynamical relaxation this sphere becomes hydrostatic and quasi-homogeneous with only a weak central density enhancement. The calculations are carried out for a given volume (originally, the volume corresponding to 0.1 of the Jeans density) that is slowly contracted (at $\dot{R} = 0.01~ {\rm cm~s^{-1}}$). The subcritical hydrostatic sphere therefore shrinks on the time-scale $R/\dot{R}$, which is $\sim$ $10^{19} ~{\rm s} \sim 320 ~ {\rm Gyr}$. The calculations (with the fully dynamic equations) proceed on this time-scale through a sequence of essentially static isothermal spheres, $L \sim 10^{-10}~{L_\odot}$ until the marginally stable equilibrium sphere is reached and collapse sets in. This procedure guarantees that the initial cloud condition is accurately static and that the critical state for the detailed constitutive relations (most prominently the pressure equation of state) is determined precisely.