A&A 395, 321-338 (2002)
DOI: 10.1051/0004-6361:20021071

Gravitational collapse of nonsingular logatropic spheres

L. Di G. Sigalotti¹ - F. de Felice² - E. Sira¹

1 - Instituto Venezolano de Investigaciones Científicas, IVIC, Apartado 21827, Caracas 1020A, Venezuela
2 - Dipartimento di Fisica Galileo Galilei, Università di Padova, via Marzolo n. 8, 35131 Padova, Italy

Received 2 May 2002 / Accepted 19 July 2002

Abstract
We present the results of high-resolved, hydrodynamic calculations of the spherical gravitational collapse and subsequent accretion of nonsingular subcritical and critical A=0.2 logatropes, starting with initial configurations close to hydrostatic equilibrium. Two sequences of models with varying masses and the same central temperature $T_{\rm c}=10$ K are defined, which differ only in the fiducial value of the truncation pressure ( $p_{\rm s}/k=1.3\times 10^{5}$ cm^-3 K and $1.0\times 10^{7}$ cm^-3 K). In all cases, we follow the calculations until the central protostar has accreted 99% of the total available mass. Thus, the models may be indicative of early evolution from the Class 0 to the Class I protostellar phase. We find that the approach to the singular density profile is never entirely subsonic. In the lower $p_{\rm s}$ sequence, about 6% of the mass collapses supersonically in a $1~M_{\odot}$ sphere, while only $\sim$0.02% behaves this way in a critical ($\approx$92.05 $M_{\odot}$) logatrope. In the high $p_{\rm s}$sequence the same trend is observed, with $\sim$0.7% of the mass now infalling supersonically at the time of singularity formation in a $1~M_{\odot}$ sphere. Immediately after singularity formation, the accretion rate rises steeply in all cases, reaching a maximum value when the central protostar has accreted $\sim$40% of its final mass. Thereafter, it decreases monotonically for the remainder of the evolution. Our models predict peak values of ${\dot M}_{\rm acc}$ as high as $\sim$ $5 {-} 6\times 10^{-5}~M_{\odot}$ yr^-1 for logatropes close to the critical mass. In contrast, a subcritical $1~M_{\odot}$ logatrope reaches a maximum value of $\sim$ $8\times 10^{-7}~M_{\odot}$ yr^-1 for the lower $p_{\rm s}$ sequence compared to $\sim$ $5\times 10^{-6}~M_{\odot}$ yr^-1 for the higher $p_{\rm s}$ case. The results also imply that the accretion lifetimes are longer in logatropes with lower $p_{\rm s}$, consistent with the observational evidence that star formation in clumped regions occurs on shorter timescales compared to more isolated environments.

Key words: hydrodynamics - methods: numerical - stars: formation - circumstellar matter

1 Introduction

Despite recent progress in the field of star formation, the exact way by which stars form out of a molecular cloud is still unclear. The early stages of star formation, which observationally correspond to the gap separating starless cloud cores from young embedded protostars, bracket in time the onset of protostellar collapse within star-forming regions. A correct theory of star formation must therefore give a complete description of the gravitational collapse of molecular cloud cores. We adopt herein the terminology employed by Reid et al. (2002, hereafter RPW), in which a molecular cloud is referred to as an extended region containing gas condensations called clumps. In turn, the clumps may contain dense cores, or cloud cores, from which single and binary stars may form. The central, densest part of a core shall be denoted as a protostellar object or protostar.

The isothermal collapse of uniform-density spheres was first described by Bodenheimer & Sweigart (1968), Larson (1969) and Penston (1969), who performed hydrodynamical collapse calculations over density ranges relevant to the initial stages of star formation. Larson and Penston also derived semianalytically similarity solutions for the isothermal collapse in which both the density and velocity profiles approach invariant forms. At large radii, or at times close enough to the upper limit t=0, the Larson-Penston (LP) solution approaches a singular r^-2 density profile with a constant inflow velocity equal to $-3.3c_{\rm s}$, where $c_{\rm s}$is the isothermal sound speed. The limit t=0, also referred to as the core formation time, is exactly the time at which the density profile becomes singular everywhere. Hunter (1977) extended the LP solution through t=0 and found that for t>0, the flow approaches a free-fall collapse with the velocity varying as r^-1/2 and the density approaching an r^-3/2 profile. In a further paper, Shu (1977) derived a separate set of self-similar solutions for the isothermal collapse of a sphere. In particular, one of his solutions - the so-called expansion wave - has played a central role in the development of the standard theory of star formation (Shu et al. 1987). This solution begins at t=0 with a singular isothermal sphere (SIS) of density $\rho (r)=2c_{\rm s}^{2}(4\pi G)^{-1}r^{-2}$ in hydrostatic equilibrium, where G is the gravitational constant. The sphere first collapses at the center and then the infall gradually propagates outwards. This is precisely the inside-out collapse or expansion-wave solution. However, the SIS model is not applicable to more realistic clouds with finite central densities such as the marginally stable Bonnor-Ebert sphere (Ebert 1955; Bonnor 1956). It also predicts a time-independent central protostellar accretion rate which seems to be inconsistent with the observed luminosities of young stars (Kenyon et al. 1994). Such a constant rate is, however, too low to form massive stars on the timescales suggested by observations (Caselli & Myers 1995; McLaughlin & Pudritz 1997, hereafter MP97).

Hydrodynamic calculations of the collapse of critically stable Bonnor-Ebert spheres were made by Foster & Chevalier (1993), who showed that only the innermost 0.05% of the mass follows the LP solution just prior to the singularity formation (t=0). They also found that about 44% of the mass falls supersonically with a velocity of $-3.3c_{\rm s}$ at the time of singularity formation. One problem with this picture is that such high infall velocities has never been detected. Magnetic fields offer a way to slow down the collapse prior to singularity formation, but rapid infall may still occur thereafter. However, observational evidence in support of the isothermal models has recently been given by Bacmann et al. (2000) and Alves et al. (2001), who found, respectively, that the starless core L1544 in Taurus and the dark cloud Barnard 68 can be very well fitted by the structure of pressure-bounded Bonnor-Ebert spheres. On the other hand, there is strong evidence that the nonthermal component of the total velocity dispersion dominates over the thermal one in massive cores (Caselli & Myers 1995). Significant nonthermal motion has also been detected in low-mass cores (Fuller & Myers 1992). The properties of such cores cannot be explained by a simple isothermal equation of state (EOS).

As an alternative to the isothermal EOS, McLaughlin & Pudritz (1996, hereafter MP96) introduced a pure logatropic EOS in which the pressure depends on density logarithmically. MP96 showed that the logatropic EOS provides the best available phenomenological description of the internal structure and average properties of both molecular clouds and cores of low and high mass. In particular, a logatrope is consistent with the observed internal velocity-dispersion profiles of cloud cores and clumps from a variety of environments. This relation is meant to account for all contributions of the total gas pressure, including the effects of disordered magnetic fields (MHD turbulence). MP96 and MP97 also derived solutions for the equilibrium and self-similar collapse of logatropic spheres. In the limit of large radii, or small positive times, they found an entire family of singular solutions. At $t=-\infty$, the cloud is nonsingular with $\rho\propto r^{-1}$ in the outer layers. As it collapses, it approaches a singular profile $\rho\propto r^{-1}$ everywhere by t=0. At this time, the solution coincides with the analytical singular solution to the equation of hydrostatic equilibrium for a self-gravitating logatrope (MP96). In the singular logatropic sphere (SLS), the inflow velocity vanishes in accordance with the physical picture at t<0 of quasistatic evolution towards the singular profile. This situation is in clear contrast with the Bonnor-Ebert sphere, where the flow approaches the $\rho\propto r^{-2}$ profile with supersonic velocities. In the opposite limit of small radii, or large positive times (t>0), the density and velocity near the center tend to $\rho \propto r^{-3/2}$ and $v\propto r^{-1/2}$, respectively, similar to the LP solution. The unstable singular equilibrium with $\rho\propto r^{-1}$everywhere is the analogous of Shu's (1977) unstable SIS solution, while the stable nonsingular equilibrium solution is the counterpart of the hydrostatic Bonnor-Ebert sphere.

Recently, RPW performed three-dimensional hydrodynamic calculations of the collapse of both singular and nonsingular logatropic spheres with a choice of the kinetic temperature and surface pressure appropriate for isolated star formation (MP97). They were able to follow the collapse and subsequent accretion phase of only low-mass logatropes because of limited spatial resolution in their calculations. In this paper, we extend the work of RPW by conducting one-dimensional (1D), hydrodynamic simulations with high spatial resolution to follow the collapse and subsequent accretion of massive, nonsingular logatropes using the same fiducial central temperature and truncation pressure employed by MP97 and RPW. A set of calculation models with a much higher surface pressure representative of clustered star formation is also presented. In Sect. 2, we outline some relevant aspects of the logatropic collapse theory introduced by MP96 and MP97. The initial conditions and methods employed in our model calculations are described in Sects.  3 and 4. Section 4 also contains the results of some comparison test cases, including the collapse of a singular $1~M_{\odot}$logatrope. The results of our collapse calculations for both subcritical and critical nonsingular logatropes are described in Sects. 5 and 6. Finally, in Sect. 7 we compare our predictions with observations and discuss our results, while Sect. 8 contains the conclusions.

2 Theory of logatropic collapse

In this section we briefly review the theory of MP96 and MP97 to facilitate discussion throughout the paper.

2.1 The logatropic gas sphere

MP96 introduced the "pure'' logatropic EOS

$$p=p_{\rm c}\left[1+A\ln\left(\frac{\rho}{\rho _{\rm c}}\right)\right] ,$$ (1)

where $p_{\rm c}$ and $\rho _{\rm c}$ are, respectively, the central cloud pressure and density, and A is a positive, free parameter. This relation was formulated to provide a unified treatment of both molecular clouds and cores of low and high mass. In particular, it guarantees that the total velocity dispersion $\sigma ^{2}=p/\rho$ increases with decreasing density and increasing radius, in accordance with observations that indicate that the line widths are essentially thermal near the center of real cloud cores and nonthermal on larger scales. In this way, the central velocity dispersion in a logatrope is identified with the thermal value

$$\sigma ^{2}_{\rm c}=\frac{p_{\rm c}}{\rho _{\rm c}}=\frac{kT_{\rm c}}{\mu m_{\rm H}} ,$$ (2)

where k is the Boltzmann constant, $T_{\rm c}$ is the central temperature of the cloud core and $\mu$ is the mean molecular weight. Note that when $\rho =\rho _{\rm c}$, Eq. (1) reduces to $p_{\rm c}=\sigma ^{2}_{\rm c}\rho _{\rm c}$. The sound speed in a logatrope is given by

$$c^{2}_{\rm s}=\frac{{\rm d}p}{{\rm d}\rho}=\frac{Ap_{\rm c}}{\rho}$$ (3)

except at the center, where $c^{2}_{\rm s}=\sigma ^{2}_{\rm c}$. Thus, it is precisely at the center that a correspondence between the isothermal and logatropic EOS is established. The best fit of Eq. (1) to the observed internal velocity-dispersion profiles for low mass (Fuller & Myers 1992) and high mass (Caselli & Myers 1995) cores was obtained for $A=0.20\pm 0.02$ (MP96). For the calculations of this paper we use A=0.2.

Once the value of A is specified, we can determine the equilibrium structure of a logatropic sphere by using Eq. (1) along with the differential equations governing the motion of a self-gravitating fluid. Under the assumption of spherical symmetry, these are: the equation of continuity

$$\frac{\partial\rho}{\partial t}+\frac{1}{r^{2}}\frac{\partial (r^{2}\rho v)} {\partial r}=0,$$ (4)

Euler's equation of motion

$$\frac{(\partial\rho v)}{\partial t}+ \frac{1}{r^{2}}\frac{\pa... ...{\partial p}{\partial r}-\rho\frac{\partial\Phi}{\partial r} ,$$ (5)

and Poisson's equation

$$\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2} \frac{\partial\Phi}{\partial r}\right)=4\pi G\rho,$$ (6)

where $\rho =\rho (r,t)$ is the mass-density, v=v(r,t) is the fluid velocity, p=p(r,t) is the gas pressure and $\Phi =\Phi (r,t)$ is the gravitational potential.

In a state of exact hydrostatic equilibrium, the velocity v and the temporal derivatives in Eqs. (4) and (5) vanish. Using Eqs. (1) and (6), the equation for the balance of forces can be solved analytically to yield the singular density profile

