Astronomy&
Astrophysics

A&A 526, A139 (2011)
DOI: 10.1051/0004-6361/201015734
c ESO 2011
Protostellar collapse of rotating cloud cores
Covering the complete first accretion period of the stellar core
J. Schönke and W. M. Tscharnuter
Zentrum für Astronomie (ZAH), Institut für Theoretische Astrophysik (ITA), Universität Heidelberg, Albert-Ueberle-Str. 2,69120 Heidelberg, Germany
e-mail: schoenke@ita.uni-heidelberg.deReceived 10 September 2010 / Accepted 19 November 2010
ABSTRACT

Aims. We investigate the influence of turbulent viscosity on the collapse of a rotating molecular cloud core with axial symmetry, inparticular, on the first and second collapse phase, as well as the evolution of the second (protostellar) core during its first accretion
period. By using extensive numerical calculations, we monitor the intricate interactions between the newly formed protostar and thesurrounding accretion disk (the first core) in which the star is embedded.
Methods. We use a grid-based radiation-hydrodynamics code with a spatial grid designed to meet the high resolution required to study
the second core. The radiative transfer is treated in the flux-limited diﬀusion approximation. A slightly supercritical Bonnor-Ebert
sphere of 1 M and uniform rotation according to a fixed centrifugal radius of 100 AU serves as the initial condition without exception.
In a parameter study, we vary the β-viscosity driving the angular momentum transport.
Results. Without viscosity (β = 0), a highly flattened accretion disk forms that fragments into several “cold” rings. For β = 10−4 , a
single “warm” ring forms that undergoes collapse due to hydrogen dissociation. For β = 10−3 , ring formation is suppressed completely.
The second collapse proceeds on the local thermal timescale, which is in contrast to the current view of a generally dynamical second
collapse. During the first accretion period of the second core, the first core heats up globally and, as a consequence, a nearly sphericaloutflow occurs, destroying the structure of the former accretion disk completely. Finally, for β = 10−2 , we see the classical dynamicalsecond collapse and a shorter but more rapid accretion phase. The impact on the surrounding accretion disk is even more pronounced.
We follow the resulting massive outflow up to several kyr after the second collapse, where the central parts (R < 0.7 AU) are now cut
out and replaced with an appropriate inner boundary condition. Matter is found to turn back to the center at a radius of 500 AU after
about 4 kyr and to reach the protostar again after approximately 7 kyr. The results suggest that the star formation process consists ofshort and rapid accretion phases (lasting on the order of 100 yr) between long and quiet outflow periods (lasting several kyr).
Key words. stars: formation – stars: winds, outflows – accretion, accretion disks – hydrodynamics – protoplanetary disks

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) thathave been constantly refined (see, e.g., Wuchterl & Tscharnuter
2003), until now no multi-dimensional numerical study has beenable to follow the evolution of an unstable rotating molecularcloud 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. Someof these investigations have at least realized to follow the second collapse (forming the protostellar core), but either did notreach 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 amuch 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 thesecond collapse.
Three-dimensional calculations covering the first 50 yr of thesecond core formation were carried out by Bate (2010) to illustrate the impact of the protostar on the first core. We commenton 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 andhence 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.

Article published by EDP Sciences

A139, page 1 of 14

A&A 526, A139 (2011)

2. Methods
2.1. Basic equations

The equations of mixed-frame flux-limited diﬀusion radiationhydrodynamics as derived by Krumholz et al. (2007) are used tomodel the collapse flow. The system of conservation laws thenreads∂ρ+ ∇ · (ρ u) = 0,∂t∂( ρ u)+ ∇ · (ρ u u) = −∇p + ∇ · Q − ρ∇Φ − λ∇E,∂t∂e
+ ∇ · (e u) = −p∇ · u + Q : ∇u
∂t


v2
+ cκρ E − aT 4 + λu · ∇E − f κρ E,
c


cλ
∂E
+ ∇ · ( f E u) = ∇ ·
∇E
∂t
κρ


v2
− cκρ E − aT 4 − λu · ∇E + f κρ E,
c

(1)
(2)

(3)
Fig. 1. The Rosseland mean opacity κ(ρ, T ).

(4)

where ρ, ρu, 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ΔΦ = 4πGρ.

