© ESO, 2011

1. Introduction

The theory of star formation in the original and holistic sense, which treats the forming protostar as a real and structured object, remains poorly understood. Apart from early spatially one-dimensional calculations (beginning with Larson 1969) that have been constantly refined (see, e.g., Wuchterl & Tscharnuter 2003), until now no multi-dimensional numerical study has been able to follow the evolution of an unstable rotating molecular-cloud core to the end of the main accretion phase of the protostar after a few free-fall timescales for the original cloud.

In the past few years, several three-dimensional studies, either grid-based (e.g., Saigo et al. 2008) or developed using SPH methods (e.g., Walch et al. 2009), have been able to investigate the fragmentation of the first core (consisting of molecular hydrogen) and the transport of angular momentum. Some of these investigations have at least realized to follow the second collapse (forming the protostellar core), but either did not reach a quasi-static state for the stellar core, cf. Stamatellos et al. (2007), Saigo et al. (2008), or could follow the post-collapse evolution for only a short period of time (about 60 d, cf. Bate 1998). Two-dimensional studies were also able to reach this point with a much higher spatial resolution at the center and a more realistic equation of state (Tscharnuter 1987; Tscharnuter et al. 2009). Going beyond radiation hydrodynamics, Tomida et al. (2010) investigated the additional influence of magnetic fields on the structure and evolution of the first core, but did not cover the second collapse.

Three-dimensional calculations covering the first 50 yr of the second core formation were carried out by Bate (2010) to illustrate the impact of the protostar on the first core. We comment on these results and compare them with ours, considering in particular the rather coarse spatial resolution of Bate’s simulations.

In this paper, we present the detailed results of grid-based radiation hydrodynamical calculations of the collapse of a rotating supercritical Bonnor-Ebert sphere of 1   M_⊙, assuming rotational symmetry. In a parameter study, we investigate the influence of the β-viscosity (controlling the turbulent viscosity and hence the angular momentum transport) on the collapse evolution. The calculations cover the first and second collapse and follow the evolution of the newborn protostar and its surroundings for up to 240 yr, which corresponds to several thousand dynamical timescales of the central object.

2. Methods

2.1. Basic equations

The equations of mixed-frame flux-limited diffusion radiation hydrodynamics as derived by Krumholz et al. (2007) are used to model the collapse flow. The system of conservation laws then reads where ρ, ρv, e, and E are the densities of mass, momentum, internal energy, and radiation energy, respectively, so the Eqs. (1)–(4) describe the balance of these four quantities. The gravitational potential Φ in Eq. (2) is given by the solution of Poisson’s equation (5)Comments on the other variables appearing in Eqs. (1)–(4) are given in this and the next two sections Sects. 2.2 and 2.3. We use the well-known flux-limiter (6)derived by Levermore & Pomraning (1981) to model the transition from optically thick regions (s → 0, λ → 1/3) to thin ones (s → ∞, λ → 1/s). Furthermore, there is a factor (7)appearing in connection with the radiation energy density and of order unity (more precisely, ranging from 1 to 4/3 from the thin to thick limit). For more detailed information on the radiative terms in our equations, we refer to Krumholz et al. (2007).

2.2. Angular momentum transport and viscosity

Transport of angular momentum is one of the key problems when modeling protostellar collapse. There are essentially three mechanisms known to accomplish this transport:

• 1.

  Gravitational torques, e.g., exerted by spiral or bar-likestructures in the protostellar disk.

• 2.

  Turbulent viscosity based on hydrodynamically driven turbulence of the flow.

• 3.

  Turbulent viscosity caused by the magneto-rotational instability in magnetized disks.

Needless to say, all these transport processes are inherently three-dimensional. Hence, assuming rotational symmetry (as done in this work, cf. Sect. 2.5), it seems a bit “exaggerated” to perform an elaborated modeling of the respective mechanisms. That is why we adopted the β-viscosity prescription for the turbulent viscosity (Duschl et al. 2000) (8)where r is the radial coordinate, ω the local angular velocity, and β = 10^{−4   ···−2} the inverse of the critical Reynolds number (with β = 10^-3 being the “canonical” value). The β-viscosity is a mechanical analog and an extension of the α-viscosity (Shakura & Sunyaev 1973) modeled according to the thermal properties of the disk.

In addition we use the artificial tensor viscosity introduced by Tscharnuter & Winkler (1979) to broaden the emerging shock fronts. The respective kinematic viscosity reads (9)where l is on the order of the local width of the grid. The effect of ν_a is restricted to the (r,θ)-plane, since owing to the axial symmetry assumption no shock fronts can appear in the azimuthal direction. If we set β = 0, too, local conservation of angular momentum is retained. Using Eqs. (8) and (9), we can write down the viscous stress tensor (10)(with the unit tensor I) used in Eqs. (2) and (3).

2.3. Equation of state and opacity