$$\rho (r)=\rho _{\rm c}\sqrt{\frac{2A}{9}}\left(\frac{r}{r_{0}}\right)^{-1}= \sqrt{\frac{Ap_{\rm c}}{2\pi G}}r^{-1},$$ (7)

where $r_{0}=3\sigma _{\rm c}/(4\pi G\rho _{\rm c})^{1/2}$ is the characteristic or scale radius of the spherical core. This solution, which is the analogous of the SIS solution, defines an unstable equilibrium for the SLS. A second solution can also be obtained which represents a pressure-bounded sphere of finite central density. This solution can only be found through numerical integration of the hydrostatic equation (see Sect. 3) and is the logatropic analogous of the Bonnor-Ebert sphere.

MP97 also derived self-similar solutions for the collapse of logatropic spheres in terms of the similarity variable x=r/(a) and the dimensionless density, velocity and mass variables

$$\alpha (x) 4\pi Gt^{2}\rho (r,t),$$ =

$$\frac{v(r,t)}{a_{t}},$$ u(x) = (8)

$$\frac{GM(r,t)}{a^{3}_{t}t},$$ m(x) =

respectively, where M(r,t) is the total mass within a radius r at time t and a_t is a normalization velocity given by $a_{t}=[Ap_{\rm c}(4\pi Gt^{2})]^{1/2}$. Here, $p_{\rm c}$ refers to the central pressure in the initial, precollapse core ($t=-\infty$). For large radii, or small positive times, MP97 obtained analytically a family of singular solutions having the form

$$\rho (r,t)\to C\sqrt{\frac{Ap_{\rm c}}{4\pi G}}r^{-1}\hspace{... ...\frac{C}{2}\left(1-\frac{2}{C^{2}}\right)a_{t}\hspace{0.1 cm}.$$ (9)

The core starts to collapse at $t=-\infty$ with a density profile at large radii given by Eq. (9). For inflow v<0, and hence $C>\sqrt{2}$. When $C=\sqrt{2}$, the similarity solution reduces to the analytic singular profile (7) with v=0. This means that the core is essentially at rest at the time of singularity formation (t=0). On the other hand, at small radii ($r\to 0$), or large times ( $t\to\infty$), the solution becomes

$$\rho (r)\to\frac{{\dot M}}{4\pi (2GM_{0})}r^{-3/2}; \hspace{0.2 cm} v(r)\to (2GM_{0})^{1/2}r^{-1/2},$$ (10)

where M₀=M(0,t) is the mass of the central accreting protostar and ${\dot M}$ is the accretion rate. These solutions are only valid near the center and after the singularity formation time. The r^-3/2 and r^-1/2 scalings for $\rho$ and v, respectively, are similar to the LP solution and Shu's (1977) SIS collapse. The more significative difference between the SIS and SLS solutions lies in the form of the central accretion rate. From Eqs. (8), it follows that $M(0,t)\propto t^{4}$, yielding a central accretion rate ${\dot M}(0,t)\propto t^{3}$. In particular, the time-varying central protostellar accretion rate takes the form

$${\dot M}_{\rm acc}={\dot M}(0,t)=4(Ap_{\rm c}4\pi G)^{3/2}t^{3} \frac{m_{0}}{G},$$ (11)

where m₀ is the reduced (dimensionless) mass of the collapsing core. For the singular logatrope ( $C\to\sqrt{2}$), $m_{0}=6.67\times 10^{-4}$(MP97). The primary function of m₀ for a given SLS collapse is just to set the central mass accretion rate.

The collapse of the SLS is qualitatively similar to that of the SIS model. In both cases, the solution is characterized by an expansion wave propagating from the origin (r=0). That is, the mass shells of a collapsing SIS or SLS do not all fall in together, as in homologous collapse, but in sequence, from the inside out. Each collapsing shell removes the pressure support from the next, so the instantaneous locus where collapse starts is an expansion wave moving outwards at the sound speed. The fractional mass of the central accreting protostar varies as

$$\frac{M_{0}}{M_{\rm tot}}=\frac{\pi ^{4}}{16\sqrt{2}}\left( \frac{t}{t_{\rm ff}}\right)^{4}m_{0},$$ (12)

where $M_{\rm tot}$ is the total mass of the collapsing core and $t_{\rm ff}=(3\pi /32G\rho _{\rm ave})^{1/2}$ is the mean free-fall time of the initial ($t=-\infty$) core defined in terms of the average density. The expansion wave reaches the outer surface of the collapsing core at

$$t_{\rm ew}=\frac{4\sqrt{2}}{\pi}t_{\rm ff}\approx 1.8t_{\rm ff}$$ (13)

after which the total mass of the core is collapsing inwards. At this point, the similarity solution of MP97 predicts that the central protostar contains only 3% of the total available mass. Thereafter, the expansion wave is reflected and driven back into the core as a compression wave. Hence further accretion will no longer obey Eqs. (11) and (12) as was demonstrated by the numerical calculations of RPW. We defer discussion of this point to Sect. 4.5, where the collapse of a singular $1~M_{\odot}$logatrope is presented as a comparison test.

2.2 Recent observational evidence

There is some direct evidence that the internal structure of turbulent cores matches the predictions of logatropic equilibrium spheres. This justifies in part the use of numerical hydrodynamic calculations to examine the collapse of nonsingular logatropes of both low and high mass. For example, André et al. (1996) and Ward-Thompson et al. (1999) reported radial density gradients approaching $\rho\propto r^{-2}$ for $r\sim 5000 {-} 15~000$ AU, and flattening to $\rho\propto r^{-1}$ for $r\leq 2500 {-} 5000$ AU in individual cores. Convincing evidence has been given by Colomé et al. (1996) and Henning et al. (1998), who detected radial density gradients consistent with $\rho\propto r^{-n}$ for $n\approx 0.75 {-} 1.5$ in the envelopes of massive Herbig Ae/Be stars. Envelopes matching the density profile of collapsing logatropic spheres have also been reported by Osorio et al. (1999), who studied the dust thermal spectra of a number of hot cloud cores, which are thought to be the sites of massive star formation. More recently, van der Tak et al. (2000) studied the structure of the envelopes around 14 massive young stars, finding density profiles consistent with $n\approx 1.0 {-} 1.5$. Imaging of submillimeter dust emission from the envelopes of individual young protostars in Taurus and Perseus by Chandler & Richer (2000) also revealed density variations consistent with $n\approx 1.0 {-} 1.5$ for the youngest sources (Class 0 protostars) in their sample.

The above observational results point towards a prevalence of nonthermal support in massive young stellar objects, yielding power-law density variations that match the internal structure of both equilibrium (n=1.0) and collapsing (n=1.5) logatropic spheres. Thus, in contrast with most previous studies, the present logatropic collapse model calculations will apply to massive star formation.

3 Initial equilibrium models and physical parameters

The properties of A=0.2 pressure-bounded logatropes have been calculated by MP97, assuming fiducial values of the central temperature ( $T_{\rm c}=10$ K) and surface pressure ( $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K) appropriate for isolated star formation. For convenience in establishing the initial conditions and parameters used in our collapse calculations, we have recalculated the internal structure of A=0.2 pressure-truncated logatropes.

3.1 Stable nonsingular equilibria

The initial core model is chosen to be a logatrope close to hydrostatic equilibrium. For pressure-bounded spheres, the core properties are set by the central temperature $T_{\rm c}$, surface pressure $p_{\rm s}$ and the truncation radius R at which the internal pressure equals a confining pressure. For all models, we use A=0.2, $T_{\rm c}=10$ K and $\mu =2.33$. Two sequences of nonsingular collapse models with varying core mass and truncation radius are defined which differ only in the fiducial value of the surface pressure $p_{\rm s}$. The first sequence of models use $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K as in MP97 and RPW, while the second sequence employs $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K, which is more representative of the conditions observed in regions of clustered star formation.

For any of these sequences, we solve numerically the equation of hydrostatic equilibrium, as obtained from Eq. (5) by setting terms in v and $\partial /\partial t$ to zero, coupled with the logatropic EOS (1) and Poisson's Eq. (6). Defining the dimensionless radius $\xi$ and gravitational potential $\psi$ as (see Appendix A of MP97)

$$\xi =\frac{r}{r_{0}};\hspace{0.2 cm} \psi =\frac{\Phi -\Phi _{\rm c}}{\sigma ^{2}_{\rm c}},$$ (14)

where $\Phi _{\rm c}$ denotes the gravitational potential at the core center, the equation of hydrostatic equilibrium integrates to

$$\frac{\rho}{\rho _{\rm c}}=\frac{1}{1+\psi /A},$$ (15)

and the Poisson equation transforms into the Lane-Emden type equation

$$\frac{1}{\xi ^{2}}\frac{\rm d}{\rm d\xi}\left(\xi ^{2}\frac{\rm d\psi}{\rm d\xi}\right)= \frac{9}{1+\psi /A},$$ (16)

with boundary conditions such that at $\xi =0$, $\psi =\rm d\psi /d\xi =0$ and for $\xi >0$, $\psi >0$. For given A, the solution of Eq. (16) determines the density profile of equilibrium nonsingular logatropes. The radial integration of this equation was started from $\xi$ sufficiently close to the origin using an appropriate series expansion for $\psi$ valid when $\xi\to 0$. In all cases, we defined a sphere of maximum truncation radius $\xi =\xi _{\rm s}$ and subdivided its volume into N concentric shells of uniform width $\Delta\xi$. The value of $\Delta\xi$ was chosen such that the starting radius was always at a distance <<0.01 from the origin in order to get expansion values for $\psi$ and $\rm d\psi /d\xi$ satisfying Eq. (16) within less than 10^-8%. In practice, a value of $\Delta\xi =10^{-8}$ has been used yielding a solution for $\psi$ accurate within one part in 8.

Table 1 lists the properties of the first calculated sequence of equilibria (see Table 1 of MP97) using fiducial values of $p_{\rm s}$representative of isolated star-forming regions. The properties of equilibria using $p_{\rm s}$ consistent with cluster-forming regions are given in Table 2. The entries in each table specify, starting from the first column, the dimensionless truncation radius R/r₀, the central density and pressure contrasts, the ratios of mean density and mass averaged line width to the central values, the mass of the core in solar mass units and the mean free-fall time. The last row in each table lists the properties of the critically stable core with truncation radius given by (MP96)

$$\xi _{\rm crit}=\sqrt{\frac{2A}{9}}\exp\left(\frac{1}{A}- \frac{1}{4}\right)\approx 24.367\hspace{0.1 cm}.$$ (17)

The critical radius is not affected by changing the surface pressure because for any given value of A, a pressure-bounded logatrope admits only one critical mode. MP96 also gave an expression for the critical mass of an A=0.2 logatrope, namely

$$M_{\rm crit}=\frac{250}{\alpha ^{3/2}}~M_{\odot}\left(\frac{T... ...s}}{1.3\times 10^{5}~k~ {\rm cm}^{-3} ~{\rm K}}\right)^{-1/2},$$ (18)

where $\alpha$ is a parameter that depends on the strength of any magnetic field threading the core. In the absence of any magnetic field, MP97 gives the value $\alpha =35/18$. Using this value in Eq. (18) with $T_{\rm c}=10$ K and $p_{\rm s}=1.3\times 10^{5}k$ cm^-3  K yields a critical mass $M_{\rm crit}=92.203~M_{\odot}$ slightly larger than the value of $92.05~M_{\odot}$ in Table 1. This value is consistent with having $\alpha\approx 35.04/18$ in Eq. (18). Since $M_{\rm crit}\propto p_{\rm s}^{-1/2}$, an increase of $p_{\rm s}$ by two orders of magnitude lowers the critical mass by a factor of $\sim$10 (see Table 2). Furthermore, in accordance with the stability criteria derived by MP96, nonsingular logatropes of total mass $M_{\rm tot}<M_{\rm crit}$, or equivalently of radius $R<R_{\rm crit}$, will be stable against collapse. Cores with $M_{\rm tot}=M_{\rm crit}$, or $R=R_{\rm crit}$, will be marginally stable, while for $M_{\rm tot}>M_{\rm crit}$, or $R>R_{\rm crit}$, the cores will always be unstable to collapse.

In Fig.1 we show the singular and nonsingular density profiles for the critically stable logatrope in the sequence of Table 1. We see that the nonsingular profile follows the singular behavior closely over a broad range of outer radii.

 

 

Table 1: Properties of pressure-bounded nonsingular A=0.2 logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K.

