© A. Bhandare et al. 2020

1. Introduction

Magnetized molecular clouds are known to be the birthplace of stars, which form by the gravitational collapse of dense, gaseous, and dusty cores within these clouds. Zooming in on the smallest scales in order to understand the complex physical processes such as hydrodynamics, radiative transfer, phase transition (in particular hydrogen dissociation), and magnetic fields has been challenging both theoretically and observationally (e.g. Nielbock et al. 2012; Launhardt et al. 2013; Dunham et al. 2014). Among various fundamental questions which remain to be answered are the comprehensive characterisation of the core collapse (see detailed reviews by Larson 2003; McKee & Ostriker 2007; Inutsuka 2012; Teyssier & Commerçon 2019; Pudritz & Ray 2019). It is thus crucial to unravel the role of various physical mechanisms involved during the transition of a molecular cloud core (i.e. pre-stellar core) to the hydrostatic Larson cores, to understand how stars form.

Several numerical studies using both grid-based (Bodenheimer & Sweigart 1968; Winkler & Newman 1980a,b; Stahler et al. 1980a,b, 1981; Masunaga et al. 1998; Masunaga & Inutsuka 2000; Tomida et al. 2010, 2013; Commerçon et al. 2011; Vaytet et al. 2012, 2013, 2018; Vaytet & Haugbølle 2017; Bhandare et al. 2018, and references therein) and smoothed particle hydrodynamics (SPH) methods (Bate 1998; Whitehouse & Bate 2006; Stamatellos et al. 2007; Bate et al. 2014; Wurster et al. 2018) have deduced the formation of a protostar to be a two-step process during which the non-homologous collapse leads to the presence of two quasi-hydrostatic cores, famously known as Larson’s cores (Larson 1969).

The initially isothermal, optically thin cloud core collapses under its own gravity due to efficient cooling via thermal emission from dust grains. The collapse may be initiated either by the effects of non-ideal magnetohydrodynamics (e.g. Shu et al. 1987; Mouschovias 1991) or the dissipation of turbulence, which reduces the effective sound speed in pre-stellar cores (e.g. Nakano 1998). It can also be triggered by an external shock wave crossing a previously stable cloud (Masunaga & Inutsuka 2000) or a marginally stable Bonnor–Ebert sphere may also undergo a collapse.

With time, the optical depth of the central collapsing region becomes greater than unity due to an increase in central density and the radiative cooling becomes inefficient. As the cloud core compresses, the temperature begins to rise in this dense region, which drives the first adiabatic collapse phase. This further leads to the formation of the first hydrostatic core, which eventually contracts adiabatically with an adiabatic index γ_actual ≈ 5/3, where γ_actual is the change in the slope of the temperature evolution with density (see Fig. 1). As the temperature rises (T ≳ 100 K), the rotational degree of freedom of the diatomic gas gets excited and the adiabatic index changes to ≈7/5. The process of formation and evolution of the first core lasts for about 10⁴ years, for the collapse of a 3000 au cloud core in marginal hydrostatic equilibrium, with an initial mass of 1 M_⊙ and 10 K temperature (e.g. Bhandare et al. 2018).

Fig. 1. Thermal evolution showing the first and second collapse phase for a 1 M⊙ cloud core. The change in adiabatic index γactual indicates the importance of using a realistic gas EOS.

Hydrogen molecules (H₂) begin to dissociate once the temperature reaches about 2000 K. This strongly endothermic process causes gravity to win over thermal pressure, which initiates the second collapse phase. During this phase the adiabatic index changes to ≈1.1, which is well below the stability limit of 4/3. The second hydrostatic core of atomic hydrogen is formed almost instantaneously once most of the H₂ is dissociated and is then followed by a phase of adiabatic contraction. An increase in the thermal pressure eventually halts the collapse. The second core grows in mass as the surrounding envelope continues to accrete onto the central core. A star is born once the centre of the core reaches ignition temperatures (T ≥ 10⁶ K) for nuclear fusion (i.e. deuterium and hydrogen burning).

Our goal is to perform isolated molecular core collapse simulations, which involve detailed thermodynamical modelling along with radiation transport. We include in the equation of state (EOS) the dissociation of H₂ and ionisation of atomic hydrogen and helium, and take the rotational and vibrational degrees of freedom into account. The numerical scheme and setup is described in Sect. 2.

Following the detailed analysis of the first core properties in Bhandare et al. (2018), we first discuss the second core properties for a wide range of initial cloud core masses spanning from 0.5 M_⊙ to 100 M_⊙ in Sect. 3.1. Going a step beyond the one-dimensional (1D) studies, we further investigate the collapse for the cases of 1 M_⊙, 5 M_⊙, 10 M_⊙, and 20 M_⊙ initial non-rotating cloud cores using two-dimensional (2D) radiation hydrodynamic (RHD) simulations with a resolution that has not been achieved before. We follow the evolution of the second hydrostatic core for ≥100 years after its formation, for each of these cases. Section 4 onward we focus on the results from our 2D simulations. In Sect. 4.1 we first describe the fiducial 1 M_⊙ case and discuss the properties of the second core. The highlight of our investigation is that we find indications for convection in the outer layers of the second hydrostatic core, triggered in these early stages of protostar formation. In contrast to fully convective stars, here, the energy is not generated at the stellar centre, but is provided by the accretion energy from outside the core.

Furthermore, the dependence of the second core properties on the initial cloud core mass are presented in Sect. 4.2. For each of these 2D cases, in Sect. 4.3, we discuss the occurrence of the standing accretion shock instability (SASI), which could describe the observed large-scale oscillations of the second accretion shock. We state the limitations of our method in Sect. 4.4 and provide comparisons with previous studies in Sect. 5. Section 6 summarises the results from both our 1D and 2D core collapse studies presented in this paper.

2. Numerical method

We perform 1D and 2D RHD calculations of a collapsing molecular cloud core using the (magneto) hydrodynamic code PLUTO (Mignone et al. 2007, 2012) combined with the self-gravity solver implemented and tested in Kuiper et al. (2010a, 2011). We consider the approximation of local thermodynamic equilibrium (LTE) and a two-temperature (2T) approach for the gas and radiation for gas thermodynamics. The equations of hydrodynamics that account for the mass, momentum, and energy conservation and the time-dependent radiation transport equation are described in Bhandare et al. (2018). We use the grey flux-limited diffusion (FLD) radiation transport module MAKEMAKE, the numerics of which is described in Kuiper et al. (2010b, 2019). It has been argued by Vaytet et al. (2012, 2013) that the grey method proves to be sufficient in the early evolutionary stages of the collapse and the multi-group radiative transfer may become more important only in the later stages of protostellar evolution.

Additionally, we use a gas EOS from D’Angelo & Bodenheimer (2013) to account for effects such as dissociation of H₂, ionisation of atomic hydrogen and helium, and molecular vibrations and rotations. This has been implemented in the PLUTO code by Vaidya et al. (2015) and updates to the radiation transport module to account for the realistic gas EOS are previously described in Bhandare et al. (2018).

Tabulated dust opacities are used from Ossenkopf & Henning (1994) whereas tabulated gas opacities are used from Malygin et al. (2014). In order to account for the contribution from dust dominating at low temperatures, we updated the evaporation and sublimation module to consider a time-dependent evolution of the dust as previously discussed in Bhandare et al. (2018).

We make use of a conservative finite volume approach based on second-order Godunov-type schemes, that means a shock capturing Riemann solver implemented in PLUTO to solve the equations of hydrodynamics. Our basic configuration for the flux computation consists of the Harten-Lax-Van Leer approximate Riemann solver that restores the middle contact discontinuity (hllc), a monotonised central difference (MC) flux limiter using piecewise linear interpolation and we integrate with a Runge-Kutta second order (RK2) method. On the other hand, the FLD equation is solved in an implicit way using a standard generalised minimal residual solver with approximations to the error from previous restart cycles (LGMRES). A relative convergence tolerance value of 10⁻¹⁰ in terms of temperature is used. More details about the open-source solver library PETSc (Portable, Extensible Toolkit for Scientific Computation) can be found in Balay et al. (1997, 2019a,b).