The pressure p, temperature T, and opacity κ (we use the Rosseland mean opacity) in Eqs. (2)–(4) are functions of e and ρ alone for a given chemical composition, which is chosen to be solar (Pop. I) with X = 0.708, Y = 0.275, and Z = 0.017. The constitutive equations account for the inner degrees of freedom of H₂, as well as the H₂ dissociation, the ionisation of H and He, and the influence of H⁻. The opacity contributions of dust particles are included as outlined in Tscharnuter & Gail (2007). Figure 1 shows κ(ρ,T) in the density and temperature regime relevant for the process of star formation. The opacity gap due to dust evaporation is clearly visible, as well as the “kappa mountain” primarily produced by the H⁻-ion and H-ionization.

2.4. Numerical techniques

We use a finite-volume scheme on a staggered grid in spherical coordinates (r,θ,ϕ) to discretize the system of Eqs. (1)–(5). In the (r,θ)-plane, we construct a self-similar, locally uniform, and static grid by choosing a logarithmic grid spacing in the radial direction, resulting in a gradually increasing resolution towards the center. The general discretization follows the “ZEUS-like” approach given, e.g., in Hayes et al. (2006). A two-step operator-splitting technique is applied to solve the equations. The first source step is divided further into an explicit evaluation of the momentum sources (right-hand side of Eq. (2)), followed by an implicit solution of the energy equations (Eqs. (3) and (4)). The time-implicit solution is necessary because in optically thin regions the “radiative” timestep becomes easily much shorter than the respective hydrodynamical one. In the second transport step, the intermediate values gained from the source step are used to evaluate the advection terms explicitly. A second-order van-Leer interpolation scheme is employed to carry out the advection.

The linear systems of equations, arising from the discretization of Poisson’s Eq. (5) and the implicit form of the discretized energy equations, are solved by using a direct solver (in this case the software library PARDISO¹). Direct solvers have the advantage that they give the exact solution of the linear system (in contrast to all other solvers, which always yield approximate results), thus giving rise to a more accurate and robust overall numerical scheme.

2.5. Symmetry assumptions and boundary conditions

A major restriction of this investigation is that we assume rotational symmetry, which is necessary because of today’s limitations of computer power. There are of course configurations of rotating accretion disks where rotational symmetry is broken, leading to spiral structures, bars, or even binary stars etc., but we do not address these phenomena here.

Beyond that we assume equatorial symmetry (saving half of the computation time). The stability of this symmetry is by far not as critical as the rotational one, although there is no proof of it for complex systems such as this one. Nevertheless, comparative simulations do not show any tendency of a deviation from equatorial symmetry.

All “inner” boundary conditions are determined by these symmetry assumptions, whereas the outer boundary of our computational domain is defined using additional information. We assume the outermost sphere to be fixed in space without any convective flux (i.e., no mass transfer). Radiation is streams into the surrounding medium, which is in thermal equilibrium and remains at the initial cloud temperature. The respective prescription for the radiation field at the boundary can be found in Yorke (1980).

In the last simulation of this paper (with β = 10^-2), we additionally investigate the long-term evolution of the system after the first accretion phase of the protostar (cf. Sect. 3.9). Since the timescale of the central object is by far the shortest, we can only advance the system over a reasonable time span (>10 kyr) if we exclude the central parts for the calculation. In this case, we therefore need an inner boundary condition. Since this condition is especially tailored to the system, we describe the details in the discussion of the results in Sect. 3.9.

2.6. Initial conditions

The outcome of protostellar collapse calculations is highly sensitive to the initial parameters. Since we only begin to grasp the detailed processes of star formation theoretically, a rather simple initial configuration should be used. However, the oversimplified approach of uniform density initial conditions (widely employed for several decades) should be questioned, because this object, if gravitationally unstable, can never exist in equilibrium. In contrast, we know from observations that there are real objects (such as Barnard 68) for which there is a striking agreement with uniformly rotating Bonnor-Ebert spheres (BES) (see, e.g., Alves et al. 2001 or Lada et al. 2003).

We therefore use a slightly supercritical BES² with a total mass of 1   M_⊙ and a temperature of 10   K as the initial object for all calculations. Furthermore, a uniform rotation with angular velocity ω is assumed, given in terms of the centrifugal barrier of the mass elements farthest away from the rotation axis as (11)where r_max is the radius of the outermost sphere and M the total cloud mass. The angular velocity itself is not a very convenient and descriptive parameter as it scales nonlinearly (quadratic) with the rotational energy. We choose R_c = 100   AU throughout all simulations. This corresponds to a revolution period of 2π/ω ≈ 7   Myr and a ratio of rotational to potential energy of 5.4 × 10^-3.

Although collapse calculations with uniformly rotating, marginally supercritical BES using SPH methods (cf. Walch et al. 2009) have already been carried out, there is still no “initial standard configuration” for this kind of numerical study. There is no doubt that an agreement on this issue, as well as reasonable parameter studies, would lead to a clearer understanding of the great variety of phenomena related to the formation of stars and protostellar disks.

3. Results and discussion

3.1. Test calculations