(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


1
1
|∇E|
λ=
coth s −
,
s=
,
(6)
s
s
κρE
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


3 1
1
f = −
coth s −coth s(7)
2 2
sappearing in connection with the radiation energy density and oforder 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).

perform an elaborated modeling of the respective mechanisms.
That is why we adopted the β-viscosity prescription for the turbulent viscosity (Duschl et al. 2000)
νβ = β r2 ω,(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 amechanical analog and an extension of the α-viscosity (Shakura& Sunyaev 1973) modeled according to the thermal propertiesof 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
νa = −l2 min (∇ · u, 0) ,(9)where l is on the order of the local width of the grid. The eﬀect 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


1
(10)
Q = 2 ρ (νβ + νa ) ∇u − (∇ · u) I3(with the unit tensor I) used in Eqs. (2) and (3).

2.2. Angular momentum transport and viscosity

Transport of angular momentum is one of the key problems
when modeling protostellar collapse. There are essentially threemechanisms known to accomplish this transport:
1. Gravitational torques, e.g., exerted by spiral or bar-like structures 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 (asdone in this work, cf. Sect. 2.5), it seems a bit “exaggerated” toA139, page 2 of 14

2.3. Equation of state and opacity

The pressure p, temperature T , and opacity κ (we use theRosseland mean opacity) in Eqs. (2)–(4) are functions of e andρ alone for a given chemical composition, which is chosen to besolar (Pop. I) with X = 0.708, Y = 0.275, and Z = 0.017. The
constitutive equations account for the inner degrees of freedomof H2 , as well as the H2 dissociation, the ionisation of H andHe, 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 todust evaporation is clearly visible, as well as the “kappa mountain” primarily produced by the H− -ion and H-ionization.

J. Schönke and W. M. Tscharnuter: Protostellar collapse of rotating cloud cores2.4. Numerical techniques

2.6. Initial conditions

We use a finite-volume scheme on a staggered grid in sphericalcoordinates (r, θ, ϕ) to discretize the system of Eqs. (1)–(5). In
the (r, θ)-plane, we construct a self-similar, locally uniform, andstatic 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 operatorsplitting technique is applied to solve the equations. The firstsource 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 toevaluate the advection terms explicitly. A second-order van-Leerinterpolation 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 discretizedenergy equations, are solved by using a direct solver (in this casethe software library PARDISO1 ). Direct solvers have the advantage that they give the exact solution of the linear system (incontrast to all other solvers, which always yield approximate results), thus giving rise to a more accurate and robust overall numerical scheme.

The outcome of protostellar collapse calculations is highly sensitive to the initial parameters. Since we only begin to grasp thedetailed 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 thisobject, 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 withuniformly rotating Bonnor-Ebert spheres (BES) (see, e.g., Alveset al. 2001 or Lada et al. 2003).
We therefore use a slightly supercritical BES2 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 angularvelocity ω is assumed, given in terms of the centrifugal barrierof the mass elements farthest away from the rotation axis as

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 configurationsof 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 ofthe 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 fromequatorial symmetry.
All “inner” boundary conditions are determined by thesesymmetry 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 andremains 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 canonly 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 conditionis especially tailored to the system, we describe the details in thediscussion of the results in Sect. 3.9.
1

http://www.pardiso-project.org

Rc =4
ω2 rmax,
GM(11)

where rmax is the radius of the outermost sphere and M the total
cloud mass. The angular velocity itself is not a very convenientand descriptive parameter as it scales nonlinearly (quadratic)with the rotational energy. We choose Rc = 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. Walchet al. 2009) have already been carried out, there is still no “initialstandard configuration” for this kind of numerical study. Thereis no doubt that an agreement on this issue, as well as reasonable
parameter studies, would lead to a clearer understanding of thegreat variety of phenomena related to the formation of stars andprotostellar 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, thepressure-free collapse and Sod’s shock tube, we also tested theaccuracy of the solution of Poisson’s equation and problems including the diﬀusion 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 ofthem came out as expected. One aspect which goes certainly beyond a simple test problem is the local conservation of angularmomentum 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 . The2
We calculate the critical BES for a somewhat higher temperature(13 K) than the designated initial temperature (10 K) and use the resulting density profile together with the lower temperature to get an
unstable cloud.