$\frac{R}{r_{0}}$	$\frac{\rho _{\rm c}}{\rho _{\rm s}}$	$\frac{p_{\rm c}}{p_{\rm s}}$	$\frac{\rho _{\rm ave}}{\rho _{\rm c}}$	$\frac{\sigma ^{2}_{\rm ave}}{\sigma ^{2}_{\rm c}}$	$\frac{M_{\rm tot}}{M_{\odot}}$	$\frac{t_{\rm ff}}{10^{5}~{\rm yr}}$
0.931	4.279	1.410	0.335	2.330	0.5	4.305
1.340	6.219	1.576	0.238	2.997	1.0	4.835
1.954	9.161	1.795	0.164	3.896	2.0	5.458
2.874	13.573	2.090	0.111	5.046	4.0	6.142
3.615	17.117	2.315	0.088	5.836	6.0	6.554
4.262	20.204	2.507	0.075	6.446	8.0	6.844
4.848	22.999	2.682	0.065	6.945	10.0	7.063
7.306	34.690	3.440	0.043	8.589	20.0	7.666
9.386	44.565	4.156	0.034	9.553	30.0	7.908
11.305	53.669	4.916	0.028	10.177	40.0	7.982
13.163	62.482	5.780	0.024	10.584	50.0	7.944
15.034	71.352	6.827	0.021	10.824	60.0	7.812
17.002	80.687	8.204	0.019	10.918	70.0	7.578
19.231	91.255	10.281	0.016	10.851	80.0	7.200
20.288	96.268	11.551	0.016	10.761	84.0	6.977
21.578	102.384	13.467	0.015	10.603	88.0	6.664
22.432	106.436	15.040	0.014	10.472	90.0	6.429
24.367	115.611	20.019	0.013	10.100	92.05	5.808

	Figure 1: Density profile of the A=0.2 critically stable logatrope of Table 1. Both the singular and nonsingular profiles are shown.
Open with DEXTER

 

 

Table 2: Properties of pressure-bounded nonsingular A=0.2 logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K.

$\frac{R}{r_{0}}$	$\frac{\rho _{\rm c}}{\rho _{\rm s}}$	$\frac{p_{\rm c}}{p_{\rm s}}$	$\frac{\rho _{\rm ave}}{\rho _{\rm c}}$	$\frac{\sigma ^{2}_{\rm ave}}{\sigma ^{2}_{\rm c}}$	$\frac{M_{\rm tot}}{M_{\odot}}$	$\frac{t_{\rm ff}}{10^{4}~{\rm yr}}$
3.026	14.304	2.137	0.105	5.218	0.5	7.110
4.493	21.310	2.576	0.071	6.649	1.0	7.907
6.751	32.052	3.262	0.047	8.272	2.0	8.625
12.024	57.079	5.233	0.026	10.356	5.0	9.097
17.036	80.847	8.231	0.019	10.920	8.0	8.635
24.367	115.611	20.019	0.013	10.100	10.5	6.622

3.2 Unstable singular equilibrium

As a suitable comparison test case, we present high-resolution calculations of the collapse of a singular $1~M_{\odot}$ logatrope. The same model was previously calculated by RPW using 3D hydrodynamics. This test also allows for a direct comparison with the self-similar solutions of MP97.

The singular equilibrium density profile is given by Eq. (7). This form corresponds to $C=\sqrt{2}$ and $m_{0}=6.67\times 10^{-4}$, in the terminology of MP97. Since at r=0, the density becomes infinite care must be taken in handling the central region containing the singularity. One trivial way of skipping over this problem is to assign a finite density to this region by integrating Eq. (7) over a small radius and then dividing the result by the volume of the region. A similar approach was used by RPW to truncate the singular profile in their Cartesian-coordinate based calculations. This method yields for the central finite density the expression

$$\rho _{\rm c}=\frac{3}{2}\sqrt{\frac{Ap_{\rm c}}{2\pi G}} r^{-1}_{\rm c},$$ (19)

where $r_{\rm c}$ is chosen such that $r_{\rm c}<<R$. The smaller the value of $r_{\rm c}$, the more singular will appear the profile in the proximity of the center. If the surface density $p_{\rm s}$ and radius R of the singular core are set by the corresponding values in a nonsingular sphere, use of Eq. (7) gives a relation for the central pressure of the form

$$p_{\rm c}=\frac{2\pi G}{A}\rho ^{2}_{\rm s}R^{2}.$$ (20)

This can be inserted in Eqs. (7) and (19) to completely determine the singular density distribution. Thus, using the parameters in Table 1 for a $1~M_{\odot}$ nonsingular sphere, the best fit to a corresponding $1~M_{\odot}$ SLS can be obtained by choosing $r_{\rm c}$ in Eq. (19) as small as possible. With this choice of $\rho _{\rm s}$ and R, the mass of the resulting SLS will be slightly higher than that of the corresponding nonsingular sphere because the inner singular density region contributes with a small mass excess (see Fig. 1). For a $1~M_{\odot}$ SLS, this excess is only about 1% of the total mass of the nonsingular sphere (see Sect. 4.5).

4 Numerical methods and tests

4.1 A Lagrangian-remap code

For the calculations of this paper we have used a 1D hydrodynamics code based on an extension of the Lagrangian-remap scheme originally developed by Lufkin & Hawley (1993). In its complete version, the code solves the equations of hydrodynamics for a self-gravitating, reacting fluid including the effects of viscosity, thermal conduction, cooling and heating, generation of chemical species and net reaction rate. Here we shall only briefly describe the methods employed in a simplified version of the code for solving Eqs. (4)-(6), coupled to the logatropic EOS (1).

The code solves Eqs. (4) and (5) written in Lagrangian integral form through the use of second-order, finite-difference (FD) methods on a staggered spherical mesh in which scalar quantities, such as the density and pressure, are assigned to the volume centers of concentric shells and vector quantities, such as the velocity, are centered at the surfaces between adjacent shells. Temporal second-order accuracy is also enforced by solving the Lagrangian equations in a predictor-corrector fashion. For problems involving more complicated physics, both the predictor and corrector steps rely on a multistep procedure to advance the solution. This amounts to operationally splitting the source contributions in the momentum and energy equations so that the final update after a given time step involves a sequence of separate substeps. For the present case, where the only sources in Eq. (5) are the pressure and gravitational forces, the predictor-corrector scheme reduces to a single-step procedure for solving the equation

$$\frac{\rm d}{{\rm d}t}\int _{\Delta {\bar V}}\rho v{\rm d}V= ... ...hi}{\partial r} +\frac{\partial p}{\partial r}\right){\rm d}V,$$ (21)

where the spatial integration is over the volume $\Delta {\bar V}$of the shifted vector shells. The continuity Eq. (4) reduces to the mass conservation statement

$$\frac{\rm d}{{\rm d}t}\int _{\Delta V}\rho {\rm d}V=0,$$ (22)

where now the integration is taken over the volume $\Delta V$ of scalar shells.

In this way, Eq. (21) is solved explicitly by evaluating the fluid acceleration according to

$$a^{n}_{i-1/2}=-\left(\frac{\partial\Phi}{\partial r}\right)^{... ...i-1/2}} \left(\frac{\partial p}{\partial r}\right)^{n}_{i-1/2}$$ (23)

so that new Lagrangian radii and velocities are obtained using the following simple scheme

$$r^{n}_{i-1/2}+\Delta t\left(v^{n}_{i-1/2}+ \frac{1}{2}\Delta ta^{n}_{i-1/2}\right),$$

$$v^{n}_{i-1/2}+\Delta ta^{n}_{i-1/2},$$ (24)

where the superscripts denote the time level ( $t^{n+1}=t^{n}+\Delta t$) and the subscripts are used to label the shell interfaces defined at half the distance between consecutive zone-centered positions, i.e., $r_{i-1/2}=({\bar r}_{i}~+~{\bar r}_{i-1})/2$. Because of the assumption of spherical symmetry, the Poisson Eq. (6) can be integrated once to give

$$\left(\frac{\partial\Phi}{\partial r}\right)_{i-1/2}= \frac{GM_{i-1/2}}{r^{2}_{i-1/2}},$$ (25)

where M_i-1/2 is the mass interior to r_i-1/2. In order to accommodate the logatropic EOS (1), we have modified the calculation of the pressure gradient term in Eq. (23) following the method suggested by RPW, which relies on the observation that the radial derivative of the logatropic pressure can be expressed as

$$\frac{1}{\rho}\frac{\partial p}{\partial r}=-Ap_{\rm c} \frac{\partial}{\partial r}\left(\frac{1}{\rho}\right),$$ (26)

which in FD form becomes