As we developed the numerical code from scratch, we solved numerous test problems to validate the general numerical approach. Besides several one- and two-dimensional advection tests, the pressure-free collapse and Sod’s shock tube, we also tested the accuracy of the solution of Poisson’s equation and problems including the diffusion of radiation. Since the focus of this work is on modeling real world physics and not on code development (the applied numerical methods are well-known), we will not discuss the test results here aside from mentioning that all of them came out as expected. One aspect which goes certainly beyond a simple test problem is the local conservation of angular momentum in the absence of viscosity (cf. Sect. 2.2). We address this point in the simulation for β = 0 (see Sect. 3.5).

3.2. General remarks on the simulations

In a parameter study we have investigated the influence of the β-viscosity by varying β between 0, 10^-4, 10^-3 and 10^-2. The spatial resolution was chosen to be 320 grid points in the radial and 35 grid points in the polar direction (between the rotation axis and the equator). The innermost (central) resolution is about 0.1 R_⊙, whereas the total radial extent of the computational domain is r_max ≈ 8700   AU. The ratio of the largest to the smallest grid size is not less than 2 × 10⁷! We performed comparative runs with higher spatial resolutions by taking 2 and 4 times as many grid points as in the runs shown here, resulting in no noticeable quantitative changes. Since we could not perform these higher resolution simulation for the whole first accretion period of the stellar core (because they are much more time consuming), they are not shown here.

The accuracy during the calculations was monitored by checking the total energy balance, which is given by (12)where the last term ϵ_out gives the accumulated radiative energy streaming out of the computational domain (F being the radiative flux). Ideally, ϵ_tot should be constant in time. As the only available energy source is the potential energy ϵ_pot, we can define a relative total energy error through (13)with the initial values ...(0). Throughout all simulations Δε_rel stayed below 5%, which is regarded as a reasonable value.

3.3. Classification of early evolutionary phases

As we know from 1D spherical symmetric models, we can divide the collapse evolution into three distinct phases. Phase 0 includes the quasi isothermal contraction of the Bonnor-Ebert-sphere on the free-fall time.

The beginning of phase 1 is marked by the so-called first collapse (on a dynamical timescale) and the formation of an optically thick central object, or more precise, if the relation (14)holds, i.e., if the total radial optical depths from all angular directions are at least 2/3. Relation (14) is the natural three dimensional generalization of the stellar zero age concept in Wuchterl & Tscharnuter (2003). Following this definition, we reset the time to zero when (14) holds, therefore phase 1 starts at t = 0.

If the first hydrostatic core heated up significantly (because of accretion and contraction) so that by means of the dissociation of the hydrogen molecules an instability sets in, the second collapse to protostellar scales (i.e., to the second core) can proceed, defining the start of phase 2. For practical reasons, we consider this point to have been reached when the maximum temperature in the first core reaches 2000 K.

3.4. Phase 0

The contraction of the Bonnor-Ebert sphere until the first collapse develops practically independent of the value of β, because the rotational velocities and centrifugal forces are still too small to have a noticeable effect. The whole phase 0 lasts about 160 kyr (1.2 free-fall times) which is short compared to the rotation period of 7 Myr, hence there are only slight deviations from spherical symmetry.

Figure 2 shows the evolution of several quantities during Phase 0. The central density increases by about nearly four orders of magnitude, while the whole cloud remains practically isothermal (at 10 K). Only over the last 10 kyr of phase 0 has the temperature risen from 10 to 18 K. This is a result of the small optical depth staying below 0.01 for all but the last 10 kyr of phase 0. When Eq. (14) holds for the first time, i.e., the starting point of phase 1, we reset the clock to zero as mentioned before. All ages given below correspond to this instant of time.

3.5. Simulation for β = 0

Without any transport of angular momentum, the first core evolves into a more and more flattened quasi-hydrostatic rotating disk surrounded by an accretion shock. As a result of the incoming material with increasingly larger specific angular momentum (and corresponding greater centrifugal barrier), an off-center density maximum evolves after about 5 kyr. Shortly after 14 kyr, the disk becomes dynamically unstable in a region approximately 5–6 AU from the center producing a ring-like structure. This finding is in perfect agreement with the Toomre stability analysis (cf. Toomre 1964). After 17 kyr, a second ring emerges in the central parts (<4 AU) of the disk.

The overall structure of the object at the end of the simulation after 28 kyr is depicted in Fig. 3 in the form of density and temperature distributions. The two rings are well separated by a clear gap between them. The disk has reached a radius of approximately 27 AU and already contains more than 43% of the total mass of the cloud core. The outer parts (>10 AU) are nearly isothermal owing to the small optical depths. Because of the efficient radiation of the large disk surface, the temperature stays below 1400 K, which is typical of dust evaporation. Conditions for the dissociation of hydrogen molecules are certainly not reached.

There is no doubt that the formation of rings in the accretion disk would (in general 3D) soon be followed by a symmetry breaking, i.e., the emergence of non-axisymmetric structures, e.g., such as spirals, bars, or the onset of binary formation. However, with our parameter study of β (in Sects. 3.6–3.8) we obtain an idea of how much turbulent viscosity (see Sect. 2.2) is at least needed to prevent the formation of these structures.