2.1. Initial setup

We used a stable Bonnor–Ebert (Ebert 1955; Bonnor 1956) sphere-like density profile as the initial density distribution, as also described in our previous work (Bhandare et al. 2018). Due to the fact that the hydrostatic equilibrium condition of a Bonnor–Ebert sphere does not have an analytical solution, the density profile ρ(r) has to be solved for numerically. The initial outer ρ_o and central ρ_c densities are determined as:

$$\begin{array}{cc}& {\rho }_{\mathrm{o}}=({\displaystyle \frac{1.18c_{\mathrm{s}0}^3}{M_0{G}^{3/2}}}{)}^2\hfill \\ \hfill & {\rho }_{\mathrm{c}}=14.1{\rho }_{\mathrm{o}},\hfill \end{array}$$(1)

where G is the gravitational constant, M₀ is the initial cloud core mass, and the initial sound speed c_s0 for a cloud core with radius R_cloud is computed as:

$$\begin{array}{c}\hfill c_{\mathrm{s}0}^2={\displaystyle \frac{G{M}_0}{\ln (14.1)R_{\mathrm{cloud}}}}·\end{array}$$(2)

The density contrast between the centre and edge of the sphere corresponds to a dimensionless radius of ξ = 6.45, where ξ is defined as:

$$\begin{array}{c}\hfill \xi =\sqrt{{\displaystyle \frac{4\pi G{\rho }_{\mathrm{o}}}{c_{\mathrm{s}0}^2}}}{R}_{\mathrm{cloud}}.\end{array}$$(3)

The integrated mass of the cloud core is the same as that of a critical Bonnor–Ebert sphere. The temperature T_BE for the stable sphere is computed as:

$$\begin{array}{c}\hfill T_{\mathrm{BE}}={\displaystyle \frac{\mu c_{\mathrm{s}0}^2}{\gamma \mathfrak{R}}},\end{array}$$(4)

where the mean molecular weight μ = 2.353, γ = 5/3, and ℜ is the universal gas constant.

The collapse is initiated by prescribing a thermal pressure using an initial temperature T₀ of 10 K instead of T_BE from a stable Bonnor–Ebert sphere setup, which allows gravity to dominate. The radiation temperature is set to be initially in equilibrium with the gas temperature. The dust and gas temperatures are treated as perfectly coupled throughout the simulation. We do not include any initial cloud rotation in the simulations presented herein.

2.2. One-dimensional setup

The 1D studies span a wide range of initial cloud core masses from 0.5 M_⊙ to 100 M_⊙. The spatial domain extends from 10⁻⁴ au up to 3000 au for all the different cases, which implies that for a constant initial temperature of 10 K, the central density of the different Bonnor–Ebert spheres scales as a function of initial cloud core mass (see Table 1). Convergence tests and effects of different initial cloud core properties have been discussed in our previous work (Bhandare et al. 2018).

We used 320 uniformly spaced grid cells from 10⁻⁴ au to 10⁻² au and 4096 logarithmically spaced grid cells from 10⁻² au to 3000 au, with identical sizes for the last uniform cell and the first logarithmic cell. A logarithmic binning throughout the whole domain would lead to extremely small grid cells resulting in much smaller time steps. The domain was thus scaled differently in the inner dense regions to prevent the simulations from being computationally very expensive. Our computational grid thus comprises of 4416 cells in total, covering a dynamical range of seven orders of magnitude. The smallest radial cell size Δx_min = Δr = 3.09 × 10⁻⁵ au. We used a minimum of 50 cells per Jeans length, which is estimated at the highest central density. Otherwise, we used 10³ − 10⁵ cells per Jeans length.

2.3. Two-dimensional setup

We expanded our 1D studies for a few cases of 1 M_⊙, 5 M_⊙, 10 M_⊙, and 20 M_⊙ cloud cores, thus accounting for cases from the low-, intermediate-, and high-mass regimes. For this, we adopted a 2D spherical Eulerian grid with axial and midplane symmetry. The grid comprises of 1445 logarithmically spaced cells in the radial direction extending from 10⁻² au to 3000 au, thus spanning a dynamical range covering five orders of magnitude. The logarithmic spacing increases resolution in the central parts of the computational domain. In the polar direction, we used 180 uniformly spaced cells stretching from the pole (θ = 0°) to the midplane (θ = 90°). The number of cells is tuned in order to ensure equal spatial extent in the radial and polar direction with the smallest cell size Δx_min = Δr = rΔθ = 8.77 × 10⁻⁵ au. We used a minimum of 49 cells per Jeans length, which is estimated at the highest central density. Otherwise, we used 300 − 10⁴ cells per Jeans length.

2.4. Boundary conditions

We used a reflective boundary condition at the inner radial edge for the hydrodynamics and a zero gradient condition for the radiation energy (i.e. no radiative flux can cross the inner boundary interface). At the outer radial edge, we used a Dirichlet boundary condition on the radiation temperature with a constant boundary value of T₀ and an outflow–no-inflow condition for the hydrodynamics that includes a zero-gradient (i.e. no force) boundary condition for the thermal pressure, the polar and the azimuthal velocity components. For the 2D runs, we used axisymmetric boundaries at the pole and mirror-symmetric boundaries at the equator.

3. Results: second core in spherical symmetry

In our previous work (Bhandare et al. 2018), we investigated the collapse of a molecular cloud core through the phase of the first hydrostatic core formation and tracked its evolution until the formation of the second core. In the studies presented herein, we follow the evolution of the second core for about 150 to 500 years after its formation, depending on the initial cloud core mass.

Figure 1 shows the different evolutionary stages during the collapse of a fiducial 1 M_⊙ cloud core in one of our simulation runs. The collapse of an initial isothermal molecular cloud core proceeds via a first phase of adiabatic compression and contraction, which forms the first hydrostatic core. This is then followed by a second collapse triggered due to dissociation of H₂ forming the second core, which undergoes another phase of adiabatic contraction as summarised in Sect. 1. The phase transitions at the different stages indicates the importance of using a realistic gas EOS as detailed in Sect. 2.

3.1. Dependence on initial cloud core mass

In order to investigate the dependence of the second core properties on the initial cloud core mass, we further span a wide range of masses from 0.5 M_⊙ to 100 M_⊙. Figure 2 shows the radial profiles of a variety of physical properties for all the different masses at a snapshot in time when the second core has evolved further and the first core no longer exists. Most of the material from the first core has been accreted on to the second core until this evolutionary stage. The accretion shock front in the radial velocity profile, which also coincides with the discontinuity in the density profile, defines the second core radius. This second shock is seen to be subcritical, that means the post-shock temperature is higher than the pre-shock temperature, suggesting that the accretion energy is transferred onto the second core and not radiated away (see discussions in Vaytet et al. 2013; Bhandare et al. 2018, and compare to the planetary case in Marleau et al. 2017, 2019).

Fig. 2. Radial profiles (across and down) of (a) density, (b) pressure, (c) gas temperature (=radiation temperature), (d) velocity, and (e) enclosed mass as well as the (f) thermal structure, (g) Mach number, (h) ratio of gas to ram pressure Pram = ρu2, (i) internal energy density as a function of temperature, (j) optical depth, (k) degree of ionisation (Bhandare et al. 2018, Eq. (16)), and (l) degree of dissociation (Bhandare et al. 2018, Eq. (17)) for all the cases at the final simulation snapshot (final times for different cases indicated in Table 2) in our simulations when the central density is roughly 0.5−0.8 g cm−3. Different colours indicate cloud cores with different initial masses as seen in the colour bar. Grey lines in the thermal structure plot show the 50% dissociation and ionisation curves.