$$\frac{1}{\rho ^{n}_{i-1/2}}\left(\frac{\partial p}{\partial r... ...eft(\frac{1}{\rho ^{n}_{i}}-\frac{1}{\rho ^{n}_{i-1}}\right) ,$$ (27)

where $\Delta r_{i-1/2}=r_{i}-r_{i-1}$ is the interdistance between consecutive volume shell-centers, where scalar variables are defined.

An identical sequence to (24) is employed in the corrector part with the only change being that aⁿ_i-1/2 is now replaced with its time-centered value a^n+1/2_i-1/2. Because of mass conservation during the Lagrangian step, time-centered values of the potential gradients follows directly from use of the time-centered radii r^n+1/2_i-1/2=(rⁿ_i-1/2+rⁿ⁺¹_i-1/2)/2 in Eq. (25). Similarly, time-centered values of the pressure forces are obtained from the time-centered densities using $\rho ^{n+1/2}_{i}=(\rho ^{n}_{i}\Delta V^{n}_{i})/\Delta V^{n+1/2}_{i}$, where $\Delta V$ is the volume of a particular concentric shell.

After completion of the Lagrangian step, the solution is remapped back onto an Eulerian grid, which can be either fixed or moving, by assuming piecewise-linear representations of the density and velocity to preserve the second-order accuracy of the Lagrangian solution. The remap equations for the density and velocity are constructed directly from the laws of mass and momentum conservation. A detailed account of the remap formulation is given by Lufkin & Hawley (1993) and so it will not be repeated here. Since the Eulerian grid is allowed to follow the Lagrangian distortions, the remap procedure is adaptive in nature. In order to guarantee entropy generation in the presence of shocks and spread them over a fixed number of zones, the code employs a scalar formulation for the artificial viscosity. However, for the present collapse calculations the artificial viscosity has never been used to mediate the arising shocks.

In addition to this spherical-coordinate based version of the code, a Cartesian-like version has also been written. Both versions have been tested on a variety of problems including the standard Riemann shock tube (Sod 1978), the propagation of a shock wave in both planar and spherical geometry (Noh 1987), the collapse of a pressureless sphere and of a protostellar gas cloud to stellar densities (Sigalotti & Klapp 2001) and the implosion of a neutron star (May & White 1966). In particular, the results of Noh's problem and the collapse of a neutron star are described in Sigalotti & Mendoza-Briceño (2002), where the code has been applied to hydrodynamic models of solar coronal loops.

4.2 Sink cell treatment

Under the assumption of spherical symmetry, the velocity goes to zero at the origin (r=0), which then allows the use of a reflecting inner boundary condition. Starting with a nonsingular density profile, the collapse towards singularity formation is followed with the Lagrangian version of the code, i.e., without repositioning the fluid elements after a time step. This allows to better resolve the fluid's behavior near the origin. However, close to singularity formation, the central core regions develop extremely high densities and velocities, causing a computational block when the Lagrangian version of the code is used. At this time, the innermost shells shrink to practically zero width, making the time step to become extremely small. Even if the Eulerian remap procedure is activated before singularity formation, this will not be of help because the infalling material may shock as it encounters the singularity. When this happens, high velocities develop behind the shock also leading to a drastic reduction of the time step. To avoid these problems, an inflow/outflow boundary condition is enforced away from the origin at the time of singularity formation. This is done by lopping off the smallest grid shells into a sink cell. This boundary condition, modelled after that of Boss & Black (1982), is here implemented following a treatment very similar to that of Foster & Chevalier (1993).

The use of a central sink cell effectively removes from the calculation the details of the flow around the singularity and isolates it from the rest of the grid. In this way, we can maintain a reasonably large time step and follow the accretion of the core envelope over its own collapse timescale. Because the flow across the sink cell boundary is supersonic off the grid, the inflow/outflow does not affect the calculation. Once the sink cell is activated, the subsequent collapse is followed using the Eulerian version of the code. The activation is done automatically in the course of the calculation by monitoring the size of the time step. Roughly independently of the initial spatial resolution, the time of singularity formation is achieved when the ratio of the current time step to the initial one becomes less than 10^-3. At this point, the sink cell is activated. For runs with 400 radial points, the central regions are highly resolved at the time of singularity formation so that the sink cell usually occupies a very small spherical volume with the first 10 - 20 computational shells merging into the sink region.

The mass which enters the sink is assumed to be condensed into a central point mass located at the origin. This mass will no longer interact hydrodynamically with the rest of the grid, but only gravitationally via a point mass potential (Boss & Black 1982; RPW). In this way, the gravitational acceleration term in Eq. (25) is modified according to

$$\left(\frac{\partial\Phi}{\partial r}\right)_{i-1/2}= \frac{G}{r^{2}_{i-1/2}}\left(M_{\rm acc}+M_{i-1/2}\right) ,$$ (28)

where $M_{\rm acc}$ is the mass within the sink region. Its time evolution is calculated by subtracting to the rest of the grid the mass that has crossed the sink boundary after a time step and adding it to the central point mass. For convenience, the density of material within the sink is set to a spatially uniform value and reset after each time step to the density of the first adjacent grid shell. Since the sink density determined this way is not related to the true density in that region, the pressure gradients at the sink surface are determined by linear extrapolation from the pressure gradients in the adjacent two shell surfaces. This inflow/outflow boundary condition effectively separates hydrodynamically the mass within the sink cell (protostar) from the rest of the collapsing core (envelope), while maintaining the gravitational interaction between the two regions as the dominant coupling process (Boss & Black 1982).

Activation of the sink cell before singularity formation makes no difference in the results provided that a supersonic inflow has already been established in the central regions. The sink cell can then be viewed as a useful numerical tool to alleviate the problem of increasingly small time steps, permitting in this way to record the mass accretion history of the central protostar.

4.3 Outer boundary conditions

The outer core surface is handled by means of a constant-pressure boundary condition, in which the pressure there is kept in time at a constant value equal to the initial truncation pressure $p_{\rm s}$. This boundary condition was developed by Boss & Black (1982) and thereafter used widely in many collapse calculations (e.g., Foster & Chevalier 1993; RPW). Although the velocity immediately outside the external surface is zeroed such that there is no mass entering the grid, the velocity at the core surface can vary in such a way as to keep the pressure there at a constant value. In this way, the outer boundary is strictly outflow allowing mass to leave the grid as needed. For a typical calculation, the amount of mass that leaves the grid this way is negligible and so the total mass is essentially conserved. Physically, this is consistent with having a growing central protostar accreting from a finite mass reservoir.

4.4 Stable logatropic equilibria

The ability of the code to maintain the equilibrium of a stable logatropic sphere was tested for both the Lagrangian and the Eulerian mode. Referring to the equilibrium models of Table 1, we followed the evolution of subcritical logatropes ( R/r₀<24.367) of mass 1, 10, 50 and 90 $M_{\odot}$ for more than 100 $t_{\rm ff}$. The calculations were made using initial grid resolutions corresponding to 200 and 400 uniformly spaced radial points. The $1~M_{\odot}$ logatrope with 200 zones oscillated about equilibrium with deviations from the initial density profile of $\sim$0.02% over periods of about 20  $t_{\rm ff}$ and velocities never exceeding $1.0\times 10^{-5}c_{\rm s}$. The 400 zones calculation produced oscillations in the density of $\sim$0.009% and maximum velocities less than $1.0\times 10^{-6}c_{\rm s}$. For comparison RPW, reported velocities less than $\sim$ $3.0\times 10^{-5}c_{\rm s}$ for the same model calculation at roughly comparable grid resolutions.

The code was also able to maintain the equilibrium structure of more massive logatropes for both types of calculations (Lagrangian and Eulerian). The 10, 50 and 90 $M_{\odot}$ cores with 200 radial zones, experienced oscillations about the equilibrium density with deviations of $\sim$0.06, 0.1 and 0.35%, respectively. Much lower deviations were obtained in the high-resolution (400 zones) calculations. For the more critical (90 $M_{\odot}$) case, the maximum velocities achieved never surpassed peaks of $\sim$ $4.0\times 10^{-2}c_{\rm s}$. These results show that the code is making a good job in maintaining for long times the equilibrium of stable logatropes even for masses close to the critical value.

4.5 Singular collapse of a 1 M$_{\odot }$ logatrope

In this subsection, we describe the results obtained for the collapse of a singular $1~M_{\odot}$ logatrope using the Eulerian version of the code. This test calculation provides direct comparison with the theoretical self-similar solution of MP97 and the numerical calculations of RPW.

In contrast with the nonsingular collapse cases, the singular collapse is started in the presence of a central sink cell. The initial profile is mapped onto a spherical grid of 200 uniformly spaced radial points as described in Sect. 3.2. The sink cell is chosen such that its boundary is made to coincide with the second shell interface from the origin, of radius $r_{\rm sink}=1.75\times 10^{14}$ cm and corresponding to a small fraction ( $r_{\rm sink}/R\approx 6\times 10^{-4}$) of the total truncated core radius. We assign to the sink region an initial density given by Eq. (19) with $r_{\rm c}=r_{\rm sink}$ and $p_{\rm c}$ as defined by Eq. (20). As in RPW, an initial 25% enhancement is applied to the sink density so that the core is slightly out of hydrostatic balance and begins to collapse slowly. With this choice of the initial parameters, the mass at t=0 within the sink boundary is $M_{\rm acc}=4.48\times 10^{-7}M_{\rm tot}$; a very small fraction of the total core mass ( $M_{\rm tot}\approx 1.0143~M_{\odot}$). Compared to a $1~M_{\odot}$ nonsingular logatrope, the SLS will be in excess of $\approx$0.0143 $M_{\odot}$ due to the inner singular density region.

Figure 2 shows the radial velocity profiles at different times during the initial collapse phase. The solid lines in the figure correspond to the analytic MP97's expansion-wave solution. The sequence of times was chosen to facilitate comparison with the results of RPW (their Fig. 5). The maximum deviations from the SLS collapse solutions of MP97 occur at early times (0.35 and 0.5 $t_{\rm ff}$) as shown in Fig. 2. This is a reflection of the initial density distribution near the center not being exactly singular as in the analytic case. Later on in the collapse the numerical and analytic solutions achieve a much better agreement. Also note that at $1.8t_{\rm ff}$, when the expansion wave reaches the surface of the collapsing core, the numerical results agree with the analytic predictions. Figure 3a depicts the density profiles throughout the core size for the same times of Figs. 2 and 3b shows a central amplification of the same data. We see that the numerical solution closely follows the analytic predictions of MP97 at all radii up to $1.8t_{\rm ff}$. After this time, the expansion wave leaves the core and we can no longer follow its behavior. This is more clearly seen in Fig. 2, where at $2.5t_{\rm ff}$ the numerical solution matches the analytic velocity profile only within the truncation radius. However, we can follow the subsequent collapse of the core and track the central protostellar accretion for sufficiently long times. After $1.8t_{\rm ff}$, the outer core regions become more and more depleted during further collapse. Meanwhile the inner regions achieve progressively higher densities as more matter is being condensed by the collapse itself. This behavior can be interpreted as a result of the expansion wave being reflected and then driven back into the core as a compression wave.

Figure 2: Radial velocity profiles (filled circles) compared with the analytic expansion-wave solution of MP97 (solid lines) for the collapse of a singular $1~M_{\odot}$ logatrope, as obtained using the Eulerian version of the code with 200 radial points. The sequence of times, starting from the lowest curves, is $t/t_{\rm ff}=0.35$, 0.5, 1.0, 1.25, 1.8 and 2.5. The dashed vertical line marks the truncation radius of the core.

Open with DEXTER

Figure 3: a) Density profiles through the entire radius of the core compared with the analytic expansion-wave solution of MP97 (solid lines) for the collapse of a singular $1~M_{\odot}$ logatrope, as obtained using the Eulerian version of the code with 200 radial points. In b) an amplification of the central core regions is shown for the same data. A sequence of times is presented, corresponding to $t/t_{\rm ff}=0.35$(filled cicles), 0.5 (plus signs), 1.0 (asterisks), 1.25 (open circles) and 1.8 (crosses). The dashed vertical line in a) marks the truncation radius of the core.

Open with DEXTER

	Figure 4: Time evolution of the central mass accretion for the collapse of a singular $1~M_{\odot}$ logatrope (solid line) compared with a fit to the analytic self-similar prediction of MP97 (dashed line) as obtained using Eq. (29). Each curve refers to the logarithm of the central protostellar mass in units of the total core mass.
Open with DEXTER

	Figure 5: Time evolution of the central mass accretion and central mass accretion rate for the collapse of a singular $1~M_{\odot}$ logatrope. The solid lines show the variation of $M_{\rm acc}/M_{\rm tot}$(thin line) and d $M_{\rm acc}/{\rm d}t$ in units of $M_{\rm tot}/t_{\rm ff}$ (thick line). The dashed line is the fit to the analytical prediction of MP97 given by Eq. (29).
Open with DEXTER

Figure 4 shows the time evolution of the mass contained within the sink cell, $M_{\rm acc}(t)$, in units of the total core mass. The evolution of the $1~M_{\odot}$ singular collapse was continued until about $8.41t_{\rm ff}$ by which time more than 90% of the available core mass has accreted onto the central protostar. Also shown in the figure is a two-parameter fit to the anlytical solution given by

$$\frac{M_{\rm acc}(t)}{M_{\rm tot}}=C_{1}\left(\frac{t}{t_{\rm ff}} \right)^{4}+C_{2},$$ (29)

where $C_{1}=\pi ^{4}m_{0}/(16\sqrt{2})$ is a scaling factor given by the constant coefficient in Eq. (12). The second fitting parameter $C_{2}=M_{\rm acc}(t=0)/M_{\rm tot}\approx 4.48\times 10^{-7}$ accounts for the fact that the singularity which is numerically represented by the sink cell has initially a small, but finite mass. The use of Eq. (29) tells us in which measure and how long the numerical solution follows the analytic prediction $M_{\rm acc}\propto t^{4}$. By $1.8t_{\rm ff}$, when the expansion wave leaves the core $M_{\rm acc}/M_{\rm tot}\approx 0.034$, which is very close to the analytically predicted value of 0.03 (see Sect. 2.1). Accretion onto the central protostar follows the t⁴ behavior until approximately $3.23t_{\rm ff}$ by which time about 32% of the total mass has fallen into the central sink. For comparison, RPW found that their solution followed the t⁴ behavior until $\sim$ $3.1t_{\rm ff}$, at which time the central protostar has accreted 35% of the total core mass. Figure 5 also shows the time variation of $M_{\rm acc}(t)$ along with that of the central accretion rate ${\dot M}_{\rm acc}(t)$. We see that all curves, including the fit given by Eq. (29), intersects at $\sim$ $3.23t_{\rm ff}$. After this time, the accretion rate increases a little bit more reaching a maximum value ${\dot M}_{\rm acc}\approx 0.323~M_{\rm tot}/t_{\rm ff}$by $3.295t_{\rm ff}$, and then slows down for the remainder of the evolution. As was outlined by RPW, this behavior is expected because the accreting central protostar is drawing mass from a finite reservoir. Once the expansion wave reaches the outer boundary, it will not set new mass into collapse and so the flow of inward-moving material cannot be maintained at the same rate. Further, on comparing with the results of RPW, we see that Fig. 5 reproduces very well the results in their Fig. 8. By extending the validity of the self-similar solution beyond $1.8t_{\rm ff}$, MP97 predicted that 100% of the mass in a SLS collapse would be accreted by $4.3t_{\rm ff}$. Because of the decline of the accretion rate, the present calculations show that more than twice this time is actually required by the central protostar to accrete $\sim$90% of the total mass. This result also follows the predictions made by RPW.

Figure 6: a) Density and b) radial velocity profiles in the collapse of a perturbed nonsingular $1~M_{\odot}$ core with initial parameters as given in Table 1. The sequence of times, starting from the lowest curves, is $t/t_{\rm ff}=-1.092$ (initial density profile) in a), and -0.189, $-4.13\times 10^{-2}$, $-7.40\times 10^{-3}$, $-2.68\times 10^{-3}$, $-1.17\times 10^{-3}$ and 0 (time of singularity formation) in a) and b), respectively. In a), the inclined dashed line indicates a $\rho \propto r^{-3/2}$variation and the vertical one marks the scale radius r0 of the core. The better inner resolution shown in b) is an effect of the density being staggered with respect to the velocity in the present calculations.

Open with DEXTER

5 Nonsingular collapse models: Isolated star formation