A139, page 3 of 14

A&A 526, A139 (2011)

spatial resolution was chosen to be 320 grid points in the radialand 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 rmax ≈ 8700 AU. The ratio of the largest to the smallest
grid size is not less than 2×107! 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 thestellar core (because they are much more time consuming), they
are not shown here.
The accuracy during the calculations was monitored bychecking the total energy balance, which is given by
tot = int + rad + pot + kin + out =

 



t 
1
1 2
e + E + ρ Φ + ρ v dV +F · dS dt,
2
2
rmax

Fig. 2. Evolution in Phase 0. We plot the central density ρc (upper left),
central temperature T c (lower left), total luminosity L (upper right), andcentral optical depth τc (left-hand side of Eq. (14), lower right) of thecollapsing cloud.

(12)

0

3.4. Phase 0

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Δεrel (t) =tot (t) − tot (0)| pot (t) − pot (0)|

,

(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 firstcollapse (on a dynamical timescale) and the formation of an optically thick central object, or more precise, if the relation⎛ rmax
⎞
⎜⎜⎜ 
⎟⎟⎟
⎜⎜⎜
⎟ 2
min ⎜⎜κ ρ dr⎟⎟⎟⎟ ≥⎠ 3
θ, φ ⎝(14)

0

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 thetime to zero when (14) holds, therefore phase 1 starts at t = 0.
If the first hydrostatic core heated up significantly (becauseof accretion and contraction) so that by means of the dissociationof 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 temperaturein the first core reaches 2000 K.
A139, page 4 of 14

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 eﬀect. 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 thetemperature 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 ofphase 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 theincoming material with increasingly larger specific angular momentum (and corresponding greater centrifugal barrier), an oﬀ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 andtemperature distributions. The two rings are well separated bya clear gap between them. The disk has reached a radius of approximately 27 AU and already contains more than 43% of thetotal mass of the cloud core. The outer parts (>10 AU) are nearly
isothermal owing to the small optical depths. Because of the eﬃcient radiation of the large disk surface, the temperature stays below 1400 K, which is typical of dust evaporation. Conditions forthe dissociation of hydrogen molecules are certainly not reached.

J. Schönke and W. M. Tscharnuter: Protostellar collapse of rotating cloud coresFig. 3. Meridional cross-section (left border: rotation axis; bottom line:
equator). The density (top) and temperature distribution (bottom) in thedisk are shown after approximately 28 kyr of the evolution in phase 1
for the β = 0 model.

Fig. 5. Meridional cross-sections. We displayed the density (top) andtemperature distribution (bottom) in the disk after 7.7 kyr of the evolution in phase 1 for the β = 10−4 model.

about 15 kyr, only a minor redistribution of angular momentum
is visible in a region l ≈ 10−6 −10−4 lmax . After 25 kyr, we finda clear drop in the spectrum for l  10−3 lmax , which is produced during the second ring formation at about 17 kyr. Owing
to numerical diﬀusion too much angular momentum has been
advected away from the rotation axis in this process.
3.6. Simulation for β = 10 −4

Fig. 4. Mass spectrum of the specific angular momentum for the β = 0model. The evolution of M(l) is shown for the initial BES (thick solidline), at 0 yr (dotted line), at 15 kyr (dashed line) and 25 kyr (solid line).

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 angularmomentum (l = ωr2 sin2 θ) 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 axiallysymmetric 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

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 tooutwardly directed angular momentum transport. The inner regions have higher temperatures because of viscous heating. Afterabout 7.7 kyr, the conditions for the dissociation of hydrogen arereached, 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, triggeredby the hydrogen dissociation of the outer ring parts. By meansof this event, the simulation comes to a “natural” end, since thespatial 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 ournumerical scheme can indeed resolve, we discuss them in moredetail.
A139, page 5 of 14

A&A 526, A139 (2011)

Fig. 6. Meridional cross-sections. Displayed is the density (top) andtemperature distribution (bottom) in the disk after 2.4 kyr of the evolution in phase 1 for the β = 10−3 model.

3.7.1. First core and dust evaporation

The β = 10−3 model results in a more compact and warmer innerpart 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 coeﬃcient κρ (the inverse mean free path). The resulting
quasi-isothermal region in the dust-free part leads to a “blowup” of the disk stratification starting at about 2520 yr. As thetemperature 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 beeﬃciently radiated away in this direction. Like piercing a hole
into a balloon, the blown-up disk layers then begin to subside
again.

Fig. 7. Meridional cross-sections. Plotted is the absorption coeﬃcientκρ in units of AU−1 for the β = 10−3 case. Several snapshots of theevolution of the dust evaporation process are shown at 2410 yr (upper
left), 2520 yr (lower left), 2640 yr (upper right), and 2670 yr (lower
right).

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 theform of the strong vortex with speeds up to Mach 0.8, centeredabout 2 AU above the equatorial plane. The subsonic motionsin 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
A139, page 6 of 14

Fig. 8. Meridional cross-section for β = 10−3 . The non-stationary (r, θ)flow field with additional streamlines at 2680 yr is plotted. The absolutemagnitude of the velocities is given by the color scale in terms of theMach number |u|/cs . The vector length scales with the radius and has no
physical meaning.

(triggered by the hydrogen dissociation) sets in, visible throughthe radial inflow near the rotation axis at a height of 1 to 2 AU.
The subsequent evolution of the collapse flow, documented inFig. 9 and 10, conflicts with our current understanding of the

J. Schönke and W. M. Tscharnuter: Protostellar collapse of rotating cloud coresFig. 9. Meridional cross-section of the flow field (cf. Fig. 8) at 2690 yr
for β = 10−3 . Note the smaller display detail compared to Fig. 8.
Fig. 11. Meridional cross-sections for β = 10−3 . The temperature distribution for several snapshots of the protostar formation are shown at
2700 yr (upper left), 2700.7 yr (lower left), 2750 yr (upper right), and2936 yr (lower right). Note that the spatial units are solar radii.

Fig. 10. Meridional cross-section of the flow field (cf. Fig. 8) at 2700 yr
for β = 10−3 .

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 andmore cone-like shaped (see Fig. 10). The whole process takes
place on a thermal timescale, although the exact estimate of thistimescale is diﬃcult, 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

Fig. 12. Meridional cross-sections of the density distribution during
protostar formation for β = 10−3 . Plotted are the same times as in Fig. 11(2700, 2700.7, 2750, and 2936 yr).

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.
A139, page 7 of 14

A&A 526, A139 (2011)

About 10 yr before the onset of the second collapse, the outermost parts of the first core allow the absolute ratio of totalrotational to gravitational energy of the first core to exceed thecritical value of 0.274, where it is likely that the object becomesprone to the growth of non-axisymmetric structures. Since therotation period of the outer first core (which should be related
to the typical timescale for the growth of non-axisymmetricfeatures) is about 125 yr, there is apparently not enough timeto 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 coneshaped 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 thecentral parts through the protostellar disk in a radial inward flow
along the equator. This transport occurs on the viscous timescaleand 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 thenewborn star is characterized by a highly variable flow field inthis 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 fielddevelops and both the protostar (with a radius of 2–3 R ) andthe 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 consistsof 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 scalesof about 10 AU). Figure 13 shows the disintegration of the firstcore’s global structure in terms of the absorption coeﬃcient forphase 2 of the evolution. Although most of the accretion energy
is radiated away along the rotation axis in the optically thin anddust-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 diﬀusing radiation heats the gas and allows it to expand, hence the density is
reduced (and so is the absorption coeﬃcient κρ) permitting theradiation to diﬀuse into the gas even more easily. An interestingquestion is how this result would change if we applied a more realistic model of the radiation field instead of using the simplified
diﬀusion approximation.
To which impressively large scales the expansion of the firstcore leads in only about 240 yr of the protostar’s life can beseen in Fig. 14. A nearly spherical outflow, becoming supersonicalready after about 50 yr of the protostar formation, disperses
A139, page 8 of 14