We compare the different cases at a point in time when the central densities are around 0.5−0.8 g cm⁻³ and the central temperatures are roughly 10⁵ K. We note that the evolutionary timescales for the cloud cores with different initial masses are not the same, which indicates that high-mass cloud cores collapse faster than the low-mass ones. Most significant differences due to initial cloud core masses are visible outwards from the second core. The thermal structure indicates that the initially similar isothermal cloud cores eventually heat up at different densities. This has a significant effect on the formation timescale and the lifetime of the first and second cores.

The radial temperature profile shows an off-centred peak. The location of this peak corresponds to the radial position of the plateau seen in the density profile. This off-centred peak is also seen in studies by Masunaga & Inutsuka (2000, see their Fig. 4) and Tomida et al. (2013, see their Fig. 2); and also Winkler & Newman (1980a,b), Stahler et al. (1980a). Both these studies make use of the gas EOS by Saumon et al. (1995), which accounts for the effects of Fermi energy of the (partially) degenerate electrons. Masunaga & Inutsuka (2000) suggest that the off-centred peak seen in the temperature profile is due to partially degenerate electrons in the central region. In the off-centred region, the densities are not the highest and the thermal energy dominates the Fermi energy. There, the gas pressure is therefore more sensitive to the temperature. In the innermost regions, P_gas ≪ P_deg and hence there is no temperature rise. In other words, the off-centred temperature peak, for example, in Masunaga & Inutsuka (2000) can be interpreted as a temperature depression in the centre. However, since the effects of Fermi energy are not included in the (ideal) gas EOS used in our models (D’Angelo & Bodenheimer 2013), the reason for the off-centred peak in the temperature profile at the highest mass densities must be a different one in this case. The steady increase in density towards the centre of the second core due to the self-gravity of the core results in lower fraction of thermal ionisation (displayed in Fig. 2k). The associated release of energy yields the peak in the radial temperature profile. Furthermore, in this density–temperature regime, the gas also departs from being fully thermally dissociated in hydrogen (see inset in Fig. 2l), which possibly plays a role in producing this peak.

Towards the high-mass regime, the temperature profiles of the second core becomes flat (i.e. isothermal) where the shock around 0.1−1 au is hotter than in the core centre. In case this behaviour would hold as the core evolves further (in the high-mass case, accretion happens up to the ignition of hydrogen burning and longer), this will severely affect its internal evolution. The effect of the outer boundary conditions is visible in some of the profiles. For example, the sharp discontinuity seen in the temperature profiles for the high-mass cases is because the temperature at the outer boundary is set to a fixed value of 10 K. However, this does not have any significant effects on the evolution and properties of the hydrostatic cores on the smaller scales.

The bumps seen in the radial Mach number profiles for the high-mass cases are due to the behaviour of the adiabatic index Γ₁ during the dissociation and ionisation phase. In the high-mass regime, the temperature at the bump positions in the radial temperature profiles corresponds to the required temperatures for thermal hydrogen dissociation and ionisation. The radial velocity in the inner core regions, which are in hydrostatic equilibrium, fluctuates around the zero value. This effect is visible as the noise or spikes in the radial profiles of the ratio of gas to ram pressure in Fig. 2h. The radial profiles of the optical depth indicates that the second cores are optically thick, which makes it extremely difficult to detect them observationally and trace this evolutionary stage.

The region of fully atomic hydrogen in the dissociation profile extends far beyond the second core radius. The underlying reason is that the infalling gas in front of the second core radius (indicated by the shock) is already heated to temperatures beyond the dissociation temperature. This effect becomes especially clear for cases of higher accretion rates (i.e. higher accretion energy).

3.2. Second core properties

In this section, we discuss the dependence of the second core properties on the initial cloud core mass. Figure 3 shows the spatial evolution of the second core radius over a period of time from the onset of the second core formation until the central density reaches ≈0.5−0.8 g cm⁻³. The second core radius is defined using the position of the accretion shock in the velocity profile, which is similar to the position of the discontinuity or sharp rise in the density profile. In this work, the onset of formation of the second core is defined as the time when a prominent second accretion shock is visible in the velocity profile. The central density is greater than 10⁻² g cm⁻³ at this time snapshot in our simulations.

Fig. 3. Top: spatial evolution of the second core radius as a function of time for all the collapse scenarios with different initial cloud core masses ranging from 0.5 M⊙ to 100 M⊙ as indicated in the colour bar. The evolution is traced from the onset of the second core formation until the central density reaches ≈0.5−0.8 g cm−3. The inset in the upper left zooms in on the back and forth behaviour for one of the curves. This results from the jump between the two local minima in the velocity profile of the accretion shock, which is used to define the second core radius. Bottom: comparison of the Kelvin–Helmholtz and accretion timescales. A large (small) ratio is associated with an expansion (contraction) phase of the second hydrostatic core (see top panel).

In the high-mass regime, we find that, initially, the second core gradually expands with time, thus growing in size. This initial expansion occurs since the Kelvin–Helmholtz timescale¹ is much greater than the accretion timescale during this phase, as seen in the lower panel in Fig. 3. After reaching a maximum radius, the second core undergoes a phase of contraction during which the accretion timescale dominates. We note a similar behaviour in the low-mass regime where the collapse proceeds relatively slowly.

There is a back and forth behaviour seen in the evolution of the second core radius for the intermediate- and high-mass cases (see top panel in Fig. 3). This effect results from the jump between the two close local minima in the velocity profile of the accretion shock, which is used to define the second core radius. Another contributing factor for this behaviour are the small-scale oscillations of the second accretion shock. Both these effects are resolved due to a high time resolution and do not affect the overall behaviour of the second core radius.

For a more quantitative comparison, Fig. 4 shows the second core radius, accretion rate, and accretion luminosity as a function of the enclosed mass for different initial cloud core masses. These second core properties are displayed over a period of time from the onset of the second core formation until the central density reaches ≈0.5−0.8 g cm⁻³ as listed in Table 2. The accretion rate Ṁ_sc is estimated as:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{sc}}=4\pi R_{\mathrm{sc}}^2{\rho }_{\mathrm{sc}}{u}_{\mathrm{sc}},\end{array}$$(5)

Fig. 4. Second core radius (top), accretion rate (middle), and accretion luminosity (bottom) as a function of the enclosed mass for all the collapse scenarios with different initial cloud core masses ranging from 0.5 M⊙ (blue) to 100 M⊙ (red). The evolution is traced from the onset of the second core formation until the central density reaches ≈0.5−0.8 g cm−3.

where ρ_sc and u_sc are the density and velocity at the second core radius R_sc, respectively. The accretion luminosity L_acc computed using the accretion rate and mass M_sc enclosed within the second core radius is given as:

$$\begin{array}{c}\hfill L_{\mathrm{acc}}={\displaystyle \frac{G{M}_{\mathrm{sc}}{\dot{M}}_{\mathrm{sc}}}{R_{\mathrm{sc}}}}·\end{array}$$(6)

The second core mass gradually increases as the core evolves through the expansion and contraction phases. As expected, higher initial cloud core masses lead to more massive second cores. Initially, material from the first core accretes onto the second core at a much faster rate of ≈10⁻² M_⊙ yr⁻¹. The accretion rate slows down over time and can decrease to roughly a few times 10⁻⁵ M_⊙ yr⁻¹ in the low-mass end. The accretion luminosity in the high-mass regime is much higher than in the low-mass regime. Various properties of the second core are listed in Table 2.

4. Results: a 2D view of the second core

Results from our 1D studies of formation and evolution of the second hydrostatic core, for a wide range of initial cloud core masses, were discussed in Sect. 3.1. In this section, we further expand our investigation to 2D collapse for selected initial cloud core masses of 1 M_⊙, 5 M_⊙, 10 M_⊙, and 20 M_⊙. The main aim of this study is to resolve the second core using a resolution that has not been achieved before. The details of the computational grid and the resolution are discussed in Sect. 2.3.