Figure 4 shows the mass spectrum of the specific angular momentum (l = ωr²sin²θ) in the form of the mass fraction M(l) with a maximal specific angular momentum of l. Since M(l) is theoretically constant in time for a non-viscous axially symmetric flow, we can test for the local conservation of angular momentum by monitoring possible deviations from the initial spectrum (cf. Norman et al. 1980). Fortunately, less than 1% of the total mass exhibits a significant deviation at any time. Until about 15 kyr, only a minor redistribution of angular momentum is visible in a region l ≈ 10^-6−10^-4l_max. After 25 kyr, we find a clear drop in the spectrum for l ≲ 10^-3l_max, which is produced during the second ring formation at about 17 kyr. Owing to numerical diffusion too much angular momentum has been advected away from the rotation axis in this process.

3.6. Simulation for β = 10^-4

The quite small value of β = 10^-4 produces a qualitatively different evolution than the β = 0 case. The viscous timescale τ_visc = 1/(βω) for inner disk regions (<1 AU) is comparable to the evolution timescale of the disk (<10 kyr). We therefore see a more centrally concentrated mass distribution owing to outwardly directed angular momentum transport. The inner regions have higher temperatures because of viscous heating. After about 7.7 kyr, the conditions for the dissociation of hydrogen are reached, again in a ring-like density concentration at approximately 0.35 AU.

The disk structure at this moment is shown in Fig. 5. The visible indentation in the pole region of the density distribution is a result of the thermal restructuring caused by the dust opacity gap. We discuss this phenomenon in more detail in Sects. 3.7 and 3.8. Only after some 10 yr (at about 7.73 kyr) does phase 2 begin and the ring collapse on a dynamical timescale, triggered by the hydrogen dissociation of the outer ring parts. By means of this event, the simulation comes to a “natural” end, since the spatial resolution requirements needed to follow the ring collapse are not satisfied outside the center of the grid. The remarks given at the end of Sect. 3.5 concerning the possibility of additional non-axisymmetric evolution after ring formation again apply in this calculation.

3.7. Simulation for β = 10^-3

Since the two remaining simulations lead to the formation of an object on protostellar scales in the center of the cloud, which our numerical scheme can indeed resolve, we discuss them in more detail.

3.7.1. First core and dust evaporation

The β = 10^-3 model results in a more compact and warmer inner part of the disk than the β = 10^-4 case. The maximum density and temperature remain at the very center. Thus, dust evaporation giving rise to the formation of the opacity gap already occurs after 2.4 kyr. The global density and temperature distributions of the disk at this moment are shown in Fig. 6. The process of dust evaporation has severe consequences for the general structure and the further evolution of the central disk regions. Figure 7 illustrates the expansion of the opacity gap in terms of the absorption coefficient κρ (the inverse mean free path). The resulting quasi-isothermal region in the dust-free part leads to a “blow-up” of the disk stratification starting at about 2520 yr. As the temperature gradient at the border of the opacity gap becomes too large, these layers become dynamically unstable and convection sets in at 2640 yr (Fig. 7 upper right panel). The blow-up of the polar disk regions is now clearly visible. After about 2670 yr (Fig. 7 lower right panel), an optically thin (dust-free) channel stabilized itself around the rotation axis, allowing the heat to be efficiently radiated away in this direction. Like piercing a hole into a balloon, the blown-up disk layers then begin to subside again.

3.7.2. Second collapse and protostar formation

The flow field at 2680 yr (only 10 yr after the last picture in Fig. 7) is shown in Fig. 8 and marks the beginning of phase 2. The accretion onto the blown-up first core is ongoing with highly supersonic velocities (>Mach 5). The massive convection resulting from the general restructuring caused by dust evaporation (see previous section) is impressively demonstrated in the form of the strong vortex with speeds up to Mach 0.8, centered about 2 AU above the equatorial plane. The subsonic motions in the (optically thick) disk are also influenced by the flow near the opacity gap and display quite erratic behavior in this phase. Amid this very non-stationary situation, the second collapse (triggered by the hydrogen dissociation) sets in, visible through the radial inflow near the rotation axis at a height of 1 to 2 AU. The subsequent evolution of the collapse flow, documented in Fig. 9 and 10, conflicts with our current understanding of the second collapse to protostellar scales, which is assumed to occur on the local dynamical timescale (in this case within 1 yr). We instead see a kind of “crawling” accretion shock, moving with approximately constant velocity during a period of about 20 yr towards the center. The accretion velocities continuously increase with time and the supersonic infall region becomes more and more cone-like shaped (see Fig. 10). The whole process takes place on a thermal timescale, although the exact estimate of this timescale is difficult, because the complex geometrical and thermal structures do not simply allow us to define an object that we can assign both a thermal energy and luminosity. However, we can exclude the viscous timescale (1/(βω) ≈ 200   yr), which is one order of magnitude longer than the collapse phase. The main reason why the dynamical collapse cannot proceed as in the well-known spherical case, is the strongly disturbed force distribution due to centrifugal forces in the rotating system.