Fig. 13. Meridional cross-sections. Plotted is the absorption coeﬃcientκρ in units of AU−1 for the β = 10−3 case. Several snapshots of theerosion of the optically thick parts of the first core are shown at 2700 yr(upper left), 2800 yr (lower left), 2850 yr (upper right), and 2936 yr(lower right).

the material roughly homogeneously. After 2936 yr (the end of
the simulation), the outward travelling shock wave has reached aradius of more than 80 AU. A very diﬀuse, low-mass disk with aradius 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 returnto the center can be generally confirmed, since the total binding
energy of the system is, of course, still negative. How the processof re-accretion appears in some detail will be investigated in thelast 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 diﬀer 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 andwe assume an explicit prescription of turbulent viscosity, a morerobust 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 inphase 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) isconnected 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 tothe first stages of the main accretion onto the second core with a

J. Schönke and W. M. Tscharnuter: Protostellar collapse of rotating cloud coresFig. 15. Evolution of several quantities during phase 2 of the calculation for β = 10−3 . Upper left: total luminosity (with magnified inlay).
Lower left: relative subsonic mass and “stellar” (ρ > 0.1 kg m−3 ) mass
(with magnified inlay). Upper right: maximum (and virtually identicalcentral) density. Lower right: maximum and central temperature.
Fig. 14. Meridional cross-sections. Plotted is the density distribution(note the larger spatial scale) for the β = 10−3 case. Plotted are several snapshots of the radial outflow of matter previously forming thefirst core. Times are the same as in Fig. 13 (2700 yr, 2800 yr, 2850 yr,
and 2936 yr).