The 2D simulations presented herein are basically a scaled-up version of the 1D runs; cloud rotation is not included as the cloud core is initialised as being at rest. We use the same initial conditions as discussed in Sect. 2.1 with a fixed outer radius of 3000 au and a constant initial temperature of 10 K. We account for the effects of self-gravity, radiation, and phase transitions on the evolution of these pre-stellar cores.

4.1. Evolution of a fiducial 1 M_⊙ pre-stellar core

In this section, we first focus on the fiducial 1 M_⊙ case, which evolves through the two stages of the first and second collapse. Figure 5 shows temperature snapshots at different stages of the evolution zooming into a 3000 au cloud core down to 10⁻² au (i.e. sub-au scales), thus covering five orders of magnitude in spatial scale.

Fig. 5. 2D temperature snapshots zooming down to sub-au scales showing the evolution of a 1 M⊙ cloud core with an initial temperature of 10 K and outer radius of 3000 au. The velocity streamlines in black indicate the infalling material and also the mixing within the second core at the last snapshot. The solid and dashed white contours indicate the first and second accretion shocks, respectively. Note the different spatial scales from top left to bottom left.

We follow the evolution of the second core for ≈312 years after its formation where the onset t = 0 is defined when a prominent second accretion shock is seen in the velocity profile. The central temperature at this stage is around 40 000 K. The gas and radiation temperatures are equal everywhere in our simulations.

Figure 6 shows the radial profiles of different properties of the cloud core at the final time snapshot (312 years after second core formation) of our simulation. The gradient from light to dark blue covers the polar angle range from the midplane (θ = 90°) to the pole (θ = 0°), respectively. We also compare this to the results from our 1D collapse studies of a 1 M_⊙ cloud core, which is indicated by the dashed red line in all the subplots in Fig. 6. The drop seen in the innermost part of the radial temperature profile is dependent on the inner radius and hence is affected by the inner boundary conditions. However, as discussed further in Appendix A, we confirm that this does not affect the second core properties nor violates energy conservation at the inner boundary.

Fig. 6. Radial profiles (across and down) at 312 years after formation of the second core, formed due to the collapse of a 1 M⊙ cloud core with an outer radius of 3000 au and an initial temperature of 10 K. The different subplots show the radial profiles (across and down) of (a) density, (b) pressure, (c) gas temperature, (d) radial velocity, and (e) enclosed mass as well as the (f) thermal structure, (g) Mach number, (h) ratio of gas to ram pressure Pram = ρu2, (i) internal energy density as a function of temperature, (j) optical depth, (k) Rosseland mean opacity, and (l) dissociation fraction. The colour gradient from light to dark blue spans the polar angle from the midplane (θ = 90°) to the pole (θ = 0°). The grey lines in the thermal structure plot show the 50% dissociation and ionisation curves. The radial profiles from the 1D collapse simulation for the same initial conditions and resolution are over-plotted as a dashed red line in all the subplots.

As expected, the same initial conditions lead to similar evolution in 1D and 2D. In both cases, the cloud core has evolved through the two phases of the first and second collapse, until a stage where the first core is no longer present and only the second core is visible, as indicated by the accretion shock in the radial velocity profile.

The four panels in Fig. 7, showing the 2D view of the second hydrostatic core (zooming into the inner 0.5 au), indicate the Mach number, density, temperature, and entropy structure of the second core. The infalling gas flow and internal mixing are indicated by the black velocity streamlines. The white contour in Fig. 7 (panel a) indicates Mach number equal to 1.0. This sets a clear separation between the supersonic outer region and the subsonic second core. This transition is also seen as the strong jump in the radial Mach profile in Fig. 6 at the second core radius. The contour lines in the density panel are labelled with density values at the different radial positions, marking an increase towards the centre, which confirms the non-homologous behaviour of the collapsing cloud core. Similar behaviour is seen in the temperature panel with the core centre having the highest value.

Fig. 7. 2D view of the second hydrostatic core at 312 years after its formation as a result of the collapse of a 1 M⊙ cloud core with an initial temperature of 10 K. The four panels show the (a) Mach number, (b) density, (c) temperature, and (d) entropy within the inner 0.5 au of the 3000 au collapsing cloud core. The velocity streamlines in black indicate the material falling on to the second core and the mixing inside the convective second core. The white contour in panel a showing Mach = 1.0 separates the sub- and supersonic regions. The different contour lines in panel b mark the increase in density towards the centre. When displayed in Adobe Acrobat, it is possible to switch to view the properties of the first core at the snapshot of the onset of the second core formation. A movie of the entire collapse is available online at https://keeper.mpdl.mpg.de/f/f04abdeabdf3472fb56d/.

We show the entropy panel again in Fig. 8 and plot it there with Line Integral Convolution to highlight the turbulent features within the second core. The entropy is calculated by the Sackur–Tetrode equation, which is consistent with the D’Angelo & Bodenheimer (2013) EOS used in the simulation and takes into account the molecular, atomic, and ionised hydrogen as well as the contribution from the electrons. It thus represents a straightforward extension of the expressions in Berardo et al. (2017, Appendix A), and details will be given in Marleau et al. (in prep.). We verified that the two agree well where ionisation (not included in Berardo et al. 2017) is not important.

Fig. 8. Left: line integral convolution visualisation of the entropy behaviour indicating the presence of eddies within the second core. Shown here is the inner 0.5 au of the 3000 au collapsing 1 M⊙ cloud core at 312 years after the formation of the second core. Right: radial entropy profile for the 1 M⊙ case, within the inner 1.0 au of the 3000 au collapsing cloud core at 312 years after the formation of the second core. The vertical grey dashed line indicates the radius of the second core. The colour gradient from light to dark blue spans the polar angle from the midplane (θ = 90°) to the pole (θ = 0°). The dashed red line shows the radial profile from the 1D collapse simulation for the same initial conditions and resolution, which by definition omits the effect of convection.

The entropy gradient is also seen in the radial profile in Fig. 8, zoomed into the inner 1 au. A convective instability is known to occur when a lower-entropy fluid lies above a higher-entropy fluid as seen in the region below the second accretion shock or the second core radius (see also left panel in Fig. 8). We interpret our results to mean that the instability is generated by the shock and grows radially inwards as the second core evolves over time (visualised in the movie²). The entropy is generated at the accretion shock and this yields a gradient from the high entropy at the shock towards the second core interior. Thus, convection is triggered in the outer layers of the second core at the accretion shock and drives the eddies inwards.

In Fig. 8, we also compare the radial entropy profile from our 1D (dashed red line) and 2D studies. Since there is no convection in 1D, entropy is generated and increases at the shock position. In comparison, since energy generated at the accretion shock is transported due to convection in the 2D case, the entropy profile flattens.

In order to further investigate this behaviour within the second core, we compare in Fig. 9 the actual ratio of the temperature and pressure gradients ∇_act to the adiabatic gradient ∇_ad, which is the gradient at constant entropy. In a (quasi-)static fluid, convective motions are expected if:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mathrm{act}}\equiv ({\displaystyle \frac{\mathrm{d}lnT}{\mathrm{d}lnP}}{)}_r>({\displaystyle \frac{\mathrm{d}lnT}{\mathrm{d}lnP}}{)}_{\mathrm{ad}}\equiv {\mathrm{\nabla }}_{\mathrm{ad}}.\end{array}$$(7)

Fig. 9. Polar-angle averaged actual temperature gradient ∇act(r) at 312 years after the formation of the second core in the 1 M⊙ case compared to the polar-angle averaged adiabatic gradient ∇ad(ρ(r),T(r)) (see Eq. (7)). Classically, regions with ∇act >  ∇ad are convectively unstable. The vertical dashed line indicates the radius of the second core.

We calculate ∇_ad according to Hansen et al. (2004, Eq. (3.98)) using simple differentiation of the P(ρ, T) function provided by the Vaidya et al. (2015) implementation of the D’Angelo & Bodenheimer (2013) EOS.