About 10 yr before the onset of the second collapse, the outermost parts of the first core allow the absolute ratio of total rotational to gravitational energy of the first core to exceed the critical value of 0.274, where it is likely that the object becomes prone to the growth of non-axisymmetric structures. Since the rotation period of the outer first core (which should be related to the typical timescale for the growth of non-axisymmetric features) is about 125 yr, there is apparently not enough time to transport angular momentum out of the central regions of the first core and influence the second collapse. Therefore, it seems realistic that the second collapse can occur on the thermal timescale rather than the dynamical one.

Figures 11 and 12 show the temperature and density structure of the innermost parts during the formation of the protostellar core, and in addition its subsequent evolution for 236 yr. The upper left panel of Fig. 12 shows the arrival of the cone-shaped accretion front (see previous paragraph) on stellar spatial scales after 2700 yr. Only 0.7 yr later does the front rebounce at a height just about 0.4   R_⊙ above the center (Figs. 11 and 12 lower left panel) because of the newly created pressure gradients in the now fully dissociated gas. A small flare is visible in the temperature distribution generated by the complex velocity field during the rebounce. About 50 yr later (Figs. 11 and 12 upper right panel), the high density region (ρ > 0.1   kg   m^-3) has become almost spherical with a radius of approximately 1   R_⊙. The main accretion now occurs in a small zone around the rotation axis. However, a small mass contribution also reaches the central parts through the protostellar disk in a radial inward flow along the equator. This transport occurs on the viscous timescale and can be inferred indirectly in the form of the cooler equator region in the temperature distribution. Apart from the compact central object and the equator, the whole environment around the newborn star is characterized by a highly variable flow field in this phase. At the end of the simulation, 236 yr after the second collapse (Figs. 11 and 12 lower right panel), the main accretion virtually stops because of missing material (the density above the star and the disk has decreased by several orders of magnitude, see next section). A comparatively “quiescent” flow field develops and both the protostar (with a radius of 2–3 R_⊙) and the surrounding disk form a quasi-static object.

3.7.3. The protostar’s impact on the first core

The main accretion onto the second (protostellar) core consists of an intricate (non-stationary) matter flow confined to a narrow region around the rotation axis. This process has a substantial influence on the structure of the first core (with spatial scales of about 10 AU). Figure 13 shows the disintegration of the first core’s global structure in terms of the absorption coefficient for phase 2 of the evolution. Although most of the accretion energy is radiated away along the rotation axis in the optically thin and dust-free area (clearly visible in Fig. 13), the overall heating of the core regions adjacent to the “tunnel” starts to erode the whole object through a process with positive feedback: the diffusing radiation heats the gas and allows it to expand, hence the density is reduced (and so is the absorption coefficient κρ) permitting the radiation to diffuse into the gas even more easily. An interesting question is how this result would change if we applied a more realistic model of the radiation field instead of using the simplified diffusion approximation.

To which impressively large scales the expansion of the first core leads in only about 240 yr of the protostar’s life can be seen in Fig. 14. A nearly spherical outflow, becoming supersonic already after about 50 yr of the protostar formation, disperses the material roughly homogeneously. After 2936 yr (the end of the simulation), the outward travelling shock wave has reached a radius of more than 80 AU. A very diffuse, low-mass disk with a radius of only 5 AU survived around the central object. It is now evident why the density above the protostar and its disk strongly decreased and the accretion rate practically vanished: the central protostar has simply blown away all its surrounding gas masses. The question of whether or not the ejected matter will return to the center can be generally confirmed, since the total binding energy of the system is, of course, still negative. How the process of re-accretion appears in some detail will be investigated in the last simulation with β = 10^-2 (cf. Sect. 3.9).

The extreme impact of the newborn protostar on the first core was also found by Bate (2010), though the spatial scales and velocity field of the outflow differ significantly for this case. The process occurs more slowly (by a factor of about 4) in our calculation and shows no indications of a bipolar flow field, but rather a spherical one. Since the initial conditions are not identical and we assume an explicit prescription of turbulent viscosity, a more robust comparison of these results is not possible at this point.

3.7.4. Evolution of some global and midplane variables

Figure 15 shows global quantities during their evolution in phase 2. The total luminosity (Fig. 15, upper left panel) has proven to be a quite sensitive indicator of any kind of change in the cloud’s evolution. The 20-yr-long second collapse leading to the formation of the protostellar core (cf. Sect. 3.7.2) is connected with a monotonic decrease in the luminosity between 2.68 and 2.7 kyr. This is because the compressional heat goes directly into the dissociation of hydrogen during this phase. The subsequent rise and heavy fluctuation of the luminosity is due to the first stages of the main accretion onto the second core with a very complex and time varying flow field (and accretion rate) in the narrow “tunnel” along the rotation axis. Besides some more or less periodic phases with periods of about 20 yr, the magnified inlay shows a short time period of about 1 yr related to a cyclic instability at the innermost accretion front right above the protostar. Values of up to 10   L_⊙ were reached during phase 2. At the end of the simulation, the luminosity decreases rapidly because it is directly connected to the accretion rate.