very complex and time varying flow field (and accretion rate) in
the narrow “tunnel” along the rotation axis. Besides some moreor 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 theend 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 thelower left panel of Fig. 15. On the one hand the relative subsonicmass Msub , which consists of all subsonic mass elements inside
the global inflow from the cloud core3 . Otherwise, the relative“stellar” mass Mst , defined as the sum of all mass elements withdensities >0.1 kg m−3. Until the first 50 yr after the formationof the protostellar core, Msub is constantly increasing, reachingvalues of nearly 10%. As the outflow (initiated by the protostar)
becomes supersonic, Msub commences to decrease. Finally, at
2936 yr Msub (= 1.77%) is only somewhat greater than Mst (=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 Mst (defined in the same way ashere) given in Bate (2010), we find a somewhat diﬀerent evolution. The increase in Mst 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 ismuch higher in our calculation (about 6 × 10−5 M yr−1 ) thanBate’s findings of about 10−5 M yr−1 . A more detailed comparison with the protostellar masses in Bate (2010) is diﬃcult, since3
Regions at the outer boundary of the simulation are also subsonic,
and need to be excluded.

the runs of Mst given there seem to depend significantly on thespatial 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 and2.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 theaccretion front in the outermost layers of the star and the matterin 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 mostof the material at the end of the simulation (at 2936 yr), thedensity is noticeably lower than before in the region of the former first core. The outwardly traveling shock is clearly visible atabout 80 AU. In the temperature distribution, we see the abovementioned 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 eﬀective angular momentum transport) and a stronger viscous heating.
A139, page 9 of 14

A&A 526, A139 (2011)

Fig. 16. Midplane density distribution for various evolution times (β =10−3 ).

Fig. 18. Meridional cross-sections. Plotted is the absorption coeﬃcientκρ in units of AU−1 for the β = 10−2 case. Several snapshots of theevolution of the dust evaporation process are shown at 1320 yr (upper
left), 1340 yr (lower left), 1370 yr (upper right), and 1400 yr (lower
right).

Fig. 17. Midplane temperature distribution for various evolution times(β = 10−3 ).

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 diﬀerences from theβ = 10−3 calculation. Since the temperature increase is generally more rapid in this case, the blow-up of the disk stratificationand the onset of convection at the border of the opacity gap occurs far more rapidly. However, in this case there is no way forthe 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 eﬃciently as in thesame 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
A139, page 10 of 14

Fig. 19. Meridional cross-section for β = 10−2 . The non-stationary(r, θ)-flow field with additional streamlines at 1408 yr is plotted. The
absolute magnitude of the velocities is given by the color scale in termsof the Mach number |u|/cs . The vector length scales with the radius andhas no physical meaning.

of 3 (see Fig. 19). This demonstrates the massive impact that the
dust evaporation can have on the structure of the inner disk parts.