The dashed grey line in Fig. 9 indicates the radius of the second core. The region interior to this radius is convectively unstable. Due to the high resolution, the grid cell size in our simulations is roughly an order of magnitude smaller than the pressure scale height, thus allowing us to resolve the convection. A comparison between different resolutions is shown in Appendix B.1, which indicates the need to use such high resolution in order to resolve the eddies.

Convection allows mixing within a star and contributes by being an efficient means of heat transport. We find that at this stage in evolution, the energy flux is still dominated by radiation, however the convective flux can become stronger at later stages.

4.2. Dependence on initial cloud core mass

We further investigate the evolution of collapsing cloud cores with initial masses of 5 M_⊙, 10 M_⊙, and 20 M_⊙, thus covering a few cases in the intermediate- and high-mass regimes. The same initial temperature of 10 K and outer cloud core radius of 3000 au are used as in the 1 M_⊙ case. We study the effects of initial cloud core mass on the convective instability discussed in Sect. 4.1. Similar to the 1 M_⊙ case, for the 5 M_⊙, 10 M_⊙, and 20 M_⊙ runs, we find a turbulent pattern within the second core, indicated by the black velocity streamlines in Fig. 10 and seen as the eddies in Fig. 11. The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M_⊙, 10 M_⊙, and 20 M_⊙ cases, respectively.

Fig. 10. 2D view of the second hydrostatic core formed as a result of the collapse of a 5 M⊙ (top), 10 M⊙ (middle), and 20 M⊙ (bottom) cloud core with an initial temperature of 10 K. The four panels in each of the subplots show the (a) Mach number, (b) density, (c) temperature, and (d) entropy within the inner 0.5 au of an initial 3000 au cloud core. The velocity streamlines in black indicate the material falling onto the second core and the mixing inside the core. The white contour in panel a showing Mach = 1.0 separates the sub- and supersonic regions. The different contour lines in panel b show the increase in density towards the centre. The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M⊙, 10 M⊙, and 20 M⊙ cases, respectively.

Fig. 11. Line integral convolution visualisation of the entropy behaviour indicating the presence of eddies within the second core. Shown here is the inner 0.5 au of 3000 au collapsing cloud cores with different masses of 5 M⊙ (top), 10 M⊙ (middle), and 20 M⊙ (bottom). The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M⊙, 10 M⊙, and 20 M⊙ cases, respectively.

Figure 12 shows the comparison between the adiabatic index and the ratio of temperature and pressure gradients. On testing the criterion for convective instability stated in Eq. (7), we do not find a strong indication as seen in the 1 M_⊙ case. However, this may change as the second core evolves further for these cases. We are currently unable to further follow the evolution of the second core due to high computational expenses and this remains to be tested as part of future studies.

Fig. 12. Polar-angle averaged actual temperature gradient ∇act(r) compared to the polar-angle averaged adiabatic gradient ∇ad(ρ(r),T(r)) (see Eq. (7)) for the 5 M⊙ (top), 10 M⊙ (middle), and 20 M⊙ (bottom) initial cloud core masses. The indication for convective instability is not as prominent as in the 1 M⊙ case (see Fig. 9). The vertical grey dashed line indicates the radius of the second core. The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M⊙, 10 M⊙, and 20 M⊙ cases, respectively.

Figure 13 shows the evolution of the second core radius, mass, and accretion rate for all the cases with different initial cloud core mass. The time t = 0 marks the onset of the second core formation, which is indicated by a prominent second accretion shock as per the definition used herein. Higher initial cloud core mass leads to a faster collapse. Following the evolution, our results indicate that the 1 M_⊙ and 5 M_⊙ cloud cores will eventually form a protostar with mass less than 0.5 M_⊙. Stars within this mass range are known to be fully convective throughout their life.

Fig. 13. Time evolution of the second core radius (top), mass (middle), and accretion rate (bottom) for the different core collapse scenarios with initial cloud core mass of 1 M⊙ (red), 5 M⊙ (yellow), 10 M⊙ (blue), and 20 M⊙ (green). The time t = 0 marks the onset of the second core formation. The overplotted dashed lines indicate the evolution from the high-resolution 1D simulations discussed in Sect. 3.2.

4.3. Standing accretion shock instability

In this section, we further investigate the source of the turbulence visible in the post-shock regions. The standing (or spherical) accretion shock instability (SASI), known to play a crucial role in the explosion mechanism of core collapse supernovae, is seen to induce large-scale non-spherical oscillations of the shock (Foglizzo & Tagger 2000; Foglizzo 2002; Blondin et al. 2003; Foglizzo et al. 2006; Guilet & Foglizzo 2012). There are two proposed mechanisms that lead to the linear growth of this instability, one being the interplay between advected entropy–vorticity perturbations and acoustic waves (Foglizzo et al. 2007; Foglizzo 2009), whereas the second is a purely acoustic mechanism, which assumes that the trapped acoustic waves can be amplified by the shock (Blondin & Mezzacappa 2006).

During the evolution of the second core in our 2D collapse simulations, we observed some non-spherical, large-scale oscillations of the accretion shock front in the 1 M_⊙, 5 M_⊙, and 10 M_⊙ cases and comparatively small-scale oscillations of the accretion shock front for the 20 M_⊙ case. In their study, Scheck et al. (2008) have reported that large-amplitude SASI oscillations produce strong variations in the entropy, which can drive convective instability in the supernova core.

We further investigate the presence of the SASI in the case of protostellar cores and its role to generate turbulence behind the shock for all the different collapse scenarios. Large-scale, non-spherical oscillations of the second accretion shock front are indicated by the black line in Fig. 14. For a quantitative analysis of the SASI, following Scheck et al. (2008), Fig. 14 shows the advected perturbations in terms of the amplitudes of the largest modes of the spherical harmonics of the quantity A(t, r, θ) given by:

$$\begin{array}{c}\hfill A(t,r,\theta )={\displaystyle \frac{1}{\sin \theta }}{\displaystyle \frac{\partial }{\partial \theta }}(v_{\theta }(t,r,\theta )\sin \theta ).\end{array}$$(8)

Fig. 14. Time evolution of the amplitude of the dominant spherical harmonics mode of the quantity A(t, r, θ) from Eq. (8). The second core radius is shown by the black line. The time t = 0 indicates the onset of the second core formation for the different collapse scenarios. Shown here is a small interval in time for the 1 M⊙ (top left), 5 M⊙ (top right), 10 M⊙ (bottom left), and 20 M⊙ (bottom right) cases.

The term r⁻¹A is the divergence of the lateral velocity component. Several works have shown that the SASI can be measured more easily by determining A even for lower amplitudes of the instability (Scheck et al. 2008; Blondin & Mezzacappa 2006). For all the cases with different initial cloud core masses, we plot the amplitude for a small time interval since it helps to view the perturbations better. Although the amplitudes show more complex patterns than in Scheck et al. (2008, see their Fig. 12), there are some noticeable advected trajectories as well as some acoustic feedback. The important acoustic feedback timescale is given by the sound crossing time from the centre of the second core to the accretion shock and back.

We thus conclude that the SASI may not be operating as strongly as seen in the supernovae core-collapse studies and the convective instability seems to be the main source of the turbulent cells seen in the post-shock regions. However, the SASI could still be responsible for the large scale oscillations of the accretion shock front seen during the evolution of the second core. Nonetheless, it is interesting to note that the SASI can operate in different regimes.

4.4. Limitations

The simulations discussed in this paper only include the effects of self-gravity, radiation transport, dissociation, an1.5d ionisation on the core properties, but not due to rotation. Our ongoing studies which account for the effects of initial cloud rotation and magnetic fields on the hydrostatic cores and early discs will be reported in a follow-up article. Magnetic fields would indeed affect the formation and evolution of the second core. However, convective instability within the second core, at least for the low-mass end, could still be generated during the core evolution.

The 2D simulations do not stay spherically symmetric. Hence, the evolution of the convective second core will also not stay axially symmetric. Unfortunately, a 3D model achieving the same resolution as in the axial and midplane symmetric 2D simulations presented here is unfeasible at the moment.