In this section, we describe the results of the collapse of nonsingular spheres, starting with initial parameters as listed in Table 1. We define a sequence of seven model calculations with both subcritical (1, 10, 50, 60, 80 and 90 $M_{\odot}$) and critical ($\approx$92.05 $M_{\odot}$) masses. All these models have the same fiducial truncation pressure $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K and so they may be representative of star-forming cloud cores in an isolated environment. If we increase the equilibrium density profile $\rho _{0}(r)$ of a subcritical sphere of mass M by a small fraction $\delta$ such that $\rho (r)\to\rho _{0}(r)(1+\delta)$, its mass and pressure will also be enhanced by the same fractional amount. A simple linear analysis can be done to demonstrate that this small mass excess will bring the sphere out of equilibrium, causing a self-gravitational acceleration, given by $-GM\delta /r^{2}$, at every radius within the sphere. For all cases, the collapse of the nonsingular sphere is initiated by adding a 5% density enhancement ( $\delta =0.05$) to the equilibrium profile throughout the core, which is the same perturbation employed by RPW in their calculations.

5.1 Collapse of a 1 M$_{\odot }$ core

We separate the discussion of the $1~M_{\odot}$ collapse from that of the more massive cases because: (a) this model provides a further direct comparison with the results of RPW and (b) it is used to check the effects of increasing the initial spatial resolution from 200 to 400 radial zones.

Figure 7: a) Time evolution of the central mass accretion and b) central mass accretion rate in the collapse of a nonsingular $1~M_{\odot}$ core with initial parameters as given in Table 1 at high (400 radial zones, thick solid line) and low (200 radial zones, thin solid line) resolution. In a) the vertical dashed line marks the time at which the density profile becomes singular and the time is represented in terms of $T=t-t_{\rm sing}$, where $t_{\rm sing}=1.092t_{\rm ff}$. To facilitate comparison with the results of RPW, the vertical dashed lines in b) indicate the transition from the protostellar collapse phase (PC) to the Class 0 stage (t=0) and from the Class 0 to the Class I stage ( $t=7.0\times 10^{5}$ yr), respectively. Note that during the whole accretion phase the high- and low-resolution runs produce essentially the same solution.

Open with DEXTER

The details of the collapse preceding the time of singularity formation are shown in Fig. 6 for the 400 zones case. This figure depicts the shape of both the density (Fig. 6a) and radial infall velocity (Fig. 6b) profiles for a sequence of times, including that of singularity formation ( $t_{\rm sing}\approx 1.092t_{\rm ff}$). Here, we adopt the same convention of RPW, in which t=0 marks the instant at which the density profile becomes singular. In this way, all intermediate adjustments from the nonsingular to the singular configuration occur at negative times (t<0). The core collapses to progressively higher densities through a sequence of profiles with decreasing inner, flat-topped regions. The evolution of the density is qualitatively similar in form to the collapse of the Bonnor-Ebert sphere (Foster & Chevalier 1993). In the logatropic case, however, the flat-topped region approaches an approximately r^-3/2 density profile by t=0 as shown in Fig. 6a. This singular region extends in radius up to $r\approx r_{0}$, while for r>r₀ the profile matches the $\rho\propto r^{-1}$ power-law. The co-existence of two distinct density gradients at the time of singularity formation is not seen in the collapse of the Bonnor-Ebert sphere, where at t=0 the inner, flat region approaches the same r^-2 variation of the rest of the sphere (Foster & Chevalier 1993). The velocity profiles also evolve self-similarly with the peak in the infall velocity moving inwards (Fig. 6b). The run with 200 zones produced essentially the same results displayed in Fig. 6, except that within a radius of $\approx$ $3.0\times 10^{14}$ cm during the rapid transition to the SLS-like collapse, the velocities were a factor of $\sim$2 lower than in the 400 zones calculation. This clearly implies that the lack of sufficient spatial resolution in the vecinity of the singularity may lead to an underestimate of the supersonic character of the flow in pure hydrodynamic calculations.

A central sink cell is activated at the precise time of singularity formation and thereafter the evolution is continued using the Eulerian version of the code. Since prior to t=0, the calculation is made in Lagrangian form we monitor the temporal variation of $M_{\rm acc}(t)$only after the sink cell is activated (t>0). Figure 7 shows the mass accretion profile for the nonsingular $1~M_{\odot}$ collapse (Fig. 7a) along with the form of the accretion rate as expressed in physical units (Fig. 7b). These figures compare very well with the results of RPW (see their Figs. 7 and 11b, respectively). We see that the accretion rate increases steeply just after singularity formation. By the time, the sink contains $\sim$2% of the total core mass, $M_{\rm acc}(t)$grows with time more slowly. During this period ( $T=t-t_{\rm sing}\leq 3.0t_{\rm ff}$), the sink cell accretes more than 60% of the total available mass. Thereafter, the central accretion enters a slow phase of relatively steady growth in which the remaining 40% core mass is accreted. This last phase occurs on a much longer timescale. The calculation was terminated at $T=t-t_{\rm sing}\approx 8.5t_{\rm ff}$ ($\sim$ $3.6\times 10^{6}$ yr), when $\sim$99% of the total mass has been accreted by the central protostar.

5.2 Collapse of massive cores

RPW studied the collapse of nonsingular logatropic spheres for masses $\leq$5 $M_{\odot}$. They found that by increasing the subcritical mass from $1~M_{\odot}$ to $5~M_{\odot}$, the adjustment of the nonsingular configuration to a singular density profile occurs with a decreasing fractional amount of the infalling supersonic mass at t=0. They also extrapolated this trend to higher mass collapses, arguing that a critical nonsingular logatrope may approach the singular profile in an entirely subsonic manner in accordance with the theoretical predictions of MP96 and MP97. With the purpose of verifying these trends, we have performed high-resolved calculations of massive subcritical (10-90 $~M_{\odot}$) and critical ($\approx$92.05 $M_{\odot}$) spheres (see Table 1). Since the primary goal of these calculations is to study the purely hydrodynamical behavior of the collapse of massive logatropes, we obviate the effects of magnetic fields and radiative transfer (see Sect. 7 for a discussion).

 

 

Table 3: Times of singularity formation and infall mass properties for A=0.2 collapsing nonsingular logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K.

$\frac{M_{\rm tot}}{M{\odot}}$	$t_{\rm sing}$	$t_{\rm sing}$	$t_{\rm sing}$	% of $M_{\rm sup}$	$\frac{R_{\rm sup}}{R}$
	(105 yr)	($\tau$)	( $t_{\rm ff}(r_{0})$)	at t=0	at t=0
1	5.28	1.436	1.287	6.01	0.164
10	4.05	1.436	1.290	0.59	0.052
50	2.76	1.435	1.286	0.08	0.020
60	2.53	1.435	1.290	0.06	0.017
80	2.07	1.435	1.289	0.04	0.013
90	1.71	1.434	1.290	0.03	0.011
92.05	1.48	1.435	1.293	0.02	0.010

Figure 8: a) Density and b) radial velocity profiles during the collapse of a nonsingular critical ( $M=92.05~M_{\odot}$) core (Table 1) preceding the time of singularity formation (t=0). The sequence of times, starting from the lowest curves, is $t/t_{\rm ff}=-0.255$ (initial density profile) in a), and $-5.25~\times~ 10^{-2}$, $-2.26~\times~ 10^{-2}$, $-5.82~\times~ 10^{-3}$, $-1.94~\times~ 10^{-3}$, $-6.90~\times~ 10^{-4}$ and 0 (time of singularity formation) in a) and b), respectively. In a), the inclined dashed line indicates a $\rho \propto r^{-3/2}$ variation and the vertical one marks the scale radius r0 of the core. The better inner resolution shown in b) is an effect of the density being staggered with respect to the velocity in the present calculations.

Open with DEXTER

For all masses considered here, the initial collapse phase prior to singularity formation is qualitatively similar to that of the $1~M_{\odot}$ case described above. This is shown in Fig. 8, which depicts the time evolution of the density (Fig. 8a) and radial infall velocity (Fig. 8b) profiles for the collapse of the critical $92.05~M_{\odot}$ logatrope up to t=0. Table 3 lists for all models the time of singularity formation in years and in units of the sound crossing time $\tau =r_{0}/\sigma _{\rm c}$ interior to r₀ and the free-fall time $t_{\rm ff}(r_{0})$ of the initial flat-topped region (r<r₀), respectively, the percentage of the total core mass that undergoes supersonic infall at t=0 and the fractional radius containing that mass. As the mass of the logatrope is increased, the timescale of singularity formation decreases. This is understandable because cores of higher mass exhibit larger central condensations, and consequently, higher average densities. The results in Table 3 also confirm the finding of RPW that the time elapsed between the initiation of collapse and the formation of the density singularity obeys the approximate scaling $t_{\rm sing}\approx 1.435r_{0}/\sigma _{\rm c}\approx 1.29t_{\rm ff}(r_{0})$ for all masses, including the critical one. Most importantly, the amount of mass falling supersonically at t=0 as well as its volume drop by increasing the mass of the logatropic core. About 6% of the total core mass collapses supersonically in a $1~M_{\odot}$ core, while only 0.02% behaves this manner in a critical logatrope. Although this result is in agreement with the trends found by RPW, the approach to the singularity is never entirely subsonic at least in a pure hydrodynamical scenario. The Mach number of the supersonic flow also decreases as R/r₀, or equivalently $M_{\rm tot}$, increases towards the critical value. The calculations predict maximum Mach numbers between $\sim$26 (critical mass) and $\sim$116 ( $1~M_{\odot}$), which are by far much higher than those reported by RPW. An explanation for this is the much higher central resolution achieved in the present models, which allows to go much deeper into the singular region than the calculations of RPW did. However, the presence of such extremely high velocities in the supersonic region may be indicative of the physical limitations of purely hydrodynamical models. For instance, under the action of a mean magnetic field these velocities should be substantially lowered. In addition, higher densities in more massive cores lead to greater optical depths and hence to increasingly higher pressures at t=0, which should also work on reducing the Mach number to physically realistic values. Future simulations of the collapse of massive cores must therefore include the effects of magnetic fields and radiative transfer to correctly describe the collapse of logatropic configurations.

	Figure 9: a) Central mass accretion and b) central mass accretion rate expressed in physical units as a function of time for all of the nonsingular collapses of Sect. 5. In all cases, the evolution is shown up to the point where 99% of the total core mass has been accreted by the central protostar.
Open with DEXTER

The main goal of pursuing the evolution farther in time after singularity formation is to track the central protostellar accretion history. Figure 9 displays the time evolution of the central mass accretion (Fig. 9a) and central mass accretion rate (Fig. 9b) in physical units for all models. In all cases, the evolution is shown up to the point where $\sim$99% of the total mass has been accreted. The $1~M_{\odot}$core is the first to complete its accretion phase. The accretion timescale then increases with increasing mass of the core up to about $50~M_{\odot}$. For masses higher than this, the trend reverses and the accretion timescale decreases as the total core mass approaches criticality. This dependence of the overall accretion lifetime on core mass is expected because the mean free-fall time has also a similar dependence on total mass (see Table 1). This is in contrast with the singular logatropic collapse theory of MP97, where more massive stars take longer to form in the logatrope. The results imply that stars of mass $1~M_{\odot}\leq M_{\star}\leq 92~M_{\odot}$ all tend to form within 3.6- $6.6\times 10^{6}$ yr. This small spread in star-formation times is due to the t⁴ variation of the mass accretion and also to the weak dependence of the mean free-fall time on mass. In terms of Fig. 9b, the accretion phase in a logatrope can be subdivided into three main stages. The first one, which marks the transition from the collapse phase into an SLS-like accretion phase, is the shortest one and is characterized by an abrupt increase of the mass accretion rate immediately after singularity formation. Note that cores of higher mass experience a steeper accretion rate during this first stage. During the second stage, the accretion rate grows more slowly until it reaches a maximum value. At this epoch, about 40% of the total mass has been accreted by the central protostar. After this point, a third stage starts in which the accretion rate declines for the remainder of the evolution until 100% of the total core mass is condensed into the central object. The higher is the mass of the core, the larger is the size of the accretion rate. In particular, our models predict peak values of ${\dot M}_{\rm acc}$ as high as $\sim$5- $6\times 10^{-5}~M_{\odot}$ yr^-1 for cores close to the critical mass. These are about two orders of magnitude larger than the peak value achieved by the highly subcritical $1~M_{\odot}$ core ($\sim$ $8\times 10^{-7}~M_{\odot}$ yr^-1).

	Figure 10: Central mass accretion rate as a function of the central accreted mass in dimensionless units for all of the nonsingular collapses of Sect. 5.
Open with DEXTER