Two masses relative to the total mass are displayed in the lower left panel of Fig. 15. On the one hand the relative subsonic mass M_sub, which consists of all subsonic mass elements inside the global inflow from the cloud core³. Otherwise, the relative “stellar” mass M_st, defined as the sum of all mass elements with densities  >0.1   kg   m^-3. Until the first 50 yr after the formation of the protostellar core, M_sub is constantly increasing, reaching values of nearly 10%. As the outflow (initiated by the protostar) becomes supersonic, M_sub commences to decrease. Finally, at 2936 yr M_sub   (= 1.77%) is only somewhat greater than M_st   (= 1.54%), i.e., almost all of the remaining subsonic matter is inside the star. The first core has lost 82% of its mass in this process.

Compared to the value of M_st (defined in the same way as here) given in Bate (2010), we find a somewhat different evolution. The increase in M_st during the earliest phase of the formation is only about 0.2%, visible in the magnified inlay in Fig. 15 (also note the small “overshoot” there – an indication of the rebounce process), whereas Bate found a stronger first increase to about 0.7%. However, after that, the mean accretion rate is much higher in our calculation (about 6 × 10^-5   M_⊙   yr^-1) than Bate’s findings of about 10^-5   M_⊙   yr^-1. A more detailed comparison with the protostellar masses in Bate (2010) is difficult, since the runs of M_st given there seem to depend significantly on the spatial resolution of the calculation.

The upper right panel in Fig. 15 shows the evolution of the maximum density to be practically identical to the central value. The density rise over 6 orders of magnitude (between 2.68 and 2.7 kyr) proceeds far more slowly than the dynamical timescale (<1 yr), confirming the assumption of a secular instability occurring on the thermal or rather viscous timescale. The central temperature (Fig. 15 lower right panel) does not increase strongly until the main accretion starts. That after second core formation the central temperature falls more and more below the maximum value is easily explained. Since the main heating comes from the accretion front in the outermost layers of the star and the matter in this temperature and density range has the highest opacity (κ) possible at all (cf. Fig. 1), the central temperature only adjusts itself on the very long (>Myr) thermal timescale of the star. That is why the central temperature stays below the maximum value for the longest part of the formation time of a star.

The evolution of the midplane density and temperature distribution is plotted in Figs. 16 and 17. After the ejection of most of the material at the end of the simulation (at 2936 yr), the density is noticeably lower than before in the region of the former first core. The outwardly traveling shock is clearly visible at about 80 AU. In the temperature distribution, we see the above-mentioned inversion inside the protostar (i.e., a cooler center) and the pronounced “dust knee” at about 1400 K. Before the destruction, the temperature scales approximately with R^-1 in the dusty region of the first core. At the end, this temperature structure is completely lost.

3.8. Simulation for β = 10^-2

3.8.1. Dust evaporation

The phase preceding the dust evaporation is qualitatively comparable to the β = 10^-3 case, though the higher turbulent viscosity produces a more spherical central region (through effective angular momentum transport) and a stronger viscous heating.

Therefore the opacity gap occurs as early on as about 1300 yr. The following dynamical restructuring is shown in Figs. 18 and 19. There are some essential differences from the β = 10^-3 calculation. Since the temperature increase is generally more rapid in this case, the blow-up of the disk stratification and the onset of convection at the border of the opacity gap occurs far more rapidly. However, in this case there is no way for the optically thin “eddies” to stabilize a dust-free channel (cf. Sect. 3.7.1) that would relax the blow-up. The reason is the different disk stratification, which is not as thin (because the process occurs much earlier in time and the disk scale-height decreases with time) and more spherically shaped (compare with Fig. 7).

Since the heat cannot be radiated away as efficiently as in the same phase of the calculation for β = 10^-3, the blow-up of the inner disk layers continues. The character of the outflow becomes more and more that of an explosion, reaching Mach numbers of 3 (see Fig. 19). This demonstrates the massive impact that the dust evaporation can have on the structure of the inner disk parts.

3.8.2. Second collapse and protostar formation

The second collapse proceeds on the dynamical timescale (of shorter than 1 yr), a familiar finding for spherical symmetric calculations. This is plausible since much more angular momentum has been transported out of the central parts than in the β = 10^-3 case. The radial inflow during the collapse now occurs from all directions, though significantly slower along the equator.

The formation of the central protostar is depicted in Figs. 20 and 21 through temperature and density distributions. The first 5 yr after the collapse are characterized by a highly variable flow field visible through the numerous hot flares that give the impression of a “campfire” (Fig. 20 left side). By means of the decreasing accretion rate, a stationary flow field develops after less than 10 yr; a hot polar cap and a hot disk surface build up, separated by a small slit where material enters the protostar (Fig. 20 upper right panel). The flattened, stellar-scale accretion disk forming instantly after the onset of the second collapse at about 1409 yr reforms an approximately spherical central object plus a diffuse disk at the equator after only some 10 yr through angular momentum transport (cf. Fig. 21). After about 70 yr, the accretion has virtually stopped again for the same reason as in the β = 10^-3 case (see next section).