5. Comparisons with previous work

In this section we first compare the second core properties from our 1D simulations to some of the previous studies for the case of a collapsing 1 M_⊙ cloud core. In their study, Masunaga & Inutsuka (2000) use a uniform initial density profile with ρ = 1.415 × 10⁻¹⁹ g cm⁻³ and an outer radius of 10⁴ au. A few years³ after the formation of the second core when the central density reaches ≈1 g cm⁻³, they find the second core radius to be ≈4 R_⊙ and the second core mass as 0.73 M_⊙. On the other hand, at central density greater than 10⁻¹ g cm⁻³, Tomida et al. (2013) report a bigger second core radius of ≈10 R_⊙ enclosing a mass of 2 × 10⁻² M_⊙ within 0.7 years after its formation⁴. They suggest that the second core continues to expand during the main accretion phase. Their simulations are stopped once the central temperature reaches 10⁵ K. In their work, they adopt a Bonnor–Ebert sphere-like density profile with an initial central density of 1.2 × 10⁻¹⁸ g cm⁻³ and an outer radius of 8800 au.

In Fig. 15, we compare the behaviour of the density and temperature from our 1D simulations (bluish purple) to those by Vaytet & Haugbølle (2017, dashed red), at a snapshot when the central density is roughly 10⁻¹ g cm⁻³, since their method is closest to ours. Both studies use the same initial conditions for the collapse of a 1 M_⊙ cloud core with an initial temperature of 10 K, outer radius of 3000 au, and an initial Bonnor–Ebert sphere-like density profile. We note some differences in the profiles which arise due to the different gas EOS (Saumon et al. 1995 used by Vaytet & Haugbølle 2017 versus D’Angelo & Bodenheimer 2013 used in this work), opacities, and grid schemes (Lagrangian versus Eulerian). Moreover, both simulations are not compared at the same time in evolution (see also discussion in Bhandare et al. 2018, Sect. 4.5). On comparing both studies, these differences in the numerical methods also lead to discrepancies in the second core properties, for example, in the second core radius and enclosed mass. Vaytet & Haugbølle (2017) report a smaller second core radius of ≈1 R_⊙ with an enclosed mass of 2.62 × 10⁻³ M_⊙. They expect the second core to grow in size due to further heating and mass accretion from the infalling envelope. In this work, at a similar central density of 9.6 × 10⁻² g cm⁻³, we find the second core radius to be ≈4 R_⊙ with an enclosed mass of 5.12 × 10⁻³ M_⊙.

Fig. 15. Comparisons of our results for an initial 1 M⊙ cloud core indicated in bluish purple to those by Vaytet & Haugbølle (2017) shown using dashed red line. Radial entropy (top), density (middle), and temperature (bottom) profiles are compared at the time when central density ρc in both simulations reach roughly 10−1 g cm−3.

In our studies, we follow the evolution of the second core for 188.2 years after the second core formation (see Fig. 3). At the final simulation time snapshot, our results indicate that as the central density reaches 0.85 g cm⁻³ the second core radius grows to be much bigger ≈54 R_⊙ with an enclosed mass of 4.40 × 10⁻² M_⊙. As demonstrated in Fig. 3, the expansion and contraction of the evolving second Larson core is controlled by the timescale ratio of Kelvin–Helmholtz contraction versus accretion.

Figure 15 also shows the comparison of the radial entropy profile. The two peaks seen in the entropy correspond to the positions of the first and second accretion shocks in both 1D studies. In Appendix B.2, we discuss the dependence of entropy on the numerical resolution. Furthermore, in Appendix C we discuss the change in the entropy profile as the cloud core evolves beyond the formation of the second core in our 2D simulations.

Schönke & Tscharnuter (2011) followed the collapse of a 1 M_⊙ cloud core using grid-based 2D RHD simulations for up to 240 years after the formation of the second core. They included an initial uniform rotation of their cloud core and could hence also investigate the early phases of disc formation. In order to evolve the system for a longer duration, they replaced the physical domain within 0.7 au with a sink prescription once the second core reached a quasi-static state. Their main goal was to investigate the effect of hydrodynamically driven turbulence using a β-viscosity prescription. In their models with β = 10⁻³ and β = 10⁻² they have described dynamically unstable layers as a consequence of dust evaporation in the central regions within the inner 3 au. They also indicate the occurrence of convection seen via temperature gradients and the presence of a strong vortex in this innermost region.

In the 2D studies presented herein, we observe some short-lived unstable regions within the first accretion shock during the evolution of the first core but without any prominent vortices or a convective instability. However, as already discussed in Sect. 4.1, we observe convection in the outer layers of the second core, which eventually evolves to become the protostar. The ability to follow the evolution for 312 years after the formation of the second core allows us to trace the evolution of the eddies and we find the convective instability to grow radially inwards from the second shock as the second core evolves over time.

6. Summary and conclusions

The collapse of a molecular cloud core proceeds through an initial isothermal phase leading to the formation of the first hydrostatic core, which undergoes an adiabatic contraction phase. This is followed by the second collapse phase triggered by the dissociation of H₂ once the central temperature rises above 2000 K. The second hydrostatic core is formed as a result of this process once most of the H₂ is dissociated.

In this work, we investigate the gravitational collapse of molecular cloud cores using 1D and 2D RHD simulations. We include self-gravity and radiation transport. Additionally, the gas EOS takes into account rotational and vibrational degrees of freedom for the H₂ molecules, which start being excited as the cloud core transitions from being effectively monatomic to diatomic, as well as their dissociation and ionisation.

For the 1D studies, we model cloud cores with an initial constant temperature of 10 K, a fixed outer radius of 3000 au, and initial cloud core masses that span a wide range from 0.5 M_⊙ to 100 M_⊙. We further expand our collapse studies to 2D with an identical initial setup as in the 1D runs. We model 3000 au non-rotating cloud cores with initial masses of 1 M_⊙, 5 M_⊙, 10 M_⊙, and 20 M_⊙, thus covering a few cases in the low-, intermediate-, and high-mass regime. We now summarise our key findings from both the 1D and 2D studies, which focus on the formation and evolution of the second Larson core.

Our 1D studies indicate that the cloud cores with a higher initial mass, collapse faster and form bigger, more massive second cores. We describe the dependence of the second core properties such as the radius, mass, accretion rate, and accretion luminosity on the initial cloud core mass. We discuss the expansion and contraction of the evolving second Larson core, which is controlled by the timescale ratio of Kelvin–Helmholtz contraction versus accretion. The accretion rate in the high-mass regime is found to be much higher than in the low-mass end, as expected. Here we have investigated cases in which the higher-mass molecular cloud cores are gravitationally more unstable than in the low-mass regime. A parameter study using different initial cloud core temperatures and outer radii for various initial cloud masses can be found in our previous work (Bhandare et al. 2018). The results presented herein are consistent with previous core collapse studies in the low-mass regime.

Circumstellar discs form as a consequence of conservation of angular momentum around stars. Currently, evolving the disc for a longer time, that means until the Class 0 phase, is often hindered due to time step restrictions and hence most studies replace the central (second) core with a sink particle. The influence of a sink particle on disc formation has been extensively discussed in Wurster & Li (2018). Data from the 1D studies presented here can be used as a lookup-table to compute the evolution of the central object (i.e. the protostar) for a longer duration, within a sink-cell paradigm.

Using our 2D setup for the four non-rotating collapse cases, we follow the evolution of the second core for ≥100 years after its formation. For the 1 M_⊙ case, we follow the evolution of the second core for 312 years after its formation. Our 2D studies show that the accretion shock leads to a convective instability in the outer layers of the second core, which grows radially inward over time. Due to the high resolution in our simulations (≈10 cells per pressure scale height below the shock), that has not been achieved before, we can resolve the convective cells for the first time. For the 1 M_⊙ case, we find convection being driven from the accretion shock towards the interior of the second Larson core. In contrast to fully convective stars, here, the energy is not generated at the stellar centre, but is provided by the accretion energy from outside the core. Investigating the evolution from these early convective phases due to accretion up to fully convective low-mass stars due to hydrogen burning remains a challenging task for future research in stellar physics.