The dependence of the central mass accretion rate on the central accreted mass in dimensionless units is shown in Fig. 10 for all models. When expressed in dimensionless units, the accretion rate of the $1~M_{\odot}$ core looks higher than that of the more massive cores. This trend is expected because for the $1~M_{\odot}$ core to be the first to accrete 99% of its total mass (see Fig. 9a), it must effectively accrete a larger fraction of its mass compared to the more massive cores. For cores of mass $\geq$50 $M_{\odot}$ the curves in Fig. 10 overlap, with all of these cores showing a lower accretion rate compared to the 1 and $10~M_{\odot}$ cores. A feature common to all curves is a maximum when the central protostar has accreted about 40% of the available mass. For higher fractions of the accreted mass, the accretion rate declines in all cases. When the accretion rate is plotted as a function of time in dimensionless units, the variation of the accretion rate has a form very similar to that shown in Fig. 10. This suggests an approximate linear dependence of the accreted mass on time. Such rough linear variation is effectively observed when $M_{\rm acc}/M_{\rm tot}$ is plotted as a function of $t/t_{\rm ff}$ and it is maintained until the accretion rate starts to decline. This is in stark contrast with the SLS case, where $M_{\rm acc}\propto t^{4}$.

6 Nonsingular collapse models: Clustered star formation

So far we have scaled our collapse models to central temperatures and truncation pressures representative of isolated star-forming regions. We now describe the results obtained for a new set of nonsingular collapse models in which the surface pressure is chosen to be $p_{\rm s}=1.0~\times~ 10^{7}k$ cm^-3 K, in keeping with observations of clustered star-forming regions such as $\rho$ Oph (André et al. 2000). The initial parameters for these model calculations are listed in Table 2. Here, we consider spheres with 0.5, 1, 2, 5, 8 (subcritical cases) and $10.5~M_{\odot}$ (critical case). As before, in each case the collapse is initiated by adding throughout the core a 5% overdensity.

For all models, the initial collapse phase towards the singular density profile proceeds as in the previous cases. Table 4 lists the times of singularity formation and the properties of the infalling mass at t=0. An inspection of this table shows that the data along the columns follow essentially the same trends seen in Table 3. In view of this and of the qualitative similar behavior of the models with those described in the preceding section, we shall discuss the results in terms of the main differences between both sets of calculations.

One effect of increasing the surface pressure is to shorten the time of collapse to singularity formation. Subcritical cores of equal mass and varying $p_{\rm s}$ will reach the singular state more rapidly for larger values of $p_{\rm s}$. In addition, for highly subcritical cores, the percentage of mass collapsing supersonically at t=0 decreases as $p_{\rm s}$ is increased, while for critical cores the relative amount of supersonic mass is independent of $p_{\rm s}$. A similar trend is also observed for the fractional volume occupied by the supersonic flow. Furthermore, increasing $p_{\rm s}$ also results in supersonic inflow of much lower Mach numbers, with maximum values ranging from $\sim$19 ( $10.5~M_{\odot}$) to $\sim$38 ( $0.5~M_{\odot}$). Thus, in a purely hydrodynamical scenario, the general effect of increasing the bounding pressure in the collapse of a logatrope is to produce a less supersonic approach to the singularity.

 

 

Table 4: Times of singularity formation and infall mass properties for A=0.2 collapsing nonsingular logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K.

$\frac{M_{\rm tot}}{M_{\odot}}$	$t_{\rm sing}$	$t_{\rm sing}$	$t_{\rm sing}$	% of $M_{\rm sup}$	$\frac{R_{\rm sup}}{R}$
	(104 yr)	($\tau$)	( $t_{\rm ff}(r_{0})$)	at t=0	at t=0
0.5	5.17	1.436	1.287	1.53	0.085
1	4.71	1.435	1.285	0.70	0.058
2	4.18	1.434	1.289	0.29	0.037
5	3.30	1.434	1.285	0.09	0.021
8	2.63	1.433	1.291	0.05	0.015
10.5	1.68	1.433	1.288	0.02	0.010

	Figure 11: a) Central mass accretion and b) central mass accretion rate expressed in physical units as a function of time for all of the nonsingular collapses of Sect. 6. In all cases, the evolution is shown up to the point where 99% of the total core mass has been accreted by the central protostar.
Open with DEXTER

For all cases, the accretion history is displayed in Fig. 11, where the time evolution of the central mass accretion (Fig. 11a) and the central mass accretion rate (Fig. 11b) are depicted in terms of physical units. As before, the accretion has been followed to the point where $\sim$99% of the total core mass has fallen into the sink cell. It is evident from Fig. 11a that the accretion timescale is longer for the $5~M_{\odot}$ collapse case. Starting from this mass, the accretion timescale shortens towards increasing and decreasing masses. Again, this is a consequence of the mean free-fall time being longer for subcritical masses around $5~M_{\odot}$ (see Table 2). For fixed A, increasing the fiducial value of $p_{\rm s}$ has little effect on the size of the accretion rate in critical logatropes. This can be seen by comparing Figs. 9b and 11b, where the accretion rate peaks at $\sim$ $6\times 10^{-5}~M_{\odot}$ yr^-1 for both the critical 10.5 and $92.05~M_{\odot}$ logatropes. However, the same is not true for subcritical logatropes sharing the same mass. For instance, a $1~M_{\odot}$ core with $p_{\rm s}=3.5\times 10^{5}k$ cm^-3 K experiences a maximum accretion rate of $\sim$ $8\times 10^{-7}~M_{\odot}$ yr^-1 compared to $\sim$ $5\times 10^{-6}~M_{\odot}$ yr^-1 achieved in a $1~M_{\odot}$core with $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K. Thus, for comparable subcritical masses, the size of the accretion rate is smaller in cores having lower values of $p_{\rm s}$. On the other hand, by comparing Figs. 9a and 11a, we see that the accretion lifetimes are considerably longer in cores with lower $p_{\rm s}$, implying that star formation in clustered regions occurs on shorter timescales ($\sim$5.4- $7.4\times 10^{5}$ yr as predicted for these models) compared to more isolated environments.

7 Comparison with observations and discussion

We have presented high-resolved calculations of the spherically symmetric collapse and subsequent accretion of nonsingular A=0.2 logatropic cores with varied masses and values of the fiducial external bounding pressure. For all cases, the evolution has been followed to the point where nearly all of the finite core mass has fallen into the central protostar. Thus, the models are indicative of the early star-formation stages, including: (a) the protostellar collapse (PC) phase, (b) the Class 0 stage, in which the central protostar accretes $\sim$50% of the total core mass (André et al. 1993) and (c) the Class I stage, in which the protostar accretes the remainder of the final stellar mass, leading to the formation of a Class II object of mass $M_{\star}>>M_{\rm env}$. These very young stellar objects are known to be optically visible and possess emission lines characteristic of classical T-Tauri stars (Adams et al. 1987). We can therefore compare the various predicted lifetimes with those derived observationally from statistical arguments and theoretically from other proposed models.

Based on a total starless core lifetime of a few 10⁶ yr, as derived from the survey of Beichman et al. (1986), Ward-Thompson et al. (1994) estimated the lifetime of prestellar cores detected in the submm continuum to be $\sim$10⁶ yr. More recently, Lee et al. (1999) found evidence for collapse in seven out of 220 cores surveyed, implying that the PC phase lasts $\sim$3% of the total starless core phase, or $\sim$3- $10\times 10^{4}$ yr. Furthermore, Gregersen & Evans (2000) discovered six more sources undergoing collapse from the catalog of 50 starless cores listed by Beichman et al. (1986). Their results imply a lifetime for the PC phase of a few 10⁵ yr. While these estimates carry large inherent uncertainties, we note that our calculations predict lifetimes for this phase of $\sim$1.5- $5.3\times 10^{5}$ yr for the $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K nonsingular collapses (see Table 3), consistent with the estimates of Gregersen & Evans (2000). For models with $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K, the predicted lifetimes are $\sim$1.7- $5.2\times 10^{4}$ yr (see Table 4), which are more consistent with the ages suggested by Lee et al. (1999). However, the detection of much larger samples of contracting sources along with more accurate estimates of their PC lifetimes, would be required to confirm the validity of our predictions. The logatropic PC lifetimes listed in Table 3 are an order of magnitude longer than the prediction of $\sim$ $4\times 10^{4}$ yr made by Whitworth & Ward-Thompson (2001), who proposed a simple analytic model for the protostellar collapse of the isolated core L1544 based on a Plummer-like radial density profile as its initial condition.

At the end of the PC phase, the density distribution becomes singular everywhere and a central point-mass forms. This transition to a SLS-like collapse marks the beginning of the Class 0 stage and is characterized by a steep increase of the mass accretion rate. This is shown schematically in Fig. 7b for the $1~M_{\odot}$ collapse case with $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K. During this stage, the accretion rate reaches a maximum value and the protostar accretes half of the mass of the envelope on a timescale of $\sim$ $7\times 10^{5}$yr, which marks the transition to the Class I protostellar phase. Class 0 objects are good candidates of very young accreting protostars in which a hydrostatic object has already formed but not yet accreted the majority of its final mass. Most of the final stellar mass is still in the form of a dense circumstellar envelope. About 30-40 confirmed Class 0 objects have been identified up to now (André et al. 2000), mostly in regions of multiple star formation. For instance, in the $\rho$ Oph main cloud, there have been reported only two good Class 0 candidates, suggesting a Class 0 lifetime of $\sim$1- $3\times 10^{4}$ yr compared to $\sim$ $2\times 10^{5}$ yr for the Class I sources (e.g., Greene et al. 1994). It then appears on observational grounds that the lifetime of the Class 0 phase is short compared to both the PC phase and the Class I near-IR phase. Tables 5 and 6 contain the predicted lifetimes for both the Class 0 and Class I phases from the sequences of logatropic collapse calculations with $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K and $1.0\times 10^{7}k$ cm^-3 K, respectively. These values are $\sim$1-2 orders of magnitude larger than those inferred observationally for the Class 0 phase, and about an order of magnitude larger than $\sim$ $2\times 10^{5}$ yr for the Class I stage in the models with lower $p_{\rm s}$ (see Table 5). The models with $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K are more representative of cloud cores contracting in a clustered environment, and as shown in Table 6, their lifetimes for the Class I stage are comparable to the observational estimates for the $\rho$ Oph main cloud. Here, we have adopted the widely accepted convention that the borderline between Class 0 and Class I sources occurs when $M_{\star}=M_{\rm acc}\approx M_{\rm env}$, while that between Class I and Class II sources occurs when 99% of the total core mass has been delivered ( $M_{\rm env}=0.01M_{\rm tot}$), which allows for residual accretion during the (Class II) pre-main-sequence evolution (Henriksen et al. 1997). For further comparison, the analytic collapse model of Whitworth & Ward-Thompson (2001) predict that L1544 should have accumulated about 50% of its mass in $\sim$ $8\times 10^{4}$yr. In the logatropic scenario, both the Class 0 and Class I lifetimes depend weakly on the total core mass and are about an order of magnitude larger in regions of isolated star formation than in cluster-forming environments.

 

 

Table 5: Times marking the borderline between Class 0 and Class I sources ( $M_{\rm acc}=0.5M_{\rm tot}$) and between Class I and Class II sources ( $M_{\rm acc}=0.99M_{\rm tot}$) for A=0.2 collapsing logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K.

$\frac{M_{\rm tot}}{M_{\odot}}$	t(Class 0/I)	t(Class I/II)
	(106 yr)	(106 yr)
1	0.70	2.89
10	1.21	4.34
50	1.56	5.01
60	1.55	4.93
80	1.46	4.55
90	1.31	4.07
92.05	1.19	3.69

 

 

Table 6: Times marking the borderline between Class 0 and Class I sources ( $M_{\rm acc}=0.5M_{\rm tot}$) and between Class I and Class II sources ( $M_{\rm acc}=0.99M_{\rm tot}$) for A=0.2 collapsing logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K.

$\frac{M_{\rm tot}}{M_{\odot}}$	t(Class 0/I)	t(Class I/II)
	(105 yr)	(105 yr)
0.5	1.12	4.23
1	1.34	4.77
2	1.56	5.26
5	1.77	5.64
8	1.73	5.39
10.5	1.36	4.16