3.8.3. The protostar’s impact on the first core

The accretion onto the protostar has an even stronger influence on the first core than in the previous calculation, since the conversion from kinetic to thermal energy occurs on a much shorter timescale (about 10 yr), which results in a more intense and short-time accretion luminosity (see also next section).

The erosion of the first core through the same positive feedback process discussed in Sect. 3.7.3 is shown in Fig. 22 in the form of the absorption coefficient. At 1409.5 yr, the structure of the first core is still visible, including the “ice line” (a slight increase in κρ, reaching the equator at R ≈ 9   AU). The outflow remains governed by the restructuring through the dust evaporation. At 1413 yr, the luminosity of the protostar reaches a maximum, resulting in the farthest extended dust-free zone that allows the radiation to penetrate the disk even more easily. The visible wave-like structure (in Fig. 22, lower left panel⁴) inside the opacity gap develops because the dust-free region forms a kind of cavity for the thermal waves. The border of the opacity gap reflects a part of the outgoing thermal waves because of the sudden increase in the opacity. Therefore, standing thermal waves build up, which are hardly visible in the temperature distributions (being of too small amplitude) but are clearly pronounced in the absorption coefficient because of its strong dependence on temperature at the slope of the opacity gap (cf. Fig. 1). Only 30 yr after the second collapse, the structure of the former first core is completely destroyed by the impact of the protostar (Fig. 22, right side).

The global outflow produced by the heating of the central object is shown in Fig. 23 in terms of density distributions. Starting with the onset of the second collapse at about 1410 yr (Fig. 23 upper left panel), where the first core is only perturbed in the central parts, the subsequent frames document the blow-out of the first core to spatial scales of 80 AU in the pole direction and more than 40 AU along the midplane in about 70 yr. That the outflow is not exactly spherical (as in the β = 10^-3 case) may be because of the preceding super sonic blow-out produced by the dust evaporation. This gives an additional radial outward momentum in the polar region that would lead to the obtained result.

Interestingly, the velocity scales in this case are comparable to those of Bate (2010), i.e., the front of the blown-out first core reaches on average about the same radial distances in the same time for both simulations. But again, we see no indication of a bipolar outflow. As said before, because of the different initial conditions and the turbulent viscosity, we cannot compare these results in any greater depth.

3.8.4. Evolution of some global and midplane variables

As for the previous simulation, we summarize the evolution of some global quantities (Fig. 24). The total luminosity run (Fig. 24, upper left panel) confirms the picture of a short and intense accretion phase lasting only about 5 yr. A maximum value of 50   L_⊙ is reached at 1413 yr. Between 1420 and 1430 yr a constant accretion luminosity (at  ≈3   L_⊙) is visible, followed by a clear drop owing to the strongly decreasing accretion rate. Compared to the β = 10^-3 calculation (cf. Fig. 15), the luminosity evolution is much more smooth here, a consequence of the quasi-stationary accretion flow that develops directly after a short (<5 yr) and highly variable initial phase.

The relative subsonic (M_sub) and stellar (M_st) mass, introduced in Sect. 3.7.4, are depicted in the lower left panel of Fig. 24. The first small drop in M_sub at 1409 yr is related to the dynamical collapse at the very beginning of the protostar’s formation, where M_st gains about 0.7% within only 3 yr after the second collapse. At about 1412 yr, the blow-out of the first core reaches supersonic velocities and M_sub begins to decrease. At 1430 yr, virtually all material of the former first core is in supersonic radially outward motion, almost all of the remaining subsonic (= 1%) mass being contained in the protostar. This corresponds to a mass loss of about 84% for the first core. Compared to the M_st given in Bate (2010), we again find acceptable agreement at least for the highest spatial resolution given there (3 × 10⁶ particles). Unfortunately, we do not know whether Bate’s results at the highest resolution are already near a converging point, i.e., whether an even higher resolution would not yield a significantly different outcome.

The maximum (and identical central) density evolution is shown in the upper right panel of Fig. 24. The far more abrupt increase in density over six orders of magnitude in less than one year is clearly visible (compared to more than 20 yr in the simulation for β = 10^-3), again pointing out the dynamical character of the second collapse. The same is true for the maximum and central temperature evolution given in the lower right panel of Fig. 24. The reason that the central temperature falls more and more below the maximum value has already been discussed in Sect. 3.7.4.

The evolution of the midplane density and temperature distribution is depicted in Figs. 25 and 26. The density in the region of the former first core after the blow-out at the end of the calculation (at 1482 yr) is even lower than in the previous simulation (cf. Fig. 16), and this in only one fourth of the time. The temperature distribution of the protostar again shows the off-center maximum value and the thermal restructuring after the second collapse.