The origin of magnetic fields in low-mass stars is still a matter of debate. Several studies speculate that the fields are either dominated by primordial or fossil-fields or are replaced by dynamo-generated fields within the first 100 years of evolution. Since young low-mass stars are observed to have strong (> kilogauss) magnetic field strengths, the likelihood of a fossil field could be excluded for cases where the magnetic field amplitude in the second core at birth is found to be less than a kilogauss. In this work, since we have already observed convection in the outer layers of the second hydrostatic core, further evolution may enable the generation of a convective-dynamo (Chabrier & Küker 2006). This would support the interesting possibility that dynamo-driven magnetic fields may be generated during this very early phase of low-mass star formation.

We note that the simulations presented here do not account for the effects due to initial cloud rotation and magnetic fields. Even though both of these will affect the evolution of the second core and its properties, we predict that convection seen in our studies should still be generated during this collapse phase in the low-mass regime.

Acknowledgments

We thank the anonymous referee for constructive and insightful comments that helped improve the quality of this manuscript. A.B. would like to thank Benoît Commerçon and Hubert Klahr for useful discussions during this work. We would like to thank A. Mignone and collaborators for their work in developing the open-source PLUTO code. A.B. is grateful to Thomas Müller (Haus der Astronomie, MPIA) for creating the 1 M_⊙ core collapse movie available online at https://keeper.mpdl.mpg.de/f/f04abdeabdf3472fb56d/ and for his help with the line integral convolution visualisation. R.K. acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundation (DFG) under grant no. KU 2849/3-1 and KU 2849/3-2. T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. G-D.M. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1) and from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. The simulations were performed on the ISAAC cluster at the Max Planck Center for Data and Computing in Garching.

References

Appendix A: Comparisons for different inner radii

Fig. A.1. Radial profiles of the density (top), velocity (middle), and gas temperature (bottom) for an initial 1 M⊙ collapsing cloud core from the 1D simulations are shown at a time step after the second core formation. The yellow line indicates an inner radius of 10−4 au and the dashed red line indicates an inner radius of 10−2 au.

We used an inner radius R_in of 10⁻² au for the simulations discussed in this paper. We note the decrease in temperature at the inner boundary seen in Fig. 6 for both the 1D and 2D studies. Due to the high computational expenses, we currently cannot perform tests with an inner radius less than 10⁻² au in 2D. However, given that there are no significant differences in the second core properties between the 1D and 2D simulations with the same initial setup, we compare the 1D results for runs with two different inner radii of 10⁻⁴ au and 10⁻² au. We do not see a drop in innermost regions for the collapse case with R_in = 10⁻⁴ au. We thus conclude that the decrease seen in case of R_in = 10⁻² au could be a numerical artefact due to the inner boundary being much closer to the second accretion shock. Besides the temperature decrease in the innermost central region, we do not find any significant differences in the second core properties when comparing both R_in cases. We also confirm that energy conservation is not violated at the inner boundary.

Appendix B: Convergence tests

B.1. Resolution tests for two-dimensional simulations

For an initial 1 M_⊙ cloud core, we perform collapse simulations using three different resolutions shown in Figs. B.1 and B.2. In both the figures, we indicate the results at a snapshot in time when the central densities are roughly similar. Polar-angle averaged radial density, velocity, and temperature profiles for the three different resolution runs in Fig. B.1 do not show any significant differences in the second core properties.

Fig. B.1. Polar-angle averaged radial profiles of the density (top), velocity (middle), and gas temperature (bottom) for an initial 1 M⊙ collapsing cloud core at an initial temperature T0 of 10 K are shown at a time step after the second core formation. The different lines indicate the results using various grid resolutions in the 2D simulations.

Fig. B.2. Line integral convolution visualisation of the second core formed from the collapse of a 1 M⊙ cloud core at an initial temperature T0 of 10 K and an outer radius of 3000 au. Top to bottom panels: results from runs using different resolutions in an increasing order as indicated in the legends. The entropy behaviour is shown at a snapshot when the central densities are roughly similar. For the results presented in this paper, we use the highest resolution, which allows us to resolve the convective eddies within the second core.

The main aim of this test was to highlight the importance of using a higher resolution in order to better resolve the eddies indicating a convective instability as discussed in Sect. 4.1. Resolution differences within the second core are more prominent in the 2D entropy plots seen in Fig. B.2. Computational time restrictions prevent us from using an even higher resolution to test the convergence.

B.2. Dependence of entropy on resolution

In Fig. B.3 we compare, for two different resolutions, the radial density, temperature, and entropy profiles from our 1D simulations of the collapse of a 1 M_⊙ cloud core at an initial temperature T₀ of 10 K and outer radius of 3000 au. The peak in the entropy profile corresponds to the position of the second accretion shock. The two runs show convergence with no significant differences in the behaviour. This suggests that a radial resolution of 1445 cells used in this study is sufficient.

Fig. B.3. Radial profiles of the density (top), gas temperature (middle), and entropy (bottom) for an initial 1 M⊙ collapsing cloud core at an initial temperature T0 of 10 K and an outer radius of 3000 au are shown. The yellow and red lines indicate the results using two different grid resolutions in the 1D simulations.

Appendix C: Entropy

Figure C.1 shows the spatial and temporal evolution of the polar-angle averaged radial entropy profile from the 2D simulation of a 1 M_⊙ collapsing cloud core with an initial temperature T₀ of 10 K and an outer radius of 3000 au. The two peaks in the entropy profile at earlier time snapshots (dashed red and yellow lines) are seen at the positions of the first and second accretion shocks. As the cloud evolves further, material from the first core is accreted onto the second core and the first shock disappears. Figure 15 shows the comparison of the entropy profile, at an earlier time snapshot, from the 1D studies presented herein with those by Vaytet & Haugbølle (2017).

Fig. C.1. Polar-angle averaged radial entropy profiles are shown at four different time snapshots as a 1 M⊙ cloud core transitions through the formation and evolution of the first and second hydrostatic cores.

All Tables

All Figures

	Fig. 1. Thermal evolution showing the first and second collapse phase for a 1 M⊙ cloud core. The change in adiabatic index γactual indicates the importance of using a realistic gas EOS.
In the text

Fig. 2. Radial profiles (across and down) of (a) density, (b) pressure, (c) gas temperature (=radiation temperature), (d) velocity, and (e) enclosed mass as well as the (f) thermal structure, (g) Mach number, (h) ratio of gas to ram pressure Pram = ρu2, (i) internal energy density as a function of temperature, (j) optical depth, (k) degree of ionisation (Bhandare et al. 2018, Eq. (16)), and (l) degree of dissociation (Bhandare et al. 2018, Eq. (17)) for all the cases at the final simulation snapshot (final times for different cases indicated in Table 2) in our simulations when the central density is roughly 0.5−0.8 g cm−3. Different colours indicate cloud cores with different initial masses as seen in the colour bar. Grey lines in the thermal structure plot show the 50% dissociation and ionisation curves.

In the text

Fig. 3. Top: spatial evolution of the second core radius as a function of time for all the collapse scenarios with different initial cloud core masses ranging from 0.5 M⊙ to 100 M⊙ as indicated in the colour bar. The evolution is traced from the onset of the second core formation until the central density reaches ≈0.5−0.8 g cm−3. The inset in the upper left zooms in on the back and forth behaviour for one of the curves. This results from the jump between the two local minima in the velocity profile of the accretion shock, which is used to define the second core radius. Bottom: comparison of the Kelvin–Helmholtz and accretion timescales. A large (small) ratio is associated with an expansion (contraction) phase of the second hydrostatic core (see top panel).