J. Schönke and W. M. Tscharnuter: Protostellar collapse of rotating cloud coresFig. 20. Meridional cross-sections for β = 10−2 . The temperature distribution for several snapshots of the protostar formation are shown at
1409.5 yr (upper left), 1410 yr (lower left), 1420 yr (upper right), and1482 yr (lower right). Note that the spatial units are solar radii.

Fig. 22. Meridional cross-sections. Plotted is the absorption coeﬃcientκρ in units of AU−1 for the β = 10−2 case. Several snapshots of the erosion of the optically thick parts of the first core are shown at 1409.5 yr(upper left), 1413 yr (lower left), 1420 yr (upper right), and 1440 yr(lower right).

β = 10−3 case. The radial inflow during the collapse now occurs from all directions, though significantly slower along theequator.
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 thedecreasing 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 atabout 1409 yr reforms an approximately spherical central object
plus a diﬀuse disk at the equator after only some 10 yr throughangular momentum transport (cf. Fig. 21). After about 70 yr, theaccretion has virtually stopped again for the same reason as inthe β = 10−3 case (see next section).
3.8.3. The protostar’s impact on the first coreFig. 21. Meridional cross-sections of the density distribution during
protostar formation for β = 10−2 . Plotted are the same times as in Fig. 20(1409.5, 1410, 1420, and 1482 yr).

3.8.2. Second collapse and protostar formation