The timescales in Tables 5 and 6 for the transition from Class 0/I to Class I/II have been determined under the assumptions that there is no mass disruption mechanism operating during the protostellar evolution and that the final stellar mass is determined by the total core mass, which have no observational support. For instance, there is confirmed evidence that NH₃ cloud cores are much more massive than the embedded young stars formed at their centers (Jijina et al. 1999). Furthermore, it is well-known that star formation within molecular cloud cores is accompanied by some degree of outflow activity (see discussion below). Such protostellar outflows have been indicated as a mechanism to remove gas from star-forming regions (e.g., Bally et al. 1999) and from protostellar cores (e.g., Velusamy & Langer 1998; Ladd et al. 1998). Mass disruption by bipolar outflows then limits the fraction of mass that can accrete and therefore the efficiency of star formation in a protostellar core (Matzner & MacKee 2000). Thus, more accurate predictions of the Class 0 and Class I lifetimes would evidently require a solution to the problem of simultaneous infall and outflow. Most, if not all, Class 0 protostars drive highly collimated or "jet-like'' CO molecular outflows (Bachiller 1996). In contrast, the CO outflows from Class I sources are much less powerful and collimated. In particular, Bontemps et al. (1996) analyzed a set of CO(2-1) outflow data around a large sample of embedded young stellar objects, including 9 Class 0 sources and 36 Class I sources. They found that the outflow phase and the infall/accretion phase coincide for the entire sample of sources examined, thereby suggesting that outflows are directly powered by accretion. In addition, magneto-centrifugal accretion/ejection models of bipolar outflows (Ouyed & Pudritz 1997; Kudoh & Shibata 1997) predict a direct proportionality between accretion and ejection. This implies that the observed decline of outflow power must be accompanied by a corresponding decline of the mass accretion rate from the Class 0 to the Class I stage (Bontemps et al. 1996). The results of Bontemps et al. (1996) further indicate that the ejected mass ${\dot M}_{\rm jet}$ declines from $\sim$10 $^{-6}~M_{\odot}$ yr^-1 for the youngest Class 0 protostars to $\sim$ $2\times 10^{-8}~M_{\odot}$ yr^-1 for the most evolved Class I objects. Realistic jet models give the relation ${\dot M}_{\rm jet}/{\dot M}_{\rm acc}\sim 0.1 {-} 0.3$ between the accretion and ejection rates (André et al. 2000). This implies that ${\dot M}_{\rm acc}$ decreases from $\sim$0.3- $1\times 10^{-5}~M_{\odot}$ yr^-1 for Class 0 objects to $\sim$0.7- $2\times 10^{-7}~M_{\odot}$ yr^-1 for evolved Class I protostars.

 

 

Table 7: Maximum values of the mass accretion rate along with its corresponding values at the borderline between Class 0 and Class I sources ( $M_{\rm acc}=0.5M_{\rm tot}$) and between Class I and Class II sources ( $M_{\rm acc}=0.99M_{\rm tot}$) for A=0.2 collapsing logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K.

$\frac{M_{\rm tot}}{M_{\odot}}$	$\max ({\dot M}_{\rm acc})$	${\dot M}_{\rm acc}$(Class 0/I)	${\dot M}_{\rm acc}$(Class I/II)
	($M_{\odot}$ yr-1)	($M_{\odot}$ yr-1)	($M_{\odot}$ yr-1)
1	$7.8\times 10^{-7}$	$7.3\times 10^{-7}$	$2.0\times 10^{-8}$
10	$5.2\times 10^{-6}$	$5.0\times 10^{-6}$	$1.3\times 10^{-7}$
50	$2.3\times 10^{-5}$	$2.2\times 10^{-5}$	$5.6\times 10^{-7}$
60	$2.8\times 10^{-5}$	$2.7\times 10^{-5}$	$6.8\times 10^{-7}$
80	$4.1\times 10^{-5}$	$3.9\times 10^{-5}$	$9.8\times 10^{-7}$
90	$5.1\times 10^{-5}$	$4.9\times 10^{-5}$	$1.2\times 10^{-6}$
92.05	$5.8\times 10^{-5}$	$5.5\times 10^{-5}$	$1.4\times 10^{-6}$

 

 

Table 8: Maximum values of the mass accretion rate along with its corresponding values at the borderline between Class 0 and Class I sources ( $M_{\rm acc}=0.5M_{\rm tot}$) and between Class I and Class II sources ( $M_{\rm acc}=0.99M_{\rm tot}$) for A=0.2 collapsing logatropes with $T_{\rm c}=10$ K and $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K.

$\frac{M_{\rm tot}}{M_{\odot}}$	$\max ({\dot M}_{\rm acc})$	${\dot M}_{\rm acc}$(Class 0/I)	${\dot M}_{\rm acc}$(Class I/II)
	($M_{\odot}$ yr-1)	($M_{\odot}$ yr-1)	($M_{\odot}$ yr-1)
0.5	$2.6\times 10^{-6}$	$2.5\times 10^{-6}$	$6.9\times 10^{-8}$
1	$4.7\times 10^{-6}$	$4.4\times 10^{-6}$	$1.2\times 10^{-7}$
2	$8.5\times 10^{-6}$	$8.1\times 10^{-6}$	$2.2\times 10^{-7}$
5	$2.0\times 10^{-5}$	$1.9\times 10^{-5}$	$5.1\times 10^{-7}$
8	$3.4\times 10^{-5}$	$3.2\times 10^{-5}$	$8.4\times 10^{-7}$
10.5	$5.8\times 10^{-5}$	$5.5\times 10^{-5}$	$1.4\times 10^{-6}$

In Tables 7 and 8 we list the maximum values of ${\dot M}_{\rm acc}$ along with its values at the borderlines between Class 0 and Class I sources ( $M_{\rm acc}\approx 0.5M_{\rm tot}$) and between Class I and Class II sources ( $M_{\rm acc}\approx 0.99M_{\rm tot}$), as predicted by our logatropic collapse models. In all cases, ${\dot M}_{\rm acc}$ reaches a peak value when the protostar contains $\sim$40% of its final mass. We see that the predicted ${\dot M}_{\rm acc}$ agrees with the observational estimates for cores with $M_{\rm tot}<10~M_{\odot}$ in the sequence of calculations with $p_{\rm s}=1.3\times 10^{5}k$ cm^-3 K (Table 7) and for cores with $M_{\rm tot}\leq 2~M_{\odot}$ in the higher $p_{\rm s}$ sequence (Table 8). For logatropes of higher mass in both sequences, the accretion rates are higher than the observational estimates by factors of $\sim$2-7. Furthermore, the logatropic models predict values of ${\dot M}_{\rm acc}$ which are factors of $\sim$40 larger for Class 0 sources than for Class I sources, implying much higher accretion luminosities for the former than for the latter. This factor is 4 times higher than that proposed by Henriksen et al. (1997). However, both factors seem to be in apparent contradiction with observations which suggest that Class 0 sources are not significantly over-luminous compared to the more evolved Class I objects. In particular, the embedded source sample of Bontemps et al. (1996) indicates an estimated ratio of the average bolometric luminosities between Class 0 and Class I sources of $\sim$1.6. However, this estimate can be affected by the observational uncertainty of the typical luminosity of Class 0 sources due to the rarity of these objects. Recently, Eisner et al. (2002) estimated the mass ejection rate from a high-mass protostar in W51-IRS 2 due to a bipolar outflow associated with it to be $\geq$ $10^{-6}~M_{\odot}$ yr^-1, which according to the relationship ${\dot M}_{\rm jet}/{\dot M}_{\rm acc}\sim 0.1 {-} 0.3$ would correspond to an accretion rate $\geq$10 $^{-5}~M_{\odot}$ yr^-1; a value which is consistent with our predicted Class 0 lifetimes for high-mass logatropes (see Tables 7 and 8).

We have considered a set of calculations starting with initial conditions appropriate for contracting cores within a high-pressure environment ( $p_{\rm s}=1.0\times 10^{7}k$ cm^-3 K; see Table 2). Such a high surface pressure is consistent with the expected average pressure in the $\rho$ Oph main body region. That is, given the $\sim$1 pc extent of the $\rho$ Oph central region, the average number density is expected to be $n_{\rm H_{2}}\sim 3\times 10^{4}$ cm^-3 and so the corresponding average pressure is $p/k\sim 2\times 10^{7}$ cm^-3 K. Johnstone et al. (2000) identified 55 independent cores in the central region of the $\rho$ Oph molecular cloud and fitted them to Bonnor-Ebert spheres, suggesting internal temperatures of 10 - 30 K and surface pressures in the range between $\sim$10⁶ and $\sim$10⁷k cm^-3 K for the 55 cores sampled. As was argued by RPW, these surface pressures depend on the assumption that all of these cores could be well-fit by Bonnor-Ebert spheres, and hence $p_{\rm s}=1.0~\times~ 10^{7}k$ cm^-3 K cannot necessarily be considered reflective of the values that would be required to truncate nonsingular logatropic spheres in a cluster-forming environment. In the lack of convincing evidence confirming this point, the present calculations represent a step forward in studying the protostellar evolution of logatropes under conditions appropriate for star formation in clustered environments.

The primary goal of these calculations is to study the dynamical aspects of logatropic collapse, and so we have obviated the effects of rotation, magnetic fields and radiative transfer. While logatropic models are able to fit the observed line width-size data well (MP96), they do not specify the physical basis of the nonthermal motions. The coupling of the magnetic field to the neutrals seems to be necessary to keep the amplitude of the observed nonthermal motions. Estimates of the ionization in molecular cloud cores suggest that ion-neutral coupling is sufficient to support MHD waves, even in the densest cores (Myers & Khersonsky 1995). Further, the relatively shallow dependence of the kinetic pressure on radius, $p\sim r^{n}$, with 0<n<1, from the line width-size relations, combined with the evidence for similarity of magnetic and kinetic energy density in some clouds (Myers & Goodman 1988), suggest that the magnetic field must also vary slowly within dense cores. If the line width-size relations originate in nonlinear MHD turbulence, it will take substantial computational power to simulate these phenomena for a compressible, self-gravitating medium in three dimensions. The magnetic field offers a way to slow down the prestellar collapse, leading to a much gentler approach towards the singular density profile compared to a pure hydrodynamical scenario, where unreasonably high supersonic infall velocities may eventually develop in a small central region. However, it is important to stress that the logatrope already incorporates a mean magnetic field as a virial parameter and so it accounts for some measure of both magnetic and turbulent support (MP97, RPW).

The inclusion of rotation is a further step towards improving the picture of logatropic collapse. In particular, this will ensure that protostellar accretion does proceed in an anisotropic manner, i.e., through the formation of a circumstellar accretion disk. Such a process is also likely to affect the spherically symmetric timescales listed in Tables 5 and 6, which will still apply to the limiting case of very slow rotation. A consistent treatment of contracting, rotating logatropes will demand calculating appropriate equilibria from which to start the collapse. The calculation of such rotating equilibrium sequences is mandatory for future work wherein we will generalize the present approach to the two-dimensional collapse of rotating, logatropic cores. Furthermore, a careful treatment of radiative transfer is essential for realistic simulations of the formation of massive stars. Because of their high luminosities, we expect that radiative acceleration will contribute significantly to the dynamical evolution of high-mass cores. Moreover, evidence for high-luminosity far-IR sources, which are known to be suspected embedded young OB stars, shows that they have powerful bipolar outflows associated with them (e.g., Shepherd et al. 2000). Such massive outflows are probably powered by disk accretion and, similarly to their low-mass counterparts, the flow energetics appear to scale with the luminosity of the source (Richer et al. 2000). Recently, Yorke & Sonnhalter (2002) have calculated the collapse of slowly rotating, massive molecular cores, including the effects of frequency-dependent radiation transfer and grain opacity contribution. Although their initial conditions correspond to a r^-2 SIS model, they predict a time-dependent accretion which follows qualitatively the behavior shown in Figs. 9 and 11. For the particular case of a $60~M_{\odot}$core model, they report a maximum value of the accretion rate of $\sim$ $2\times 10^{-3}~M_{\odot}$ yr^-1 which is achieved after $\sim$ $8\times 10^{3}$ yr in the collapse. This value of ${\dot M}_{\rm acc}$is significantly higher compared to the predictions of our logatropic collapse models. In the same way, the timescale of $\sim$ $8\times 10^{3}$ yr is much shorter than those marking the temporal position of the ${\dot M}_{\rm acc}$ peaks in Figs. 9 and 11. Clearly these results show that radiative transfer may strongly influence the evolution of high-mass cores, and also limit the amount of mass that will end into the form of stars.

8 Concluding remarks