In the text

	Fig. 4. Second core radius (top), accretion rate (middle), and accretion luminosity (bottom) as a function of the enclosed mass for all the collapse scenarios with different initial cloud core masses ranging from 0.5 M⊙ (blue) to 100 M⊙ (red). The evolution is traced from the onset of the second core formation until the central density reaches ≈0.5−0.8 g cm−3.
In the text

Fig. 5. 2D temperature snapshots zooming down to sub-au scales showing the evolution of a 1 M⊙ cloud core with an initial temperature of 10 K and outer radius of 3000 au. The velocity streamlines in black indicate the infalling material and also the mixing within the second core at the last snapshot. The solid and dashed white contours indicate the first and second accretion shocks, respectively. Note the different spatial scales from top left to bottom left.

In the text

Fig. 6. Radial profiles (across and down) at 312 years after formation of the second core, formed due to the collapse of a 1 M⊙ cloud core with an outer radius of 3000 au and an initial temperature of 10 K. The different subplots show the radial profiles (across and down) of (a) density, (b) pressure, (c) gas temperature, (d) radial velocity, and (e) enclosed mass as well as the (f) thermal structure, (g) Mach number, (h) ratio of gas to ram pressure Pram = ρu2, (i) internal energy density as a function of temperature, (j) optical depth, (k) Rosseland mean opacity, and (l) dissociation fraction. The colour gradient from light to dark blue spans the polar angle from the midplane (θ = 90°) to the pole (θ = 0°). The grey lines in the thermal structure plot show the 50% dissociation and ionisation curves. The radial profiles from the 1D collapse simulation for the same initial conditions and resolution are over-plotted as a dashed red line in all the subplots.

In the text

Fig. 7. 2D view of the second hydrostatic core at 312 years after its formation as a result of the collapse of a 1 M⊙ cloud core with an initial temperature of 10 K. The four panels show the (a) Mach number, (b) density, (c) temperature, and (d) entropy within the inner 0.5 au of the 3000 au collapsing cloud core. The velocity streamlines in black indicate the material falling on to the second core and the mixing inside the convective second core. The white contour in panel a showing Mach = 1.0 separates the sub- and supersonic regions. The different contour lines in panel b mark the increase in density towards the centre. When displayed in Adobe Acrobat, it is possible to switch to view the properties of the first core at the snapshot of the onset of the second core formation. A movie of the entire collapse is available online at https://keeper.mpdl.mpg.de/f/f04abdeabdf3472fb56d/.

In the text

Fig. 8. Left: line integral convolution visualisation of the entropy behaviour indicating the presence of eddies within the second core. Shown here is the inner 0.5 au of the 3000 au collapsing 1 M⊙ cloud core at 312 years after the formation of the second core. Right: radial entropy profile for the 1 M⊙ case, within the inner 1.0 au of the 3000 au collapsing cloud core at 312 years after the formation of the second core. The vertical grey dashed line indicates the radius of the second core. The colour gradient from light to dark blue spans the polar angle from the midplane (θ = 90°) to the pole (θ = 0°). The dashed red line shows the radial profile from the 1D collapse simulation for the same initial conditions and resolution, which by definition omits the effect of convection.

In the text

	Fig. 9. Polar-angle averaged actual temperature gradient ∇act(r) at 312 years after the formation of the second core in the 1 M⊙ case compared to the polar-angle averaged adiabatic gradient ∇ad(ρ(r),T(r)) (see Eq. (7)). Classically, regions with ∇act >  ∇ad are convectively unstable. The vertical dashed line indicates the radius of the second core.
In the text

Fig. 10. 2D view of the second hydrostatic core formed as a result of the collapse of a 5 M⊙ (top), 10 M⊙ (middle), and 20 M⊙ (bottom) cloud core with an initial temperature of 10 K. The four panels in each of the subplots show the (a) Mach number, (b) density, (c) temperature, and (d) entropy within the inner 0.5 au of an initial 3000 au cloud core. The velocity streamlines in black indicate the material falling onto the second core and the mixing inside the core. The white contour in panel a showing Mach = 1.0 separates the sub- and supersonic regions. The different contour lines in panel b show the increase in density towards the centre. The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M⊙, 10 M⊙, and 20 M⊙ cases, respectively.

In the text

Fig. 11. Line integral convolution visualisation of the entropy behaviour indicating the presence of eddies within the second core. Shown here is the inner 0.5 au of 3000 au collapsing cloud cores with different masses of 5 M⊙ (top), 10 M⊙ (middle), and 20 M⊙ (bottom). The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M⊙, 10 M⊙, and 20 M⊙ cases, respectively.

In the text

Fig. 12. Polar-angle averaged actual temperature gradient ∇act(r) compared to the polar-angle averaged adiabatic gradient ∇ad(ρ(r),T(r)) (see Eq. (7)) for the 5 M⊙ (top), 10 M⊙ (middle), and 20 M⊙ (bottom) initial cloud core masses. The indication for convective instability is not as prominent as in the 1 M⊙ case (see Fig. 9). The vertical grey dashed line indicates the radius of the second core. The plots are shown at 128 years, 91.4 years, and 86.4 years after the second core formation, for the 5 M⊙, 10 M⊙, and 20 M⊙ cases, respectively.

In the text

Fig. 13. Time evolution of the second core radius (top), mass (middle), and accretion rate (bottom) for the different core collapse scenarios with initial cloud core mass of 1 M⊙ (red), 5 M⊙ (yellow), 10 M⊙ (blue), and 20 M⊙ (green). The time t = 0 marks the onset of the second core formation. The overplotted dashed lines indicate the evolution from the high-resolution 1D simulations discussed in Sect. 3.2.

In the text

Fig. 14. Time evolution of the amplitude of the dominant spherical harmonics mode of the quantity A(t, r, θ) from Eq. (8). The second core radius is shown by the black line. The time t = 0 indicates the onset of the second core formation for the different collapse scenarios. Shown here is a small interval in time for the 1 M⊙ (top left), 5 M⊙ (top right), 10 M⊙ (bottom left), and 20 M⊙ (bottom right) cases.

In the text

	Fig. 15. Comparisons of our results for an initial 1 M⊙ cloud core indicated in bluish purple to those by Vaytet & Haugbølle (2017) shown using dashed red line. Radial entropy (top), density (middle), and temperature (bottom) profiles are compared at the time when central density ρc in both simulations reach roughly 10−1 g cm−3.
In the text

	Fig. A.1. Radial profiles of the density (top), velocity (middle), and gas temperature (bottom) for an initial 1 M⊙ collapsing cloud core from the 1D simulations are shown at a time step after the second core formation. The yellow line indicates an inner radius of 10−4 au and the dashed red line indicates an inner radius of 10−2 au.
In the text

	Fig. B.1. Polar-angle averaged radial profiles of the density (top), velocity (middle), and gas temperature (bottom) for an initial 1 M⊙ collapsing cloud core at an initial temperature T0 of 10 K are shown at a time step after the second core formation. The different lines indicate the results using various grid resolutions in the 2D simulations.
In the text

Fig. B.2. Line integral convolution visualisation of the second core formed from the collapse of a 1 M⊙ cloud core at an initial temperature T0 of 10 K and an outer radius of 3000 au. Top to bottom panels: results from runs using different resolutions in an increasing order as indicated in the legends. The entropy behaviour is shown at a snapshot when the central densities are roughly similar. For the results presented in this paper, we use the highest resolution, which allows us to resolve the convective eddies within the second core.

In the text

	Fig. B.3. Radial profiles of the density (top), gas temperature (middle), and entropy (bottom) for an initial 1 M⊙ collapsing cloud core at an initial temperature T0 of 10 K and an outer radius of 3000 au are shown. The yellow and red lines indicate the results using two different grid resolutions in the 1D simulations.
In the text

	Fig. C.1. Polar-angle averaged radial entropy profiles are shown at four different time snapshots as a 1 M⊙ cloud core transitions through the formation and evolution of the first and second hydrostatic cores.
In the text