Online Material

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

Boundary conditions

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

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

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

       

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

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

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

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

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

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

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

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

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

  

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

----------

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

Appendix B: Numerical techniques and improvements

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

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

B.1 Spatial resolution

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

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

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

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

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

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

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

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

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

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

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

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

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

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

B.2 Advection

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

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

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

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

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

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

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

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

B.3 Hydrostatic parts and total energy

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

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

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

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

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

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

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

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

B.4 Discretization remarks

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

B.5 Constructing the initial conditions

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