3.9. Subsequent evolution of the blown-out core for β = 10^-2

At the end of the β = 10^-2 calculation, after only about 70 yr of protostar formation, the key question is what the further evolution of the ejected matter from the former first core will be.

3.9.1. Technical realization of an inner boundary

The protostar is now in a quasi-static state, since the accretion rate has decreased almost to zero. This gives us the opportunity to cut out the inner region containing the star. At a radius of 0.7 AU, we define a polar-angle-dependent mixed boundary condition consisting of a supersonic radial inflow condition beginning at the rotation axis and a given pressure condition near the equator (coming from the quasi-static structure of a diffuse, almost massless accretion disk re established after the first core blow-out). The mass distribution inside the central “hole” is assumed to be spherical. The radiative flux at the boundary (even though small now) is calculated from energy balance considerations of the masses accreted by the hole (using a mean accretion radius of the known protostar scales). A deeper discussion of this artificial boundary would lengthen the section unreasonably, hence for more details we refer to Schönke (2010)⁵.

3.9.2. Results

Figures 27 and 28 show the process of matter ejection and subsequent re-accretion during a time span of 56 kyr. After 2 kyr, the blown-out material has been compacted to form a shell of material that is still expanding supersonically (Fig. 27 upper left panel, compare with the previous state in Fig. 23 lower right panel at 1482 yr). At 4 kyr, the shell has reached the turning point with a maximum radial extent of about 500 AU, the density between the shell and the innermost parts has decreased substantially to values comparable with the initial BES (Fig. 27 lower left panel). After 5–6 kyr, the inner shell parts begin to fall back to the center. Interestingly, a large fraction of the material first concentrates around the rotation axis and then falls back in a thin tube directly above the center (Fig. 27 lower right panel). At 8 kyr, we see the subsequent bidirectional radial inflow along the rotation axis and from the equatorial plane (Fig. 28 upper left panel). What follows is the formation of a new diffuse core with scales of about 100 AU that re-arranges itself with time to a widely extended accretion disk spanning a radius of more than 400 AU after 56 kyr (Fig. 28 right side).

At 13 kyr, all of the blown-out material from the former first core has been re-accreted, and after that, again only “fresh” material from the outer cloud parts is accreted. At the end of the calculation (after 56 kyr), 56% of the total mass is already inside the “hole” (i.e., the protostar) and 13% inside the new accretion disk, so more than 69% of the total mass has reached the central parts in only one third of the initial BES free-fall time since the beginning of phase 1.

Whether or not the re-accretion of the blown-out material occurs indeed in such a “quiet” way as found here, cannot be surely decided, because we do not model the true accretion process onto the protostar, which could possibly lead to another blow-out. The main intention of this additional long-term calculation was to determine the typical timescale of the ejection process (about several kyr) and investigate the way in which the material flows back to the center.

4. Conclusions

We have uncovered several phenomena related to the protostellar collapse that have never been previously discussed. The variation in the turbulent viscosity and therefore the effectivity of angular momentum transport obviously leads to significantly different evolutionary scenarios.

The calculations generally indicate that the phase of dust evaporation inside the first core is able to initiate a major dynamical restructuring. Depending on the state of the first core, the dust-free (low opacity) regions can either find an optically thin path to the core’s surface (providing an effective cooling

mechanism, cf. the β = 10^-4 and 10^-3 case) or the new thermal structure produces a rapid blow-up of the upper – core layers (cf. the β = 10^-2 case). However, the dust evaporation process has a decisive influence on the evolution of the first core and therefore on the way the second collapse proceeds. This is important and, to our knowledge, has not been pointed out so far in the literature.

A completely new finding is the possibility that the second collapse does not necessarily take place on the short dynamical timescale, but rather on the much longer thermal timescale (cf. the β = 10^-3 case). Which timescale is relevant depends on the geometry of the first core and the force distribution during the phase of hydrogen dissociation. The more spherical the collapsing region the more dominant the dynamical timescale, i.e., whether the angular momentum transport is effective enough to reduce the ratio of centrifugal to gravitational forces sufficiently, the well-known dynamical collapse proceeds (cf. the β = 10^-2 case).

It has been clearly illustrated that the newborn protostar has a huge impact on its surroundings (i.e., the first core). Driven by the stars’ thermal pressure⁶ during a short (on the order of 100 yr) accretion phase, the first core is easily blown out and destroyed, since it is centrifugally supported. The qualitative agreement with the results obtained by Bate (2010) shows that we have to dismiss the picture of a protostellar accretion disk that simply evolves uninfluenced by the star.

The overall process of star formation may be thought of as a highly non-stationary and cyclic one. The short and luminous accretion phases would then alternate with long (on the order of several kyr) blow-out phases characterized by extremely low luminosities. Future investigations will need to clarify whether this is really the case.

---

Acknowledgments

This work has been supported by the Deutsche Forschungsgemeinschaft within the Forschergruppe 759 “The Formation of Planets: The Critical First Growth Phase”. We thank the referee for the constructive and thoughtful comments.

References

All Figures