The second collapse proceeds on the dynamical timescale (ofshorter 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

The accretion onto the protostar has an even stronger influenceon 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 andshort-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 coeﬃcient. At 1409.5 yr, the structure of the first core is still visible, including the “ice line” (aslight 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 amaximum, resulting in the farthest extended dust-free zone that
A139, page 11 of 14

A&A 526, A139 (2011)

Fig. 24. Evolution of several quantities during phase 2 of the calculation
for β = 10−2 . Upper left: total luminosity. Lower left: relative subsonicmass and “stellar” (ρ > 0.1 kg m−3 ) mass. Upper right: maximum (andvirtually identical central) density. Lower right: maximum and central
temperature.
Fig. 23. Meridional cross-sections. Plotted is the density distribution(note the larger spatial scale) for the β = 10−2 case. Plotted are several snapshots of the radial outflow of matter previously forming thefirst core. Times are (1410 yr, 1430 yr, 1450 yr, and 1482 yr).

allows the radiation to penetrate the disk even more easily. The
visible wave-like structure (in Fig. 22, lower left panel4 ) inside
the opacity gap develops because the dust-free region forms akind of cavity for the thermal waves. The border of the opacity gap reflects a part of the outgoing thermal waves becauseof 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 coeﬃcient 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 theformer 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. Startingwith the onset of the second collapse at about 1410 yr (Fig. 23
upper left panel), where the first core is only perturbed in thecentral 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. Thatthe outflow is not exactly spherical (as in the β = 10−3 case)may be because of the preceding super sonic blow-out producedby 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 corereaches on average about the same radial distances in the sametime for both simulations. But again, we see no indication of abipolar outflow. As said before, because of the diﬀerent initialconditions and the turbulent viscosity, we cannot compare theseresults in any greater depth.
4This phenomenon is visible in previous figures of the absorption coeﬃcient, too, but is most conspicuous here.

A139, page 12 of 14

3.8.4. Evolution of some global and midplane variables

As for the previous simulation, we summarize the evolutionof some global quantities (Fig. 24). The total luminosity run(Fig. 24, upper left panel) confirms the picture of a short andintense 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 ofthe quasi-stationary accretion flow that develops directly after ashort (<5 yr) and highly variable initial phase.
The relative subsonic (Msub ) and stellar (Mst ) mass, introduced in Sect. 3.7.4, are depicted in the lower left panel ofFig. 24. The first small drop in Msub at 1409 yr is related tothe dynamical collapse at the very beginning of the protostar’s
formation, where Mst gains about 0.7% within only 3 yr after the second collapse. At about 1412 yr, the blow-out of thefirst core reaches supersonic velocities and Msub 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 Mst given in Bate (2010), we again find acceptable agreement at least for the highest spatial resolution given
there (3 × 106 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 notyield a significantly diﬀerent outcome.
The maximum (and identical central) density evolution isshown 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 characterof the second collapse. The same is true for the maximum andcentral temperature evolution given in the lower right panel ofFig. 24. The reason that the central temperature falls more andJ. Schönke and W. M. Tscharnuter: Protostellar collapse of rotating cloud coresFig. 25. Midplane density distribution for various evolution times (β =10−2 ).
Fig. 27. Meridional cross-sections for the β = 10−2 case. Plotted is thelong time evolution of the blown-out first core in form of density distributions at 2 kyr (upper left), 4 kyr (lower left), 5 kyr (upper right), and6 kyr (lower right).

3.9.1. Technical realization of an inner boundary

Fig. 26. Midplane temperature distribution for various evolution times(β = 10−2 ).

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 regionof 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 oﬀ-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.

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 radiusof 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 diﬀuse,
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 (eventhough 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 ofthis artificial boundary would lengthen the section unreasonably,
hence for more details we refer to Schönke (2010)5.
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 ofmaterial that is still expanding supersonically (Fig. 27 upper left
panel, compare with the previous state in Fig. 23 lower rightpanel 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 backto the center. Interestingly, a large fraction of the material first
concentrates around the rotation axis and then falls back in a
5

The Ph.D. thesis can be downloaded from:http://archiv.ub.uni-heidelberg.de/volltextserver/
volltexte/2010/10720/pdf/dissertation.pdf
A139, page 13 of 14

A&A 526, A139 (2011)

Fig. 28. Continuation of Fig. 27 at 8, 10, 20, 56 kyr.

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 diﬀuse corewith scales of about 100 AU that re-arranges itself with time toa 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 firstcore has been re-accreted, and after that, again only “fresh” material from the outer cloud parts is accreted. At the end of thecalculation (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 thebeginning 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 blowout. 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 eﬀectivity 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 opticallythin path to the core’s surface (providing an eﬀective coolingmechanism, 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 importantand, to our knowledge, has not been pointed out so far in theliterature.
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 thecollapsing region the more dominant the dynamical timescale,i.e., whether the angular momentum transport is eﬀective enoughto reduce the ratio of centrifugal to gravitational forces suﬃciently, 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 pressure6 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 wehave 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 asa highly non-stationary and cyclic one. The short and luminous
accretion phases would then alternate with long (on the orderof several kyr) blow-out phases characterized by extremely low
luminosities. Future investigations will need to clarify whether
this is really the case.
Acknowledgements. This work has been supported by the Deutsche
Forschungsgemeinschaft within the Forschergruppe 759 “The Formationof Planets: The Critical First Growth Phase”. We thank the referee for theconstructive and thoughtful comments.

References
Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159
Bate, M. R. 1998, ApJ, 508, L95
Bate, M. R. 2010, MNRAS, 404, L79
Duschl, W. J., Strittmatter, P. A., & Biermann, P. L. 2000, A&A, 357, 1123
Hayes, J. C., Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188
Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626
Lada, C. J., Bergin, E. A., Alves, J. F., & Huard, T. L. 2003, ApJ, 586, 286
Larson, R. B. 1969, MNRAS, 145, 271
Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
Norman, M. L., Wilson, J. R., & Barton, R. T. 1980, ApJ, 239, 968
Saigo, K., Tomisaka, K., & Matsumoto, T. 2008, ApJ, 674, 997
Schönke, J. 2010, Ph.D. Thesis, Universität Heidelberg
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
Stamatellos, D., Whitworth, A. P., Bisbas, T., & Goodwin, S. 2007, A&A, 475,
37
Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2010, ApJ, 714, L58
Toomre, A. 1964, ApJ, 139, 1217
Tscharnuter, W. M. 1987, A&A, 188, 55
Tscharnuter, W. M., & Gail, H.-P. 2007, A&A, 463, 369
Tscharnuter, W. M., & Winkler, K. 1979, Comp. Phys. Commun., 18, 171
Tscharnuter, W. M., Schönke, J., Gail, H., Trieloﬀ, M., & Lüttjohann, E. 2009,
A&A, 504, 109
Walch, S., Burkert, A., Whitworth, A., Naab, T., & Gritschneder, M. 2009,
MNRAS, 400, 13
Wuchterl, G., & Tscharnuter, W. M. 2003, A&A, 398, 1081
Yorke, H. W. 1980, A&A, 86, 286
6

A139, page 14 of 14

Radiation pressure is negligible here.