In this paper, we have presented solutions for the spherical collapse and subsequent accretion of nonsingular (subcritical and critical) logatropic spheres, starting with initial conditions close to hydrostatic equilibrium. Two sequences of calculations have been conducted by increasing the fiducial truncation pressure from $p_{\rm s}/k=1.3\times 10^{5}$ cm^-3 K, representative of the conditions in relatively isolated star-forming regions, to $p_{\rm s}/k=1.0\times 10^{7}$ cm^-3 K, as would be appropriate to clustered environments. The main conclusions of this study can be summarized in the following points:

(i) The initial collapse phase preceding the time (t=0) of singularity formation is qualitatively similar for all models. At t=0, the density profile becomes singular everywhere, matching an approximate r^-3/2variation for all radii within $r\approx r_{0}$. For r>r₀, however, the density follows the initial $\rho\propto r^{-1}$ power-law. Along each sequence of calculations, the time of singularity formation decreases with increasing mass of the core. For logatropic spheres of mass ranging from 1 (highly subcritical) to $\approx$92.05 $M_{\odot}$(critical mass) in the lower $p_{\rm s}$ sequence, the predicted lifetimes for the prestellar collapse phase are $\sim$1.5- $5.3\times 10^{5}$ yr, compared to $\sim$1.7- $5.2\times 10^{4}$ yr for the higher $p_{\rm s}$sequence involving masses from 0.5 to $\approx$10.5 $M_{\odot}$ (critical case).

(ii) The approach to the singular density profile never occurs in an entirely subsonic fashion, as was suggested by RPW. Referring to the lower $p_{\rm s}$ models, about 6% of the total core mass collapses supersonically in a $1~M_{\odot}$, while only $\sim$0.02% behaves this manner in a critical ($\approx$92.05 $M_{\odot}$) logatrope. In the higher $p_{\rm s}$ sequence the same trend is observed, with $\sim$0.7% of the mass now infalling supersonically at t=0 in a $1~M_{\odot}$ sphere. Thus, the effect of increasing the fiducial truncation pressure is to produce a gentler adjustment from the nonsingular to the singular collapse stage. The time elapsed between the initiation of collapse and the singularity formation obeys the approximate scaling $t_{\rm sing}\approx 1.435r_{0}/\sigma _{\rm c}\approx 1.29t_{\rm ff}(r_{0})$ (first noted by RPW) for all masses, including the critical one.

(iii) Immediately after singularity formation, a phase of vigorous accretion follows in which the accretion rate rises steeply, reaching a maximum value precisely when the central protostar has already accreted $\sim$40% of its final mass. Thereafter, it decreases steadily for the remainder of the evolution. In all cases, we have continued the calculations until the core has delivered 99% of its mass to the central protostar. Thus, the models may be indicative of the early star formation stages, including the Class 0 ( $M_{\rm acc}=0.5M_{\rm tot}$) and the Class I ( $M_{\rm acc}=0.99M_{\rm tot}$) accretion phase.

(iv) In particular, our models predict peak values of ${\dot M}_{\rm acc}$as high as $5 {-} 6\times 10^{-5}~M_{\odot}$ yr^-1 for cores close to the critical mass, independently of the fiducial value of the truncation pressure. These are from $\sim$1 to 2 orders of magnitude larger than the peak values achieved by a highly subcritical $1~M_{\odot}$ core ($\sim$ $8\times 10^{-7}~M_{\odot}$ and $\sim$ $5\times 10^{-6}~M_{\odot}$ for the lower and higher $p_{\rm s}$ cases, respectively). At the borderline indicative of the transition between Class 0 and Class I sources, the logatropic models predict ${\dot M}_{\rm acc}\sim 0.07 {-} 5.5\times 10^{-5}~M_{\odot}$ yr^-1for the lower $p_{\rm s}$ sequence and ${\dot M}_{\rm acc}\sim 0.25 {-} 5.5\times 10^{-5}~M_{\odot}$ yr^-1 for the higher $p_{\rm s}$ cases. Similarly, the corresponding values at the borderline indicative of the transition between Class I and Class II protostars are ${\dot M}_{\rm acc}\sim 0.02 {-} 1.4\times 10^{-6}~M_{\odot}$ yr^-1 and ${\dot M}_{\rm acc}\sim 0.07 {-} 1.4\times 10^{-6}~M_{\odot}$ yr^-1 for the lower and higher $p_{\rm s}$ sequences, respectively. These results imply a factor of ${\dot M}_{\rm acc}$ between Class 0 and Class I of $\sim$40, which is at least 4 times higher than that inferred observationally.

(v) For masses in the range $1~M_{\odot}\leq M_{\rm tot}\leq 92.05~M_{\odot}$(lower $p_{\rm s}$ sequence), the accretion timescale increases with mass up to about $50~M_{\odot}$. For $M>50~M_{\odot}$, the trend reverses and the accretion timescale decreases as the total mass core approaches criticality. Similarly, for the higher $p_{\rm s}$ sequence with core masses between 0.5 and $\approx$10.5 $M_{\odot}$, the accretion timescale is longer for subcritical masses around $5~M_{\odot}$. In particular, our models predict Class 0 lifetimes of $\sim$0.7- $1.6\times 10^{6}$ yr (for lower $p_{\rm s}$) and $\sim$1.1- $1.8\times 10^{5}$ yr (for higher $p_{\rm s}$) and Class I lifetimes of $\sim$2.9- $5.0\times 10^{6}$ yr (for lower $p_{\rm s}$) and $\sim$4.2- $5.6\times 10^{5}$ yr (for higher $p_{\rm s}$), implying that star formation in clustered regions occurs on shorter timescales compared to more isolated environments. The predicted lifetimes for the higher $p_{\rm s}$ sequence are comparable with the observationally inferred Class I lifetime of $\sim$ $2\times 10^{5}$ yr in the $\rho$ Oph cluster-forming cloud.

(vi) More realistic simulations of the formation of massive stars must include the effects of radiative transfer. Because of the high luminosities of young embedded massive stars, we expect that radiative acceleration will contribute significantly to the dynamical evolution of high-mass cores, probably producing larger sizes of the maximum accretion rate, shorter lifetimes and smaller final stellar masses than predicted by the present purely hydrodynamic logatropic collapse models. An improved picture of the logatropic collapse must then necessarily require the inclusion of the effects of radiative transfer along with the effects of rotation and magnetic fields.

Acknowledgements

We thank the referee for a number of very helpful comments and suggestions that have improved the quality of the manuscript. The calculations of this paper were made using the facilities of the Laboratory of Computational Physics at the IVIC Center of Physics.

References

 

• Adams, F. C., Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788 In the text NASA ADS  
• Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159 In the text NASA ADS  
• André, P., Ward-Thompson, D., & Motte, F. 1996, A&A, 314, 625 In the text NASA ADS  
• André, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122 In the text NASA ADS  
• André, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, ed. V. G. Mannings, A. P. Boss, & S. S. Russell (Tucson: Univ. Arizona Press), 59 In the text  
• Bacmann, A., André, P., Puget, J.-L., et al. 2000, A&A, 361, 555 In the text NASA ADS  
• Bachiller, R. 1996, ARA&A, 34, 111 In the text NASA ADS  
• Bally, J., Reipurth, B., Lada, C. J., & Billawala, Y. 1999, AJ, 117, 410 In the text NASA ADS  
• Beichman, C. A., Myers, P. C., Emerson, J. P., et al. 1986, ApJ, 307, 337 In the text NASA ADS  
• Bodenheimer, P., & Sweigart, A. 1968, ApJ, 152, 515 In the text NASA ADS  
• Bonnor, W. B. 1956, MNRAS, 116, 351 In the text NASA ADS  
• Bontemps, S., André, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858 In the text NASA ADS  
• Boss, A. P., & Black, D. C. 1982, ApJ, 258, 270 In the text NASA ADS  
• Caselli, P., & Myers, P. C. 1995, ApJ, 446, 665 In the text NASA ADS  
• Chandler, C. J., & Richer, J. S. 2000, ApJ, 530, 851 In the text NASA ADS  
• Colomé, C., di Francesco, J., & Harvey, P. M. 1996, ApJ, 461, 909 In the text NASA ADS  
• Ebert, R. 1955, Z. Astrophys., 37, 217 In the text NASA ADS  
• Eisner, J. A., Greenhill, L. J., Herrnstein, J. R., Moran, J. M., & Menten, K. M. 2002, ApJ, 569, 334 In the text NASA ADS  
• Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303 In the text NASA ADS  
• Fuller, G. A., & Myers, P. C. 1992, ApJ, 384, 523 In the text NASA ADS  
• Greene, T. P., Wilking, B. A., André, P., Young, E. T., & Lada, C. J. 1994, ApJ, 434, 614 In the text NASA ADS  
• Gregersen, E. M., & Evans, N. J. 2000, ApJ, 538, 260 In the text NASA ADS  
• Henning, T., Burkert, A., Launhardt, R., Leinert, Ch., & Stecklum, B. 1998, A&A, 336, 565 In the text NASA ADS  
• Henriksen, R., André, P., & Bontemps, S. 1997, A&A, 323, 549 In the text NASA ADS  
• Hunter, C. 1977, ApJ, 218, 834 In the text NASA ADS  
• Jijina, J., Myers, P. C., & Adams, F. C. 1999, ApJS, 125, 161 In the text NASA ADS  
• Johnstone, D., Wilson, C. D., Moriarty-Schieven, G., et al. 2000, ApJ, 545, 327 In the text NASA ADS  
• Kenyon, S. J., Gomez, M., Marzke, R. O., & Hartmann, L. 1994, AJ, 108, 251 In the text NASA ADS  
• Kudoh, T., & Shibata, K. 1997, ApJ, 476, 632 In the text NASA ADS  
• Ladd, E. F., Fuller, G. A., & Deane, J. R. 1998, ApJ, 495, 871 In the text NASA ADS  
• Larson, R. B. 1969, MNRAS, 145, 271 In the text NASA ADS  
• Lee, C. W., Myers, P. C., & Tafalla, M. 1999, ApJ, 526, 788 In the text NASA ADS  
• Lufkin, E. A., & Hawley, J. F. 1993, ApJS, 88, 569 In the text NASA ADS  
• McLaughlin, D. E., & Pudritz, R. E. 1996, ApJ, 469, 194 (MP96) In the text NASA ADS  
• McLaughlin, D. E., & Pudritz, R. E. 1997, ApJ, 476, 750 (MP97) In the text NASA ADS  
• Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364 In the text NASA ADS  
• May, M. M., & White, R. H. 1966, Phys. Rev., 141(4), 1232 In the text  
• Myers, P. C., & Goodman, A. A. 1988, ApJ, 329, 392 In the text NASA ADS  
• Myers, P. C., & Khersonsky, V. K. 1995, ApJ, 442, 186 In the text NASA ADS  
• Noh, W. F. 1987, J. Comput. Phys., 72, 78 In the text  
• Osorio, M., Lizano, S., & D'Alessio, P. 1999, ApJ, 525, 808 In the text NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1997, ApJ, 482, 712 In the text NASA ADS  
• Penston, M. V. 1969, MNRAS, 144, 425 In the text NASA ADS  
• Reid, M. A., Pudritz, R. E., & Wadsley, J. 2002, ApJ, in press [arXiv:astro-ph/0201261] v1 16 Jan. 2002 (RPW) In the text  
• Richer, J. S., Shepherd, D. S., Cabrit, S., Bachiller, R., & Churchwell, E. 2000, in Protostars and Planets IV, ed. V. G. Mannings, A. P. Boss, & S. S. Russell (Tucson: Univ. Arizona Press), 867 In the text  
• Shepherd, D. S., Yu, K. C., Bally, J., & Testi, L. 2000, ApJ, 535, 833 In the text NASA ADS  
• Shu, F. H. 1977, ApJ, 214, 488 In the text NASA ADS  
• Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 In the text NASA ADS  
• Sigalotti, L. di G., & Klapp, J. 2001, IJMP D, 10(2), 115 In the text  
• Sigalotti, L. di G., & Mendoza-Briceño, C. A. 2002, A&A, submitted In the text  
• Sod, G. A. 1978, J. Comput. Phys., 27, 1 In the text  
• van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J. II, & Blake, G. A. 2000, ApJ, 537, 283 In the text NASA ADS  
• Velusamy, T., & Langer, W. D. 1998, Nature, 392, 685 In the text NASA ADS  
• Ward-Thompson, D., Scott, P. F., Hills, R. E., & André, P. 1994, MNRAS, 268, 276 In the text NASA ADS  
• Ward-Thompson, D., Motte, F., & André, P. 1999, MNRAS, 305, 143 In the text NASA ADS  
• Whitworth, A. P., & Ward-Thompson, D. 2001, ApJ, 547, 317 In the text NASA ADS  
• Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846 In the text NASA ADS