A grid of onedimensional lowmass star formation collapse models^{⋆}
Centre for Star and Planet Formation, Niels Bohr Institute and Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5–7, 1350 Copenhagen K, Denmark
email: neil.vaytet@nbi.ku.dk
Received: 26 January 2016
Accepted: 12 October 2016
Context. Numerical simulations of star formation are becoming ever more sophisticated, incorporating new physical processes in increasingly realistic setups. These models are being compared to the latest observations through stateoftheart synthetic renderings that trace the different chemical species present in the protostellar systems. The chemical evolution of the interstellar and protostellar matter is very topical, with more and more chemical databases and reaction solvers available online to the community.
Aims. The current study was developed to provide a database of relatively simple numerical simulations of protostellar collapse as a template library for observations of cores and very young protostars, and for researchers who wish to test their chemical modelling under dynamic astrophysical conditions. It was also designed to identify statistical trends that may appear when running many models of the formation of lowmass stars by varying the initial conditions.
Methods. A large set of 143 calculations of the gravitational collapse of an isolated sphere of gas with uniform temperature and a BonnorEbertlike density profile was undertaken using a 1D fully implicit Lagrangian radiation hydrodynamics code. The parameter space covered initial masses from 0.2 to 8 M_{⊙}, temperatures of 5−30 K, and radii 3000 ≤ R_{0} ≤ 30 000 AU.
Results. A spread due to differing initial conditions and optical depths, was found in the thermal evolutionary tracks of the runs. Within less than an order of magnitude, all first and second Larson cores had masses and radii essentially independent of the initial conditions. Radial profiles of the gas density, velocity, and temperature were found to vary much more outside of the first core than inside. The time elapsed between the formation of the first and second cores was found to strongly depend on the first core mass accretion rate, and no first core in our grid of models lived for longer than 2000 years before the onset of the second collapse.
Conclusions. The end product of a protostellar cloud collapse, the second Larson core, is at birth a canonical object with a mass and radius of about 3 M_{J} and 8 R_{J}, independent of its initial conditions. The evolution sequence which brings the gas to stellar densities can, however, proceed in a variety of scenarios, on different timescales or along different isentropes, but each story line can largely be predicted by the initial conditions. All the data from the simulations are publicly available.
Key words: stars: formation / stars: protostars / stars: lowmass / hydrodynamics / radiative transfer / gravitation
The figures and raw data for every simulation output can be found at this address: http://starformation.hpc.ku.dk/gridofprotostars. Copies of the outputs, as well as Table C.1, are also available in the form of static electronic tables at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/598/A116
© ESO, 2017
1. Introduction
The generally accepted view that stars form by gravitational condensation of diffuse interstellar matter is very old; however, it is only in recent decades that we have gained physical understanding of the process, thanks to ever more precise observations and complex numerical models (McKee & Ostriker 2007; André et al. 2014). Stars form within turbulent molecular clouds when density fluctuations are large enough to become gravitationally unstable and begin to collapse onto themselves. The details of this process remain largely unclear, as many different physical mechanisms operate over a very wide range of spatial scales. The progenitorial density enhancements (also known as parent clouds) are initially optically thin and contract isothermally, while the compression heating is radiated away. As the density of the gas inside the cloud rises, the efficiency of the radiative cooling drops until it can no longer counterbalance the compressive calefaction. The gas evolution becomes adiabatic, and a hydrostatic body, known as the first Larson core after the pioneering numerical studies of Larson (1969), is formed at the centre of the collapsing system. It accretes material from its surrounding envelope, and the sustained increase in mass, density, and temperature eventually triggers the dissociation of H_{2} molecules, when temperatures exceed ~2000 K. This leads to the second phase of collapse because of the endothermic nature of the dissociation process. The collapse ceases when most or all of the H_{2} molecules have been split, at which point a much more dense and compact hydrostatic core (Larson’s second core) is formed (Larson 1969; Masunaga & Inutsuka 2000; Vaytet et al. 2013a). The temperature inside the second core continues to rise until the ignition of nuclear reactions: the young star is born.
This multiscale, multiphysics problem is very difficult to tackle numerically; the parent cloud scales ~10^{4} astronomical units (AU), while the protostar measures only 10 AU at birth. An intricate interplay operates between largescale environmental factors, which regulate the supply of mass, angular momentum, and magnetic flux (Li et al. 2015), and smallscale processes close to the star, which control the evolution and dynamics of protostellar systems (Hennebelle & Charbonnel 2013). Numerical models are becoming more sophisticated every year, incorporating new physical processes (Commerçon et al. 2010; Federrath et al. 2014; Bate et al. 2014; Tomida et al. 2015; González et al. 2015; Masson et al. 2016; Hopkins & Lee 2016; Wurster et al. 2016) in increasingly realistic setups (Krumholz et al. 2012; Seifried et al. 2012; Li et al. 2014; Nordlund et al. 2014; Padoan et al. 2014; Federrath 2015; Matsumoto et al. 2015)^{1}. The comparison between these models and the latest observations has also reached new heights in recent years, using radiative transfer tools and observatory simulators to produce astonishing synthetic maps and spectra (Commerçon et al. 2012a,b; Frimann et al. 2016a,b; Visser et al. 2015; Kataoka et al. 2015; Seifried et al. 2016). This analysis is heavily focused on the emission from the dust and the different chemical species present in the protostellar systems, and determining which structures they trace (discs, outflows, protostar). While the chemical abundances in these studies are derived chiefly from the dust and/or gas temperatures, much effort is currently devoted to including timedependent nonequilibrium chemical evolution in such simulations (Grassi et al. 2014, 2016; Kim et al. 2014; Hincelin et al. 2016; Dzyurkevich et al. 2016). Comprehensive databases of interstellar chemical reaction rates, such as the UMIST database (McElroy et al. 2013) or kida (Wakelam et al. 2015) for instance, are now also publicly available to the community.
In this context, the current study was developed to provide a database of relatively simple numerical simulations of protostellar collapse for researchers who wish to test their chemical modelling under dynamic astrophysical conditions. Instead of producing a collection of large complex turbulent 3D simulations, which have a prohibitively high computational cost and can be very cumbersome to analyse and interpret, the aim was to use a simple spherically symmetric setup and varying the initial parameters. The study extends the simulation set in Vaytet et al. (2013a) with 143 new calculations of the gravitational collapse of an isolated sphere of gas with uniform temperature and a BonnorEbertlike profile. The calculations solve the equations of radiation hydrodynamics with self gravity, using a nonideal gas equation of state and an uptodate table of dust and gas opacities. The code uses a Lagrangian mesh which follows the evolution of each piece of fluid throughout the entire simulation, making it ideal for tracing the chemical evolution of gas elements. This paper describes the numerical method used and the initial parameters for the different runs, and also reports on statistical trends that emerge from running a large set of collapse models. The masses and radii of the first and second Larson cores, among others, are given and the connection between the different initial parameters and final core properties are investigated.
2. Numerical method and initial conditions
The SINERGHY1D code was used to solve the equations of radiation hydrodynamics; it is a 1D fully implicit Godunov Lagrangian code, described in Vaytet et al. (2012, 2013a). It uses the M_{1} closure to model the radiative transfer (Levermore 1984; Dubroca & Feugeas 1999) and the grid comprises 4096 cells (see Appendix A). The mesh is fully Lagrangian, no regridding scheme is employed in this work, contrary to what has been done in some of our previous studies. To have a higher resolution towards the centre of the mesh, the cell size is increased progressively towards larger radii. The size of cell i + 1 is defined by dr_{i + 1} = (1 + α) dr_{i}, where we have chosen α = 8 × 10^{4}. The integration timestep was chosen so that the variations from one time step to the next of each primary hydrodynamic (mesh position, velocity and total energy) and radiative (energy density and flux) variable was limited to a maximum of 10%. This ensured rapid convergence of the NewtonRaphson iterations in our timeimplicit scheme.
The code incorporates the gas equation of state of Saumon et al. (1995), and its extension to low densities (see Vaytet et al. 2013a), for a mixture of hydrogen and helium, with a helium mass concentration of 0.27. The interstellar dust and gas opacities were reused from Vaytet et al. (2013a). These comprise, at low temperatures (below 1500 K), the opacities of Semenov et al. (2003) as the absorption coefficients are dominated by the one percent in mass of dust grains. For temperatures between ~1500−3200 K, the dust grains are rapidly destroyed and molecular gas opacities prevail. These were obtained from the model of Ferguson et al. (2005). Finally, at temperatures above ~3200 K, there are no more molecules and the atomic opacities from the OP project (Badnell et al. 2005) were employed. We also made use of the opacity averaging scheme presented in Vaytet et al. (2013b) to transition between Planck and Rosseland opacity in regions of different optical thickness.
Fig. 1 Model parameter space for 6 different radii: R_{0} = 3000,5000,7500,12 000,20 000,30 000 AU. The hatched area represents clouds that are stable against gravitational collapse, while the colours code for different values of ϵ for clouds in the unstable zone. The ϵ = 1 separating line between hatched and coloured areas is defined in Eq. (5). Each run number links to the corresponding publicly available dataset. (Interactive image with hyperlinks available in the PDF version.) 

Open with DEXTER 
The initial setup for the dense cloud collapse was slightly modified compared to Vaytet et al. (2013a), and somewhat similar to Tomida et al. (2013). A BonnorEbert (Bonnor 1956; Ebert 1955) like density profile was adopted, with a density contrast between the centre and the edge of the sphere ρ_{c}/ρ_{0} = 14.1, corresponding to a dimensionless radius of ξ = 6.45. To fully define the gas density, fixed values for the sphere’s temperature T_{0} and radius R_{0} were chosen. The density at the centre of the cloud for a critical BonnorEbert (BE) sphere is then found via (1)where k_{B} is Boltzmann’s constant, μ represents the mean atomic weight (= 2.31), m_{H} the hydrogen atom mass, and G the gravitational constant. This is used to scale the dimensionless BE density profile to physical densities. Integrating the mass inside the cloud gives the same result as the wellknown critical BE mass (2)where c_{0} is the isothermal sound speed of the gas. Finally, to obtain a given initial cloud mass M_{0}, we divide the density everywhere in the cloud by the ratio of BE mass to cloud mass ϵ = M_{BE}/M_{0}, which determines whether a cloud is stable (ϵ> 1) or unstable (ϵ< 1) to gravitational collapse. Combining Eqs. (1) and (2) with the definition of ϵ, we see that the central density is independent of the temperature, only depending on the average density in the cloud, and is given by (3)The cloud’s freefall time, defined as (4)gives a good estimate of the time needed to form the first Larson core. The formation of the second Larson core will follow almost “instantaneously” (on a timescale of ~100−1000 yr). The radiation temperature is initially in equilibrium with the gas temperature and the radiative flux is set to zero everywhere. The boundary conditions are reflexive at the centre of the grid (r = 0). At the outer edge of the sphere, the hydrodynamical variables have a zerogradient boundary condition, while the radiative variables are fixed to their initial values inside the ghost cells. The equations of radiative transfer were integrated over all frequencies (grey approximation) since including frequency dependence only yields small differences for a much increased computational cost (Vaytet et al. 2012, 2013a).
The parameter space was chosen to cover a large range of ϵ between 0 and 1, with different initial cloud masses, temperatures and radii. The cloud masses range from 0.2 M_{⊙} to 8 M_{⊙}, the temperatures cover 5−30 K, while six separate radii were chosen between 3000 and 30 000 AU. The parameter space is illustrated in Fig. 1, and all the run parameters are listed in Table C.1. Only gravitationally unstable clouds (ϵ< 1) were selected. The dividing line between stable and unstable clouds, for a given radius is defined by (5)where C_{BE} ≃ 1.18 is a constant. Each simulation was stopped shortly after the formation of the second Larson core. We were not able to follow the evolution of the system for very long beyond this point, even with our timeimplicit scheme. The integration timesteps are limited to ensure small enough relative variations (<10%) of the hydrodynamic and radiative variables from one time step to the next. Having a Lagrangian mesh with two more or less stationary accretion shocks implies that pieces of fluid must travel through the shocks, and as they do, they suffer large increases in density, temperature and radiative flux. This has the consequence of keeping the integration timestep very low (especially after the second core has formed, as it has very high accretion velocities), making any lengthy integrations into the main accretion phase very costly. Other authors (Foster & Chevalier 1993; Masunaga & Inutsuka 2000, for instance) have commonly switched to a Eulerian scheme once the core is formed, to follow the system’s evolution on much longer timescales. Rewriting an Eulerian version of our code to simulate the main accretion phase is deferred to future work.
3. Results
For the details on the evolution of a typical run, which forms a first and second Larson core, the reader is refered to the previous paper by Vaytet et al. (2013a) which provides an extensive discussion on the evolution and structure of spherically symmetric collapsing clouds. In the present work, we shall limit ourselves to studying trends in the results from the multiple simulations performed, attempting to extract statistically meaningful relationships between initial parameters and final core properties.
3.1. Central densities, temperatures and entropies
Figure 2a shows the thermal evolution (temperature as a function of density) for the central cell in the computational grid. All the tracks follow a similar evolution. They first go through an initial isothermal phase, where the low optical depth allows compression heating to be radiated away. As the density rises, the efficiency of the radiative cooling drops and the system begins to heat up, following an adiabat with constant entropy (see panel b). This adiabatic phase is refered to as the formation of the first Larson core.
Fig. 2 Top row: evolution of the central temperature (a) and entropy (b) as a function of density for all 143 runs. The lines have been colourcoded according to the cloud’s initial temperature T_{0} (see legend in panel b)). Bottom row: c) evolution of the central temperature as a function of density for runs 30, 3, and 140 (see legend). d) Optical depth of the systems as a function of radius for runs 30, 3, and 140 (same colours as in panel c)), when the central density has just reached log (ρ_{c}) = −12.75 g cm^{3} (solid), and when the central temperature has reached log (T_{c}) = 1.2 K (dashed). 

Open with DEXTER 
Fig. 3 Evolution track in the densitytemperature plane for 10 different pieces of fluid in runs a) 69, b) 86, c) 71, and d) 115. In all panels, 10 cells were selected and followed throughout the simulation. The number of the cell in the mesh is given for each track in the bottomright inset, along with the cell’s radius at the start of the simulation (“s”) and at the end (“e”). While the cell numbers are the same in all panels, the cell radii change. The colour codes for the time in kyr; the colourbar is different in all panels, but they all display the final 700 yr of evolution. The pink dots show the position of each piece of fluid at the end of the simulation. The dashed grey line traces the final (ρ,T) profile in the entire grid, which connects the pink dots. The inset in the top left corners gives a reminder of the parameters of each run. 

Open with DEXTER 
The 143 runs show a wide spread in the thermal tracks they follow during this first adiabatic phase; its origin is twofold. First, as shown from the colour coding in the top row of Fig. 2, a vertical spread is seen to be the direct consequence of the different initial temperatures for the various runs. In the first adiabatic phase, the effective ratio of specific heats of the gas γ_{eff} changes from 1 (isothermal collapse) to ~7/5. This sets the slope of the tracks, which is approximately the same for all runs. This implies that a system starting with a higher initial temperature T_{0}, as it becomes optically thick, will follow an adiabatic track that is shifted upwards in the vertical direction. Moreover, the Rosseland mean opacity of the dust in the low temperature region behaves like κ_{R} ∝ T^{2}. The runs with high initial temperatures will therefore have higher optical depths even when the density is the same, and will consequently become adiabatic earlier, attaining higher entropies (as seen in Fig. 2b). In addition, a second spread is also observed in runs with identical T_{0}, this time in the horizontal direction. To determine the origin of this spread, we select three runs (30, 3, and 140) with T_{0} = 10 K but different initial masses and radii. The thermal evolutions of these three runs is displayed in Fig. 2c. Runs 30 and 3 begin along the same isothermal path, but start heating up at different central densities, following separate tracks in the first adiabatic phase. This first differentiation appears to largely determine the subsequent evolution of the system, as the runs never rejoin the same tracks during the rest of the simulations. The cloud begins to heat up once the radiative cooling is no longer efficient compared to the compressional heating (see the discussions in Masunaga et al. 1998; Masunaga & Inutsuka 1999). Figure 2d shows the optical depth in the systems at two different times: when the central density has just reached log (ρ_{c}) = −12.75 g cm^{3} (solid lines), and when the central temperature hits 16 K (dashed lines). It is clear that, for the chosen central density (when the spread in temperature is noticeable), the higher the initial mass, the higher the optical depth, and hence the more the cloud has heated up. The dashed curves illustrate that the runs show the same temperature rise for approximately the same optical depth. Comparing to run 140 supports further the statement that the evolution is governed by the optical depth. Run 140 has the same mass as run 30, but a radius 10 times larger, greatly diminishing its optical depth. In fact, it turns out to have an optical depth very similar to run 3 (see panel d), and both runs have matching evolutions. In our case, isothermality breaks down around an optical depth of unity or slightly above, but we note that, as pointed out by Masunaga & Inutsuka (1999), this does not by any means constitute a general rule. The determining factor is the balance between radiative cooling and compression heating, which is not always ruptured for an optical depth of unity. To estimate the density at which the cloud begins to depart from isothermality in its early evolution, we refer the reader to the exhaustive discussion in Masunaga et al. (1998).
The first spread in thermal tracks at the start of the first adiabatic phase appear to also set the entropy level at the centre of the system for the remainder of its evolution. Indeed, Fig. 2b shows that after a first phase of decrease during the isothermal collapse, the central entropy remains largely constant for ρ_{c}> 10^{12} g cm^{3} (it should be noted that while the entropy at the centre of the system decreases, the entropy of the system as a whole increases). The spread in entropy at the centre of the protostar, at the end of the simulation, is quite significant compared to the average differential the central piece of fluid suffers during the entire simulation (~25% of the average total differential). This implies that the different protostellar seeds formed in our simulation will have different early evolutions, as small differences in entropy can cause large variations in luminosity (see for example Stahler & Palla 2005).
3.2. Evolutionary tracks for single pieces of fluid
As this work was motivated by the evolution of chemical species in circumstellar matter going through various stages of protostellar evolution, we now turn to following the thermal evolution of other parts of the collapsing system, instead of just the centre. As many chemical evolution codes have no spatial dimensions and simply follow the chemical reactions and abundances in time, they could very easily use the thermal properties of a piece of fluid inside the mesh as a function of time as an input. And as our mesh is purely Langrangian, following a piece of fluid means selecting a single cell in the grid and tracing its time evolution. We plot in Fig. 3a the evolution in the (ρ,T) plane of 10 different pieces of fluid in a “standard” run 69. It has an initial cloud mass of 2 M_{⊙} , temperature of 10 K, and an intermediate radius of 7500 AU. The 10 pieces of fluids were selected to illustrate different behaviours, and the cell number of each piece of fluid can be found in the lowerright inset in Fig. 3a. We see that the cells close to the centre of the system (tracks 1 and 2) follow very closely the types of evolution shown in the previous section. Their temperatures increase smoothly as they are already beyond the first and second core accretion shocks when these are formed, and thus do not undergo rapid changes in density/temperature. Track 3 illustrates a piece of fluid that has just passed through the second core shock, showing a sharp rise in both density and temperature at late times. Tracks 4 and 5 both depict fluid parcels that have passed through the first core accretion shock, but have done so at different times. In this diagram, a shock is identified by a discontinuity in density. Indeed, purely hydrodynamic shocks would also exhibit a temperature discontinuity, but strong supercritical radiative shocks (such as the first core accretion shock) have equal upstream and downstream temperatures. A shock is then identified because it implies strong compression, but not necessarily a jump in temperature. Track 4 shows a sharp (horizontal) density jump around a density of 10^{12} g cm^{3}, while the discontinuity in track 5 lies closer to 10^{13} g cm^{3}. As we will see further down (Fig. 5), the first core has a size of about 2 AU at the time of formation, which corresponds to the density of the jump in track 4 (this can be seen in the online data), and subsequently grows in size as it accretes material from the surrounding envelope. By the time track 5 passes through the shock, the core’s radius is located around 8 AU which now corresponds to a density of 10^{13} g cm^{3}. In these two cases, the rises in temperature just upstream of the shock are due to radiative heating from the first core, which heats up the outer envelope. This is further illustrated by the final (ρ,T) profile at the end of the simulation represented by the grey dashed line. Fluid parcel 6 has just gone through the accretion shock, whereas parcels 7 and 8 only get heated via radiation. It is important to note here, especially for chemical evolutions, that the inner parts of the system evolve much more rapidly that the outer envelope. Indeed, we can see that over the timeframe of the simulation, the inner parts reach protostellar densities and temperatures, while the outer gas stays below 100 K. This will inevitably lead to very contrasting chemical evolutions, depending on which piece of fluid is considered.
The evolution will also differ when different initial parameters are used. Figure 3b shows the evolution of pieces of fluid located at the same positions in the mesh, but for a simulation with a much smaller freefall time (run 86 has the same initial radius and temperature than run 69, but four times the initial mass). In this case, the most striking difference is that only the very central cells reach protostellar densities and temperatures. Tracks 2 and 3 have not yet reached the second core, while tracks 5 to 7 are still all outside of the first core at the end of the simulation. The collapsing gas goes through the entire first adiabatic phase, second collapse and second adiabatic phase much faster than in run 69 (in run 86, most of the curve is coloured with shades of red and orange, meaning that this evolution takes place within the very last quarter of the time colour bar, for t ≥ 18 kyr). This explains why the cells in the outer envelope remain at a low temperatures and densities in this run; the gas simply has not had enough time to travel inwards and heat up before the second core was formed and the simulation was stopped. As seen in panel a, it does take some time for the temperature of fluid parcels 7 and 8 to rise. Panel c illustrates the same evolution as in panel a but for run 71 that has an initial temperature of 20 K. While the outer cells, which sit at a temperature twice as high as in run 69, may see different chemical reactions dominating, the evolution of the rest of the system is very similar to run 69. We note, however, that the simulation end time has changed significantly in run 71 compared to run 69 due to the additional pressure support in the hotter cloud at the start of the simulation, showing that for marginally stable BE clouds, the freefall time is only an order of magnitude estimate for the collapse time, as it does not consider pressure support in the cloud (runs 69 and 71 have identical freefall times). Finally, panel d shows the fluid parcels’ evolution in a case where the initial density is a factor of almost 20 lower. The freefall time is about four times as long, and it will take the gas much longer to reach the high end of the (ρ,T) space. As chemical reactions are highly timedependent, this can also greatly modify the abundances of the species in the gas. Barring the different timescales, the evolution of the different fluid parcels are remarkably similar to those in run 69. Chemical abundances are of prime importance in protostellar evolution, as they regulate the gas hydrodynamics (through the heat capacity, line cooling, and chemical energy), the radiation (through the gas and dust opacities) and the magnetic field (through resistivities).
3.3. Radial profiles
Fig. 4 Radial profiles of the gas density (a), velocity (b), temperature (c), and luminosity (d) for all 143 runs. The different colours represent different initial temperatures T_{0} (see legend). The vertical dotted lines mark the positions of the average first and second core radii. 

Open with DEXTER 
We now take a look at the radial profiles of various quantities and how they vary in the different runs. Figure 4 shows the density, velocity, temperature and luminosity profiles for all 143 runs. As in Fig. 2, the different colours represent different initial cloud temperatures. We first notice that, considering the range of initial parameters used, all the profiles are remarkably similar, especially inside the first core radius. The average first and second core radii are marked by the vertical grey dotted lines, and we see very little spread around these values. This result was already found in Vaytet et al. (2013a), and the present study comfirms the statement accross a much larger variety of simulations. There is however some spread in the different profiles in the protostellar envelope, outside of the first core accretion shock. The initial densities and temperatures vary among the different runs, and it thus makes sense to see a spread in the outer regions of the system. This however shows that once a hydrostatic body is formed, the system “forgets” its initial conditions and all the profiles lie on top of each other. Finally, we find a spike in luminosity around 1 AU in all runs, which actually lies inside the first core and occurs at the dust destruction front, and not at the first core accretion shock. A small spike is also visible just at the second core accretion shock.
3.4. Properties of the first and second Larson cores
Fig. 5 First and second core mass (a) and radius (b) for all 143 runs. Each mass and radius has a starting and an end point. The starting point represents the quantity when the core has just formed, while the end point is at the end of the simulation. A small arrow in the middle of each segment points in the direction from start to end. The vertical alternating grey and white bands indicate changes in initial radius (see panel b) for values of R_{0}). 

Open with DEXTER 
Table C.1 lists the properties of the first and second Larson cores (mass, radius, etc.) at the end of each simulation. Figure 5 is a more graphical representation of masses (a) and radii (b) for the first (circles) and second (squares) cores. It brings further information as it shows the cores’ masses and radii not only at the end of the simulation, but also just when the accretion shock at the core borders have formed.
Looking at Fig. 5a, we see that, as expected, all cores are accreting material and growing in mass over time. While the first cores show a significant spread in mass at the time of formation (standard deviation of σ = 40%), their masses are much more similar at the end of the simulations (σ = 13%). In the case of the second cores, we see that all runs show very homogeneous masses in the range 2−4 × 10^{3}M_{⊙}. The jumps in final first core mass from one run to the next coincide with changes in initial radius R_{0}, illustrated by the alternating grey and white background colour bands. Panel b reveals that the simulations produce first cores with very comparable radii. All of them increase with time, starting from ~2 AU at core formation, and reaching ~10 AU at the end of the runs. This does contradict some earlier results who describe the first core as contracting with time (Masunaga et al. 1998; Tomida et al. 2013). To understand this further, we plot the time evolution of the first core radii in Fig. 6a.
The different colours represent different initial cloud temperatures T_{0} (see legend). We can see that all the runs actually display an initial phase of contraction which lasts less than about 10 yr, followed by an expansion phase. In panel b, we plot the first core radii as a function of central density, alongside the analytical predictions of Masunaga et al. (1998). It should be noted that their expressions have been rescaled to match the runs with T_{0} = 5 K; the scaling constant of 5.3 AU in their paper has been changed to 2.5 AU ^{2}. The Masunaga et al. laws actually reproduce the slope of the initial contraction phase very well, and runs with the same initial temperature appear to evolve along the same powerlaw early on. The spread in initial core radius for varying initial cloud temperature goes in the right direction, but is larger than the one observed in the simulations. While the initial contraction phase is in good agreement with previous studies, to understand the subsequent evolution of the first cores further, we isolate run 69, which exhibits a behaviour typical of the majority of the models. This run has been highlighted in panels a and b with a thick black solid curve.
In Fig. 6c, we plot the different pressures on each side of the shock, which enter the momentum conservation RankineHugoniot relation for a radiating fluid (see Mihalas & Mihalas 1984, Eq. (104.9)); namely the gas internal pressure P^{gas} (blues), the ram pressure P^{ram} = ρu^{2} (red/orange), and the radiation pressure P^{rad} (greens). The subscripts 1 and 2 represent quantities upstream and downstream of the shock, respectively. We first notice that the radiation pressure on either side of the shock is orders of magitude lower than the other quantities, and does not play a significant role in the momentum and energy exchanges. An accretion shock forms when the gas pressure cancels or exceeds the ram pressure of the infalling material. We can see that the postshock gas pressure (bold blue) and the preshock ram pressure (bold red) are comparable at early times (t< 10 yr). In fact, is slightly higher than , which leads to an early contraction of the first core. For convenience, we have plotted in panel c the difference between the downstream gas pressure and the upstream ram pressure (we note from panel d that is always at least a factor of 3 below , and will therefore have virtually no impact on the analysis). Eventually, as the gas pressure inside the core increases, it overcomes about 10 yr after core formation. This is followed by a phase where is noticeably higher than , causing the core to expand. The expansion is halted when the gas and ram pressure become comparable once again, and the radius stabilizes to ~ 6 AU (small bumps in this and other runs show oscillations around hydrostatic equilibrium). While the radiative pressures do increase during this expansion phase, they remain much smaller than the hydrodynamic quantities. The radiative pressure is only important in very strong radiative shocks, but the radiative flux, which enters the energy conservation RankineHugoniot condition, is important in all radiative shocks (see Mihalas & Mihalas 1984, Sect. 104 and Eq. (104.10)). Indeed, we also show for relative comparison in Fig. 6d the radiative flux (normalised with the local gas velocity for dimensional homogeneity) upstream (black) and downstream (grey) of the accretion shock. The radiative fluxes are clearly seen to increase as the shock moves outwards, illustrating the fact that radiative effects become important when the shock reaches larger radii (the spikes in the curves correspond to sign inversions). As the shock moves outwards, the temperature of the envelope falls (see the temperature profiles in Fig. 4c). This means that the dust mean opacity, which scales approximately with T^{2}, also declines, enabling radiation to escape once the opacity just outside of the accretion shock is small enough.
Another way of looking at the radiative efficiency of the accretion shock is to examine the shock compression ratio, ℛ_{sim} = ρ_{2}/ρ_{1}. We compare in panel e the compression ratio in the simulation (black line) to one predicted from the pre and postshock pressures for a purely hydrodynamical gas ℛ_{hyd} (red). Because of radiative effects, we expect ℛ_{sim} ≥ ℛ_{hyd}. At early times, ℛ_{sim} is comparable to its hydrodynamic counterpart. It appears that the shock moves outwards until it becomes radiatively efficient, where the compression ratio begins to climb above 5, and the core expansion stalls. Towards the end of the simulation, a second period of expansion is observed. We believe this is due to a combination of two effects. The first comes from a sustained increase in temperature, and hence pressure, inside the core (as it continues to accrete material) which implies that (bold cyan) rises once again above (bold red), as seen in panel c around a time of 200−300 yr. The evaporation of dust in the centre of the system where temperatures exceed 1500 K could also contribute, as it allows radiation from the hot central region to raise the temperature in the first core. In addition, the first core heats up the envelope (see for instance track 7 in Fig. 3a), which leads to a new rise in opacity and optical depth of the envelope, reducing the ability of the accretion shock to radiate accreted energy; the fall in compression ratio at late times attests to that fact.
Fig. 6 a) First core radius as a function of time, where t = 0 corresponds to the time when the first core accretion shock is formed. The different colours represent different initial cloud temperatures T_{0} (see legend). The thick black line is for run 69. b) First core radius as a function of central density (solid lines), alongside the (rescaled) analytical predictions of Masunaga et al. (1998, dashed lines). The colour coding is the same as in panel a). c) Difference between the downstream gas pressure and the upstream ram pressure. d) Downstream gas pressure (blue), upstream ram pressure (red) and radiation pressures (greens) just upstream (subscript 1) and downstream (subscript 2) of the first core accretion shock, in run 69. The upstream radiative fluxes, scaled with the local gas velocity for dimensional coherence, are also plotted for relative comparison in the total energy budget (black). e) First core accretion shock compression ratio for run 69 measured in the simulation (black) and predicted for a nonradiative shock (red). 

Open with DEXTER 
Returning to Fig. 5, we see that the second core radii are profoundly alike across the board, with an average final size of 4 × 10^{3} AU. We however note that while the first core radii increase with time, the second core radii decrease. Larson (1969), Stahler et al. (1986) and Tomida et al. (2013), to name a few, have all reported second Larson cores which grow in size over time. We believe that our present results are not inconsistent with these other works, as we have stopped the simulations at very early times, while the second core is still settling. It is very possible that the second cores behave like the first cores, with an initial contraction phase, and a subsequent increase in size due to heating and/or mass accumulation (which also depletes the accreting envelope).
3.5. First core formation time and freefall times
In this section, we take a closer look at the collapse time of the cloud, that is the time it takes to form a first Larson core, which we shall call t_{core1}. In comparison, the second core forms almost instantaneously (as we will see in the next subsection) after the formation of the first core, and the full simulation end time is essentially identical to the first core formation time. We stated in Sect. 2 that the freefall time gives a good estimate of the time it takes for the cloud to collapse and form a protostar. We have however found this not to hold in many cases, especially in clouds with a significant amount of thermal support. In Fig. 7, we compare the first core formation times to the freefall times computed using the cloud’s inital central density ρ_{c} (t_{ff1}, blue squares). While many points lie on the ideal dashed line, there are also a significant number of symbols close to the critical BEmass at the highϵ end which are at a large distance from it, with up to an order of magnitude in disagreement. An alternate definition of the free fall time, based on the cloud’s initial mass and radius (t_{ff2}, green triangles), shows a similar trend, but lies below a ratio of 1 for most of the points, which can be easily understood given the proportionality of the two (cf. Eq. (3)). The definition of the freefall time only considers gravity and ignores pressure forces from the gas, and it is thus not too surprising that in marginally unstable BE spheres, it does not perfectly predict the collapse time.
Fig. 7 Ratio of the time it takes in each run until the first core is formed (t_{core1}) to various collapse time estimates, as a function of ϵ. The blue squares represent the freefall time computed using the cloud’s inital central density ρ_{c} (t_{ff1}), while the green triangles show the results from an alternate definition of the freefall time based on the cloud’s average density (t_{ff2}). The black stars represent the cloud’s initial sound crossing time (t_{sound}). Finally, the best estimate is obtained with a modified version of t_{ff1} which takes into account how unstable the cloud initially is to gravitational collapse by introducing the ϵ parameter (t_{ff3}). 

Open with DEXTER 
A timescale which does consider pressure forces is the sound crossing time of the cloud t_{sound} = R_{0}/c_{0}, and this is represented by the black stars in Fig. 7. This timescale also fails to correctly estimate the first core formation time. However, it is obvious that the first core formation time depends on the ϵ parameter. Indeed, looking at the t_{ff1} estimates, clouds with low epsilon (i.e. very unstable initial conditions) lie right on the dashed line, when gravity is the dominant force. When clouds have an ϵ close to unity, where pressure support and gravity forces are comparable, it takes longer for the collapse to initiate. We found that correcting the t_{ff1} to now read (6)gives a very good estimate of the first core formation time; this is represented by the red circles in Fig. 7. The freefall timescale is dominated by the first stages of the contraction, as the collapse time becomes shorter and shorter the higher the density gets. It therefore makes sense that, at early times, the effective gravitational pull is given by that which is unbalanced by the pressure gradient. One would then expect the correct density to enter the formula for the freefall time to be the differential density compared to the critical BE density: ρ_{c}−ρ_{c,crit} = ρ_{c}−ρ_{c}ϵ, which is exactly what we find in Eq. (6). We can also see that as ϵ goes to 0, we recover the original definition t_{ff1}, while as ϵ reaches unity, t_{ff3} goes to infinity and the cloud does not collapse.
3.6. First core lifetimes
Fig. 8 a) First core lifetime τ_{fc} as a function of ϵ measured in the simulations (circles). The dashed and thin solid lines represent our analytical predictions (see text). b) First core lifetime as a function of mass accretion rate for all runs. The black dashed line portrays a simple estimate from dimensional arguments (see text). c)–i) Ratio of the mass accretion rate measured in the simulation to various predictions Ṁ_{pred} from numerical and similarity solutions, as a function of ϵ. c) The hearts represent the similarity solution of Larson (1969) and Penston (1969), the stars are the solution of Shu (1977), while the triangles are from the results of Hunter (1977). In panels d) to f), the Ṁ_{pred} estimate is obtained from dividing the initial cloud mass M_{0} by the freefall times t_{ff1}, t_{ff2}, and t_{ff3}. g) The hearts represent a modified version of Penston’s formula, which is equivalent, barring a scaling factor, to the prediction using the cloud sound crossing time (diamonds). h) The circles represent the first best fitting formula Ṁ_{A} we have used to predict the lifetimes. i) The squares show the second mass accretion rate estimate Ṁ_{B} we have selected for our analysis. The colours in all panels code for different initial temperatures (see legend in panel a)). 

Open with DEXTER 
This section focuses on the time elapsed between the formation of the first and second cores, which we will refer to as the “first core lifetime”. Once the second core is formed, the observable properties of the system such as spectral energy distributions dramatically change, due to the high temperatures inside the second core (see for example Masunaga & Inutsuka 2000)^{3}. The first core lifetimes for all the runs are displayed in Fig. 8. The top panel a presents the lifetimes as a function of the initial mass to BonnorEbert mass ratio ϵ. For ϵ ≳ 0.4, the lifetimes appear to depend only weakly on ϵ, showing a spread between 100 and 1000 yr. Below, the lifetimes drop with ϵ.
If one assumes that the first core is in quasihydrostatic equilibrium, one can see that, if it is continuously accreting material from its outer envelope, its temperature will also steadily increase. It is also safe to assume that the higher the mass accretion rate, the faster it will heat up and enter the second phase of collapse. This therefore predicts a correlation between first core lifetime and mass accretion rate, and we show this in Fig. 8b, where the tight correlation is glaring (the value for the accretion rate is taken as the average rate over the lifetime of the first core). Assuming a uniform density and temperature, the equation of hydrostatic equilibrium for a sphere is (7)Approximating dP/ dr = (P_{0}−P_{c}) /R where P_{c} is the central core pressure, P_{0} (≪ P_{c}) is the outside pressure, and R the core radius, we find (8)which we finally convert to temperature to obtain (9)If we assume the radius of the first core remains approximately constant, there is a critical mass M_{crit} that the first core must reach, so that T_{c} = 2000 K. This is further justified by the fact that all first cores appear to reach very similar masses at the end of the simulations in Fig. 5a. In fact, by setting R = 5 AU in Eq. (9), we obtain M_{crit} = 0.04 M_{⊙}, which is fully consistent with Fig. 5a. We can now define the first core lifetime as being the time it takes for the core to reach M_{crit}, having started at an initial mass M_{init} at the time of formation. One can go even further by neglecting M_{init} altogether as it is on average 10 times lower than M_{crit}, allowing us to write the very simple relation (10)where is the first core mass accretion rate. This relation (for M_{crit} = 0.04 M_{⊙}) is plotted with a dashed black line in Fig. 8b, alongside the measured first core lifetimes and mass accretion rates from the simulations. Although not perfect, the law is a decent fit to the measured lifetimes, considering the simplicity of the argument we have used.
To try and estimate the first core lifetimes from the initial conditions, we now only need to predict the first core mass accretion rates, as they are fully correlated to the lifetimes. Many authors have studied numerically and analytically the problem of an isothermal sphere of gas, collapsing under its own gravity. As discussed in Foster & Chevalier (1993), several estimates of mass accretion rates were derived by Larson (1969), Penston (1969), Shu (1977) and Hunter (1977). They are expressed as a function of the cloud’s isothermal sound speed , where K is a constant equal to 47 in the LarsonPenston solution, 0.975 for Shu, and 36 for Hunter. These values were computed at the epoch of first core formation, and are not constant in time (see Hunter 1977, for more details), but are useful to obtain order of magnitude estimates. We show these predictions in Fig. 8c, where the mass accretion rates measured in the simulations are compared to the aforementioned predictions, as a function of ϵ. All three distributions have the same shape, and show a distinct departure from the measured mass accretion rate at low ϵ. We also note here that Shu’s prediction performs much worse that the other two, with ratios always higher than 10. A different way of looking at this problem is to estimate the mass accretion rate through a very simple dimensional analysis. We suppose that the sphere of gas collapses onto itself in a freefalling manner, meaning that all the mass will reach the centre of the system in one freefall time. The mass accretion rate is then simply the mass of the cloud divided by the freefall time (11)The predictions for the three different definitions of the freefall time from Sect. 3.5 are shown in Figs. 8d–f. While the dimensional estimate works better than the predictions in panel c, we can see that a significant amount of spread above and below a ratio of unity is still present, especially in the t_{ff3} case, showing that even a good estimate of the cloud collapse time does not provide a reliable measure of the mass accretion rate. The final timescale we used in Sect. 3.5 was the sound crossing time, and the mass accretion rate obtained from Ṁ_{sound} = M_{0}/t_{sound} is shown in Fig. 8g (diamonds). This time, the distribution is close to horizontal; it only lacks a scaling constant to bring it onto the dashed line. In fact, a similar expression can be found from the Penston solution, if we take into account the effect of the ϵ parameter. Indeed, the increasing ratio with decreasing ϵ panel c prompted us to attempt to correct the formulas by incorporating ϵ in the expressions. Simply dividing the Penston expression by ϵ yields (12)which is identical to Ṁ_{sound}, except for a factor of ~20. This expression also implies that for an infinite cloud mass with a given radius, ϵ → 0 and the mass accretion rate goes to infinity. We have ultimately decided to use a scaled version of Ṁ_{sound} (which is the same as scaling the modified Penston estimate) as it showed the best agreement with the measured mass accretion rates, with a factor of 6, (13)as shown in Fig. 8h. This constant probably accounts for the fact that the mass accretion rate onto the core is not constant in time, but begins slowly, and accelerates as the density increases in a runaway fashion during the early stages of the collapse. The value of 6, however, implies our dimensional analysis agrees within an order of magnitude with the numerical models. The prediction is solid for ϵ ≤ 0.3, above which it starts to depart from the ideal dashed line. A second attempt at accurately predicting the first core mass accretion rates was made by correcting the expressions in panels d–f with the ϵ parameter. Multiplying the initial mass M_{0} by ϵ is equivalent to using the BE mass in the numerator of the Ṁ_{pred} formula, and again using a scaling constant, we found (14)to give good results, as shown in Fig. 8i. This time the agreement with the measured mass accretion rates is excellent at high ϵ, and only starts to break down for ϵ< 0.2.
Returning to the top panel of Fig. 8, we can now express the mass accretion rate, and hence the first core lifetime, in terms of ϵ. This yields when using the Ṁ_{A} estimate in Eq. (13). This relation is represented in Fig. 8a by the dashed coloured lines for different values of T_{0}. The temperature spread in lifetimes for a given ϵ is well reproduced by the analytical function. At the highϵ end, the analytical description overestimates the lifetimes, as the Ṁ_{A} prediction worsens for the most stable runs (as seen in the right hand side of panel h). Using the Ṁ_{B} expression in Eq. (14) gives which is plotted in Fig. 8a with thin solid lines. Once again the range of temperatures is adequately covered, but this time the fit is much better towards high values of ϵ, where the ϵ^{1/2} slope follows the data much more closely. The two relations are equal for ϵ_{AB} = 9π^{2}/ (25 × 8 × 2.4) = 0.185, and they are a good match to the data in their respective ϵ ranges.
3.7. Further correlations
To extract further correlations between the different simulation initial parameters and final core properties from Table C.1, we computed for each column in the table the correlation of the data in that column to all the other columns. The correlation between columns x and y is defined by (15)where the covariance is given by (16)and (17)is the weighted mean. The subscript i denotes a particular run. The weights are used to minimise the influence of outlier data points, and the method used to compute the weights is given in detail in Appendix B. This correlation matrix r_{xy} is plotted in Fig. 9 (it should be noted that the arrays x and y actually represent the logarithm of the quantities in Table C.1). Within the initial parameters, some obvious strong correlations appear, such as (ρ_{c},t_{ff1  2}) or (ρ_{c},R_{0}). It is also selfevident that the simulation end time (or t_{ff3}) should correlate with the other freefall times. Some correlations also emerge between the run number and various parameters related to the initial radius, which comes from the manner in which we have ordered the runs from lower to higher radii.
Fig. 9 Correlation matrix for the parameters and core properties in Table C.1. The correlations of the type r_{xx} in the white diagonal have been ignored. For r_{xy}< 0.5, the values are given in the cells for indicative purposes, and have intentionally been made less visible to avoid cluttering the figure with too many insignificant values. The interested reader should be able to see the numbers by zooming in enough on the figure. Each square links to a correlation plot stored in the online database. (Interactive image with hyperlinks available in the PDF version.) 

Open with DEXTER 
Among the first and second core parameters, we also unsurprisingly find that accretion luminosity is correlated to the mass accretion rate, since the accretion luminosity is computed from the mass accretion rate. It also relates to the mass and radius of the cores. As in Sect. 3.6, we find once again that the first core lifetime τ_{fc} is strongly correlated to the mass accretion rate . More interestingly, it appears to also be correlated to the first core radius. As we have seen in Sect. 3.6, the first core grows in size with time for most runs, and it seems obvious that the longer the first core lives, the bigger it will get. This correlation is actually visible in the top right corner of Fig. 6a, where the final radius R_{1} increases for longer lived cores. The luminosity at the outer radius of the simulated system L_{out} appears quite solidly related to, among others, the first core mass, as result fully expected from fundamental stellar physics where stars are considered to be radiating spheres of gas (see for example Phillips 1999, chap. 1). Fitting a power law to the luminosities yielded the relation , which is close to the value of 3.5 commonly quoted in text books.
In the second core properties, the core radius naturally correlates well to the accretion and radiative luminosities (r> 0.85) as it enters the calculation of both (L^{acc} = GMṀ/R and L^{rad} = πR^{2}F), while the core mass has relatively strong relationships with all of the second core properties (r> 0.7). The accretion luminosity and the radiative luminosity are also correlated, even though they differ by many orders of magnitude; they are most probably related via the core radius. Strangely enough, the ϵ parameter does not seem to strongly correlate, on its own, with the first and second core properties. As it was found to be an important parameter in the preceding sections, this illustrates that using correlations alone are not enough to extract a complete physical understanding of the collapse of spherical clouds.
There are many possible relations to be interpreted from this data, and going through them all would be too lengthy for the present article. Each square in the figure (coloured or grey) links to a figure on the online server which displays the x quantity plotted against the y, providing the reader with the entire collection of plots.
4. The collapse of a uniform density cloud
Throughout past studies of the collapse of an isolated sphere of gas, two main types of initial conditions have been used. The first, as was employed in this work, is a supercritical BonnorEbert sphere. The second is simply a sphere with a uniform density. There are many arguments against how a uniform density sphere could form in nature, but a such a body was found to very rapidly adopt a BonnorEbert like profile once it started to collapse “forgetting” its initial state (see for instance Larson 1969; Masunaga et al. 1998), and the subsequent evolution was deemed not to depend on the initial density profile. Some authors, including Masunaga & Inutsuka (2000), however argued that the density profile does affect the protostellar evolution, mass accretion rates and luminosities, due to different dynamics. In this section, we rerun some of our simulations by adopting a uniform density as a starting point, and compare them to our standard runs.
Fig. 10 a) Final density profiles for BonnorEbert (dashed) and uniform density (solid) setups with initial temperatures of 5 K (blue) and 30 K (red). b) Same as in a) but for the gas velocity. c) Ratio of the gas thermal pressure to ram pressure as a function of radius. The arrows point to the absence of a first Larson core in the most unstable uniform density run. 

Open with DEXTER 
We selected intermediate runs 95 to 100, which all have an initial mass of 4 M_{⊙}, initial radius of 12 000 AU, and varying temperatures. We show the results of the uniform density simulations in Fig. 10 (solid lines), alongside the original BE runs (dashed), for initial temperatures of 5 and 30 K. We first note that the initial density profile does not significantly affect the final density and velocity profiles, as well as the sizes and masses of the protostellar cores, when T_{0} = 30 K. This was found to hold true for all initial temperatures above 5 K. The case of T_{0} = 5 K is however intriguing, as it appears to be devoid of a first core accretion shock. We observe an accretion shock at the second core border, but the flow is continuous throughout the rest of the infalling envelope. There are no discontinuities in either density or velocity, as indicated by the arrows in Fig. 10. The system has proceeded directly to the second collapse, before forming a hydrostatic body during the first adiabatic phase. This occurs when the gas at the centre of the system rapidly reaches the dissociation temperature of H_{2} molecules T_{H2} = 2000 K, due to a strong compressional heating. It can also be understood in terms of when a cloud is initially highly unstable (because it has very little initial thermal support with such a low initial temperature of 5 K), the strong ram pressure due to the high infalling velocity is always above the thermal pressure inside the core, and an accretion shock never forms. This is illustrated in Fig. 10c which displays the ratio of thermal pressure P_{them} to the ram pressure P_{ram} = ρu^{2} as a function of radius at the end of the simulations. For the uniform T_{0} = 5 K run, unlike the other three runs, the thermal pressure never counterbalances the ram pressure from the infalling gas in the region where the first core accretion shock typically forms (~1–10 AU). Contrastingly, all runs have a clear accretion shock that forms at the second Larson core border (3−5 × 10^{3} AU).
After having run many more models with a uniform density profile (~90 additional runs), only the most unstable (ϵ< 0.1) setups were found to lack a first Larson core. The other runs behaved very similarly to the BE models. This does not invalidate previous studies which were using a uniform initial density, as the clouds were never unstable enough to skip the first core stage. We believe that this peculiar behaviour is because with an initial uniform density, the freefall time is the same everywhere in the cloud. This implies that all the mass should reach the centre at the same time. The freefall time considers gravity to be the only important force, ignoring gas pressure forces, which are in fact still able form a first core, even when the freefall time is homogeneous. On the other hand, when gravity strongly dominates in very unstable clouds, the strong infalling motions do not permit the formation of a first hydrostatic core. The pressure inside the second Larson core is so high that the collapse is eventually halted, but the system’s density and velocity profiles strongly resemble that of a singular isothermal sphere.
5. Comparison with observations
The initial conditions of our simulations can be compared to the observations of Könyves et al. (2015) who performed an extensive survey of dense globules^{4} in the Aquila cloud complex. We selected from their sample all the globules that were marked as “prestellar” and “protostellar”, and plotted them in Fig. 11, where we have made use of the deconvolved radii in their database. In each panel, the squares and circles represent the observations, while the grey area deliniates the region in the parameter space covered by our simulations. The empty circles symbolise observed globules with ϵ > 1, which we predict will not be collapsing under the influence of gravity alone. Their masses are below the BonnorEbert critical value, and they will only form protostars if the external (thermal or kinetic) pressure increases their density, or if they are allowed to cool down until the internal pressure gives in to gravity. Most empty circles actually lie outside of our parameter space coverage, as we have, for obvious reasons, decided to only simulate systems with a mass higher than M_{BE}.
Fig. 11 Comparison between the radii, temperatures and masses of the prestellar cores in the Aquila molecular cloud and the parameter space covered by the initial conditions of our set of numerical simulations. In each panel, the empty circles represent observed cores with ϵ > 1, while the coloured squares are for cores with ϵ < 1. The colour of the squares code for their predicted first core lifetimes. The grey area deliniates the simulation parameter space: every corner or intersection of a dotted line with either another dotted line or the thick grey boundary represents one simulation. 

Open with DEXTER 
The colours of the squares represent the estimated lifetimes of their first cores before the second collapse takes place, obtained from our analytical estimates in Eqs. (17) and (20), where we have used for ϵ<ϵ_{AB} and elsewhere. As all the coloured globules have temperatures below 20 K, we estimate that almost all first cores will live for at least 100−200 yr, but only 50% have a predicted lifetime above 500 yr. With no temperatures below 6 K, no first core will live longer than 1000 yr. We note that our predicted first core lifetimes are lower limits, as rotational and magnetic support in 3D setups will stretch the lifetimes by reducing the mass accretion rates.
There is a certain number of coloured globules which lie outside of the grey area, at the lowradius (highdensity) end of the dataset. We did not extend our parameter space too far in that direction, as we deemed initial densities above ~10^{16} g cm^{3} to be too high for realistic initial conditions. In hindsight, it appears that such high densities are observed in molecular clouds, but are probably the result of strong compression from interstellar turbulence, shocks, or colliding flows. Exploring the highmass highdensity regime with our simulations may not be the most interesting, or realistic, as turbulence most probably induces low levels of fragmentation of the globule in highmass star formation (Commerçon et al. 2011; Palau et al. 2014), which cannot be captured with a 1D model. However, the very lowmass end of the highdensity region may represent a pathway to forming brown dwarfs. It is always difficult to form brown dwarfs in spherically symmetric simulations, as the globules must have either a very low temperature, or a very high density (or both) to be able to collapse with such a small mass^{5}. Indeed, coloured globules with masses of the order of 0.1 M_{⊙} and radii below 3000 AU are absent from the observations. It is in fact believed that one actually starts with a higher mass globule which is more likely to exceed M_{BE}, and a significant portion of the initial mass is subsequently ejected through a powerful outflow, as suggested by the simulations of Machida et al. (2009). We have thus decided not to extend our parameter space to radii below 3000 AU.
We believe that our simulation set adequately covers the conditions typically found in molecular clouds (and that remain reasonable to simulate with a spherically symmetric scheme), and we have even explored temperatures higher than most of the observed globules in the Aquila cloud. Our higher temperature models will be relevant for e.g. inner Galaxy and outer Galaxy clouds that are more irradiated due to nearby massive stars (Lindberg & Jørgensen 2012). The relations found in this study can be used to make evolutionary predictions on observed prestellar systems, such as the time they will take to form a first core, and the lifetimes of these cores. The typical lifetime of the first core will also give an indication of the probability of observing a first core object in a given field. Our analytical predictions suggest that one has a better chance of observing a first Larson core in lowtemperature systems, where the time span between the first and second core formation is the longest.
6. Conclusions
An extensive set of simulations modelling the gravitational collapse of a dense cloud has been carried out. The calculations follow the initial isothermal contraction of the sphere, the transition into the first adiabatic phase, the dissociation of H_{2} molecules, the second collapse, and finally the early stages of the second adiabatic phase. They can also be used to trace the evolution of individual pieces of fluid, as they fall inwards, travel through accretion shocks, and/or get heated via radiation from the central object. The 143 runs show a wide spread in the thermal evolutionaly tracks (ρ_{c} vs. T_{c}) they follow; this is due to differing initial temperatures, and also to different optical depths which regulate the radiative cooling. The models also revealed that while some dispersion was observed in the density, velocity, or temperature radial profiles outside of the first core accretion shock, the profiles showed very little scatter inside hydrostatic body where they seemed to “forget” their initial conditions.
Considering the wide ranges of masses and densities covered by the parameter space, the simulations produced first and second cores with remarkably similar properties. The average first core mass and radius are ×10^{2}M_{⊙} and AU, while the average second core mass and radius are and ×10^{3} AU. Within less than an order of magnitude, these quantities were basically independent of the initial conditions. In our simulations, the first cores were observed to contract for a brief period before expanding until the accretion shock at their border became radiatively efficient, where the (mostly kinetic) energy of the infalling material impacting onto the core surface could be radiated away.
A curious behaviour was observed when adopting a uniform initial density profile instead of the BE solution; the most unstable runs were found to be devoid of a first Larson core. After a careful study of this phenomenon, we concluded that it was an artefact of the initial conditions where the freefall time is the same in all parts of the cloud. Along with the fact that it seems unlikely for a uniform density body to form in nature, we believe the BE simulations are a better representation of the evolution of collapsing dense clouds.
A comprehensive study of the lifetimes of the first Larson cores (i.e. the time elapsed between the formation of the first and second cores) showed that they strongly depended on the first core mass accretion rate. Two approximate predictive laws linking the initial conditions to the first core lifetimes were also derived. We note that the lifetimes given here are lowerbounds on realistic estimates as magnetic and kinetic support (rotational or turbulent) in 3D situations will slow down the core evolution, reducing the mass accretion rates. The flow of mass will also proceed via the formation of accretion discs, as opposed to covering the entire 4π steradian surface of the core. It is, however, interesting to mention that Tomida et al. (2013) found first core lifetimes to range between 700 and 1000 yr (ϵ = 0.8), while Bate et al. (2014) observed 50−100 yr (ϵ = 0.14) in their simulations, which is not too different from our 1D results, even though both works consider a 3D magnetised cloud with rotation. Magnetohydrodynamic simulations (ideal or with Ohmic dissipation) of rotating clouds probably result in similar first core lifetime because of the efficient angular momentum transport provided by magnetic braking. Without magnetic fields, the lifetime can be considerably longer. Recent nonideal MHD simulations including ambipolar diffusion (Tomida et al. 2015) show longer lifetimes, as the resistive effects suppress angular momentum transport.
The quantity ρ_{ad} has been given the fixed value of 10^{13} g cm^{3}. This actually varies between 10^{13} and 10^{14} g cm^{3} (see Fig. 2a) which will change the vertical scaling somewhat for different runs, but it will not affect the slope.
As Masunaga & Inutsuka (2000) show, it may take some time before the observable properties of the protostellar system, such as the spectral energy distributions, are altered. They predict that this does not happen immediately after the second core is formed, but well into the main accretion phase.
These starless density fluctuations in molecular clouds are often called “dense cores”, which represent an evolutionary stage analogous to the initial conditions in our setup. We shall call them “globules” here, to make a distinction between these prestellar cloud cores and the first and second Larson cores, which are much smaller and denser objects.
This is precisely the origin of the lowmass turnover in the initial mass function (Padoan & Nordlund 2002).
Acknowledgments
N.V. gratefully acknowledges support from the European Commission through the Horizon 2020 Marie SkłodowskaCurie Actions Individual Fellowship 2014 programme (Grant Agreement No. 659706). T.H. is supported by a Sapere Aude Starting Grant from the Danish Council for Independent Research. We thank B. Commerçon and J. Masson for very useful discussions during the writing of this paper. We would also like to express our gratitude to the referee, whose extremely constructive comments were much appreciated and have greatly improved this manuscript. This work made use of the astrophysics HPC facility at the University of Copenhagen, which is supported by a research grant (VKR023406) from VILLUM FONDEN. We finally thank the Service d’Astrophysique, IRFU, CEA Saclay, and the Laboratoire Astrophysique Interactions Multiéchelles, France, for granting us access to the supercomputer IRFUCOAST where some of the calculations were also performed.
References
 André, P., Di Francesco, J., WardThompson, D., et al. 2014, Protostars and Planets VI, 27 [Google Scholar]
 Badnell, N. R., Bautista, M. A., Butler, K., et al. 2005, MNRAS, 360, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R., Tricco, T. S., & Price, D. J. 2014, MNRAS, 437, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Berthon, C., Charrier, P., & Dubroca, B. 2007, J. Sci. Comp., 31, 347 [CrossRef] [Google Scholar]
 Bonnor, W. B. 1956, MNRAS, 116, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Hennebelle, P., & Henning, T. 2011, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Launhardt, R., Dullemond, C., & Henning, T. 2012a, A&A, 545, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Levrier, F., Maury, A. J., Henning, T., & Launhardt, R. 2012b, A&A, 548, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubroca, B., & Feugeas, J.L. 1999, C. R. Acad. Sci. Paris, 329, 915 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Dzyurkevich, N., Commerçon, B., Lesaffre, P., & Semenov, D. 2016, A&A, submitted [arXiv:1605.08032] [Google Scholar]
 Ebert, R. 1955, Z. Astrophys., 37, 217 [Google Scholar]
 Federrath, C. 2015, MNRAS, 450, 4035 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Schrön, M., Banerjee, R., & Klessen, R. 2014, ApJ, 790, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Foster, P. N., & Chevalier, R. A. 1993, ApJ, 416, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Frimann, S., Jørgensen, J., & Haugbølle, T. 2016a, A&A, 587, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frimann, S., Jørgensen, J. K., Padoan, P., & Haugbølle, T. 2016b, A&A, 587, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Grassi, T., Bovino, S., Haugboelle, T., & Schleicher, D. R. G. 2016, MNRAS, 466, 1259 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Charbonnel, C. 2013, EAS Pub. Ser., 62 [Google Scholar]
 Hincelin, U., Commercon, B., Wakelam, V., et al. 2016, ApJ, 822, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Hopkins, P. F., & Lee, H. 2016, MNRAS, 456, 4174 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, C. 1977, ApJ, 218, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J.h., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Könyves, V., André, P., Men’shchikov, A., et al. 2015, A&A, 584, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spectr. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Li, S., Frank, A., & Blackman, E. G. 2014, MNRAS, 444, 2884 [NASA ADS] [CrossRef] [Google Scholar]
 Li, H.B., Yuen, K., Otto, F., et al. 2015, Nature, 520, 518 [NASA ADS] [CrossRef] [Google Scholar]
 Lindberg, J. E., & Jørgensen, J. K. 2012, A&A, 548, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Machida, M. N., Inutsuka, S.i., & Matsumoto, T. 2009, ApJ, 699, L157 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masunaga, H., & Inutsuka, S.i. 1999, ApJ, 510, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Masunaga, H., & Inutsuka, S.i. 2000, ApJ, 531, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Masunaga, H., Miyama, S. M., & Inutsuka, S.i. 1998, ApJ, 495, 346 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, T., Dobashi, K., & Shimoikura, T. 2015, ApJ, 801, 77 [NASA ADS] [CrossRef] [Google Scholar]
 McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Mihalas, B. D. 1984, in Foundations of Radiation Hydrodynamics (New York: Oxford University Press), 557−585 [Google Scholar]
 Nordlund, Å., Haugbølle, T., Küffmeier, M., et al. 2014, IAU Symp., 299, 131 [NASA ADS] [Google Scholar]
 Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Haugbølle, T., & Nordlund, Å. 2014, ApJ, 797, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Palau, A., Estalella, R., Girart, J. M., et al. 2014, ApJ, 785, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Penston, M. V. 1969, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Phillips, A. C. 1999, The Physics of Stars, 2nd edn (WileyVCH) [Google Scholar]
 Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2012, MNRAS, 423, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Seifried, D., SánchezMonge, Á., Walch, S., & Banerjee, R. 2016, MNRAS, 459, 1892 [NASA ADS] [CrossRef] [Google Scholar]
 Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Stahler, S. W., & Palla, F. 2005, The Formation of Stars (WileyVCH) [Google Scholar]
 Stahler, S. W., Palla, F., & Salpeter, E. E. 1986, ApJ, 302, 590 [NASA ADS] [CrossRef] [Google Scholar]
 Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., & Masson, J. 2012, A&A, 543, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaytet, N., Chabrier, G., Audit, E., et al. 2013a, A&A, 557, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaytet, N., González, M., Audit, E., & Chabrier, G. 2013b, J. Quant. Spectr. Rad. Transf., 125, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Visser, R., Bergin, E., & Jørgensen, J. 2015, A&A, 577, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wakelam, V., Loison, J.C., Herbst, E., et al. 2015, ApJS, 217, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Resolution study
Fig. A.1 Resolution study for run 69. Panels a) and c) show the temperature and entropy, respectively, of the cell at the centre of the system as a function of density using 256, 512, 1024, 2048, 4096, and 8192 cells (see legend). In panel a), the thin solid lines are for simulations with the standard timestep calculation based on relative variations of primary variables from one time step to the next, while the thick dashed lines are for the reruns when the values of dt are enforced according to the central gas density (see text). Panel b) shows the central temperature as a function of time. Panels d) and f) show profiles of the gas density and velocity, respectively, as a function of radius, at the end of the simulations. Panel e) shows the integration timesteps in the different simulations as a function of central density. The thick grey dashed line represents the forced values of dt we have chosen for our temporal resolution reruns (see text). 

Open with DEXTER 
Appendix A.1: Spatial resolution
A resolution study was carried out to ensure the results presented in this work are numerically converged. Indeed, the evolution of the first and second Larson cores are very sensitive to numerical resolution, as enough cells are needed to adequately sample not only the Jeans length but also the accretion shocks at the core borders. A typical run was chosen (run 69) and resimulated at different resolutions. The same scheme to increase grid cell size with increasing radius was used in all runs, but to allow for a fair comparison, the value of α was changed for each resolution. Indeed, for 4096 cells, the size ratio between the innermost and outermost cells is 1.0008^{4095} ≃ 26.4. If we used the same α for 256 cells, this ratio would stand at 1.0008^{255} ≃ 1.23. By considering only two cells in the grid, the first cell has size Δr and the second has size (1 + α_{0})Δr. If one wishes to double the resolution, i.e. use 4 cells for the same mesh size and a constant size ration (1 + α_{1}), each coarse cell is split into two small cells. Let δ_{i} = 1 + α_{i}, where i = 0,1. This gives us the relations
(A.1)where dr is the size of the first cell in the finer resolution mesh. This yields , or . For our resolution study, we adopt the final relation (A.2)where N_{z} is the number of grid cells. It is now easy to see that the size ration between the innermost and outermost cells (A.3)is almost independent of N_{z}, as long as the number of cells is large enough.
Fig. A.2 Evolution sequence for two simulations with different resolutions. The first and third columns are for 256 grid cells, while 4096 cells have been used in the second and fourth columns. Panels a) to n) show radial profiles of various quantities, while panels o) and p) show a small portion of the thermal evolutions from Fig. A.1a. The different colours represent a succession of epochs to illusrate the evolution of the system around the “knee” which appears in lowresolution runs. The numbers in panels o) and p) show the order in time of the colours. In panels g) and h), dashed lines represent negative values of the radiative flux (i.e. a flux pointing towards the centre of the system). 

Open with DEXTER 
The thermal evolutions at the centre of the system (temperature as a function of density) are shown in Fig. A.1a (thin solid lines). While the first core is accreting, for central densities (in g cm^{3}) 10^{10} ≲ ρ_{c} ≲ 10^{8}, the low resolution runs display a strong deviation from adiabaticity (we will refer to this as a “knee” below). During the second collapse (10^{8} ≲ ρ_{c} ≲ 10^{4}), all runs follow very similar adiabats, but the spread is again large once the second core has formed. This causes a significant spread in second core entropy in panel c, and it is thus of prime importance that the simulations have the necessary resolution. The results appear to converge when 2048 cells and above are used.
Because the simulations are stopped as soon as a “hydrostatic bounce” occurs inside the second core (as soon as the density at the centre of the system decreases from one timestep to the next), the resolution also affects the final central densities and temperatures in the different runs. Table A.1 lists the final density, temperature and entropy at the centre of the system, for the different resolutions. It also lists in the last column the first core lifetime. This now suggests that a minimum resolution of 4096 is probably needed. This is confirmed by the temporal evolution of the central temperature in Fig. A.1b. All simulations show a first hydrostatic bounce when the first core is formed (T_{c} ~ 200−300 K), but it occurs at different times. Simulations with lower resolutions then display several more bounces before proceeding to the second collapse, while the most resolved cases reach it in a much more “direct” manner. They also get there much faster, as illustrated by the wide spread in first core lifetimes τ_{fc} (i.e. the time elapsed between the formation of the first and second cores) in the final column of Table A.1. Finally, in panels d and f, density and velocity profiles, respectively, show how well the accretion shocks at the first and second core borders are resolved. Not only is it important to have an accurate description of the evolution of the gas at the centre of the system, it is also necessary to retain a good enough resolution at the first and second core borders where the inflowing gas passes through sharp accretion shocks where density and temperature rise very abruptly.
The origin of the “knee” seen in Fig. A.1a can be understood from looking at a succession of profiles, in two runs with different resolutions. This is shown in Fig. A.2 for 256 and 4096 cells, where different snapshots in time are represented by different colours. Since runs with different resolutions collapse and evolve at slightly different rates (see Fig. A.1b for example), it does not make sense to compare profiles at exactly the same times, central temperatures or central densities. Instead, we have tried to select snapshots where the systems are approximately in the same states. The first snapshot (green) is at a time just when the first core accretion shock has formed. The second (pink) is when the gas velocity just behind the shock becomes positive and the core starts to expand. The third (blue) was taken when the core has settled and the gas velocities are once again negative everywhere in the grid. The final two snapshots (orange then black) show the subsequent evolution of the system, at times just before and just after the “knee” in the lowresolution run. We can see in panels a, b, e, f, o, p that at the time of first core formation, the highresolution run reaches higher densities and temperatures in the centre. Nevertheless, the different velocity profiles are very similar (panels c and d). We note that shortly after the formation of the first core, both high and lowresolution simulations show a small bump in gas temperature just behind the accretion shock, around 2−4 AU (best visible in the pink curve). The optical depth in both cases is well above unity (panels m and n) and one can easily approximate the radiative flux to always be pointing in the direction of a temperature gradient^{6}. This can indeed be observed in panels g and h, where the radiative flux becomes negative between the summit of the temperature bump and the temperature minimum just to the left of the bump (r~ 2 AU)^{7}. What happens next is interesting. Because the resolution is lower in the simulation, the optical depth per cell is higher than in the run by about an order of magnitude (panels k and l). The heating from the core gets trapped and the temperature in the bump rises. Contrastingly, in the highresolution run, the heat is mostly radiated ahead of the shock; we can indeed see in panel f that the bump has hardly changed temperature between the pink and blue curves, but a sizeable radiative precursor has developed ahead of the shock (the temperature has risen by about a factor 5 between 3 and 10 AU). Because the bump has not risen in temperature in run , the temperature in the centre of the system is still higher than the peak of the bump, and the radiative flux only shows a small sign inversion between 1.5−2.5 AU. In the calculation, the bump temperature exceeds the temperature at the centre of the core, and the radiative flux is now negative throughout almost the entire core (i.e. to the left of the accretion shock; dashed lines in panel g). From here on, the flux continuously travels inwards, and begins to heat the core, from outside towards the centre. Indeed, in the orange curve of panel e, we can see that the bump widens mostly in the inner direction. This heating Marshak wave travels towards smaller radii and eventually hits the centre of the system where the temperature rises abruptly. This causes the “knee” in panel o. The last snapshot (black) shows profiles just at the end of the inward heating phase, right before the temperature profile inside the first core becomes flat again. In the case of 4096 grid points, the radiative flux is never fully directed towards the centre, and the temperature at the centre rises smoothly as time progresses. This shows that it is essential to resolve the characteristic lengths of all the radiative processes in our system. Typical radiation schemes usually require the optical depth per cell to remain below unity, but the asymptotic preserving scheme of Berthon et al. (2007) used in our code allows us to go to values ten times higher. The above analysis however shows that the solver cannot go much higher than κ_{R}ρdr ~ 10.
We have shown that the results converge for 4096 cells and above, and all 143 runs in the present work were thus run using 4096 cells. We are aware that this may change for a different set of initial conditions. We have carried out our resolutions study on a couple of other runs and found that the results were robust for N_{z} = 4096. Even though we have not done this for all the 143 runs, we are confident that the conclusions reached in this work are not affected by resolution effects.
Density ρ_{c}, temperature T_{c}, and entropy S_{c} at the centre of the system, at the end of the simulations, as a function of resolution.
Appendix A.2: Temporal resolution
In our timeimplicit scheme, we limit the integration timesteps by constraining the variations from one time step to the next of each primary variable to a maximum of 10% (see Sect. 2). Using a finer mesh also implies that smaller timesteps are used, and we may be seeing a double effect in the spatial resolution study above, as the lower spatial resolution runs may in fact also suffer from a lower temporal resolution. To verify this, we have rerun the simulations by forcing the timesteps to values similar to the highest resolution run (8192 cells). The timesteps of the above runs as a function of central density are plotted in Fig. A.1e. The timesteps in all simulations follow a similar trend. They start at a very low value, set by the radiation CFL which defines the initial timestep to avoid any convergence issues at the beginning of the runs. They rise sharply until they reach about 1000 yr, where they become limited by strong variations of primary variables. The timesteps gradually fall as the system collapses. We chose a piecewise law for the value of dt as a function of central density which is represented by the thick grey dashed line in Fig. A.1e, which remains below all the curves at all times. The results of the reruns with forced timesteps are shown in Fig. A.1a (thick dashed lines). We see that the thermal evolutions are virtually identical to the original runs, and a higher temporal resolution does not alleviate the departures from adiabaticity in low spatial resolution runs.
Appendix B: Correlations weighting scheme
When computing simple linear correlations between two quantities x and y where all the points were given the same weight of 1, we found that several correlation coefficients r_{xy} were poorly reflecting (x,y) relations which appeared quite strong, apart from a few outliers. Others seemed too high for variables which did not seem connected in any way, and outliers were setting the value of r_{xy}. This is illustrated in Fig. B.1, left column, where three relations with identical unweighted correlations are displayed. It is evident from the figure that the relation in the top panel shows in fact quite a tight correlation between the outer luminosity L_{out} and the first core radius R_{1}, with a few outliers, while the middle and bottom correlations, are relatively weak. The outliers have a strong influence on the final result, and it is well known that visual inspection is usually a mandatory check when working with linear correlations. It is never trivial to robustly identify outliers from a mathematical point of view, while they are in practise quite easy to spot by eye. To avoid having to inspect the 462 correlation plots, we devised a scheme to reduce the influence of the outliers on r_{xy}.
Fig. B.1 Illustration of the weighting scheme for computing correlations between two variables x and y. The top row shows the relation between L_{out} and R_{1}, the middle row is for M_{0} as a function of run number, while the bottom row displays versus T_{0}. In the left column, the light blue dots represent the different simulations (the run numbers are indicated inside each symbol), while the coloured contours show the normalised weights given to the data points. In the right column, the blue solid line represents the cumulative distances for each run, sorted in ascending order. The red dashed line is the derivative of this distance distribution. The thick vertical black solid line marks the position of the first derivative spike above 0.075 (grey dotted line), while the black arrow reads the value of the distance d_{0} on the blue curve at the location of the derivative spike. The insets show the weighting functions w(d,d_{0}). 

Open with DEXTER 
Outliers can be defined as isolated points that are distant from the majority of the other data points, a property that we shall take advantage of in our weighting scheme. To determine if a point is distant from the rest of the group, we computed the cumulative distance of a given point i with all the other points j. We have (B.1)and (B.2)where the min and max values are used to normalise the distances. If a point lies in the middle of the cluster, many of its distances to other points will be rather small, and only a few large distances towards the outliers will be recorded. In the case of an outlier, almost all of the distances computed will be large, amounting to a very high cumulative distance. We now need to choose a weighting function which gives a large weight to points with small distances and viceversa. Many options exist, and common choices often involve 1 /r or 1 /r^{2} laws. We have found these to decrease too fast, giving too much weight to the points lying in the middle of a cluster, yielding too high correlation values computed from only a few data points from the sample. An optimal function should have an almost constant weight of 1 for most of the cluster, and a relatively sharp drop when distances become too large. The function (B.3)has the required properties. The value of a determines the steepness of the decrease, while the value of d_{0} dictates the value of d for which w starts to drop towards zero. The additional constant 0 <b< 1 makes sure that no weight ever goes to zero to avoid discarding any data point completely. The arbitrary choices of a = 10 and b = 0.05 were made after some trial and error on several correlation relations. The value of d_{0} however needs to be different for each correlation plot, as choosing a canonical value of 0.5 may cause many points to be omitted in both a correlated and noncorrelated samples. For instance, in a perfectly correlated dataset, where all the points lie along the same line, the points in the middle of the line will have maximum distances of half the line length, while points at the extremities will have maximum distances of the entire line length, and consequently a larger cumulated distance. These points will be given a very low weight even though they are perfectly aligned with the others. This may not change the overall value of r_{xy} in this idealised case, but it does limit the extent of the region of the parameter space over which r_{xy} is determined. In a slightly less idealised case, where some scatter is present in the data, points far away from the mean of the sample but still lying along the same line may in fact help strengthen the correlation, or at least increase the accuracy of r_{xy}. On the other hand, in the case of a fully uncorrelated sample such as a random circular cloud of points, points towards the edges of the circle will have a much lower weight than the ones close to the centre, for no particular reason. By choosing an arbitrary value for d_{0}, one can strongly influence the statistical results in a biased way.
The method we have devised to objectively determine the value of d_{0} for each plot is illustrated in the right column of Fig. B.1. We first computed the cumulative distances for all the points in the sample, and sorted them in ascending order; this is represented by the blue solid line in panels b, d, and f. If one or more outliers are present in the data, there will be a sharp rise in the blue curve, as d jumps from being in the cluster to lying at a large distance from the majority of the other points. The derivative of the sorted distance array is depicted by the red dashed line. We select the first derivative spike above 0.075 as the limit beyond which we are dealing with outliers (thick vertical black line). By then reading off the cumulative distance of the run responsible for the spike, we finally obtain our value for d_{0}. The small insets in panels b, d, and f show the weighting function for three different values of d_{0}.
We plot the 2D weighting function in the left panels using contours. In the top row, the weight drops moderately as one moves away from the cluster (d_{0} = 0.69), reducing the influence of the outliers in the bottom left corner of the plot. The value of  r_{xy}  has increased from 0.60 to 0.79, now revealing the apparent correlation in the data which was not obvious in the nonweighted (“old”) scheme. In the middle row, the weight falls much more rapidly with increasing d (d_{0} = 0.24) and the derivative quickly spikes above the threshold value of 0.075. This limits the effects of the outliers in the bottom left corner, lowering the value of  r_{xy}  from 0.60 to 0.38. In the bottom row, as the points are pretty much evenly spread out over the viewport, all the cumulative distances are similar, and the derivatives only exceed the threshold for the last few data points. Consequently, the value of d_{0} = 0.92 is large, most of the data points have a weight close to 1, and the value of r_{xy} is virtually unchanged. We have inspected the vast majority of our final correlation plots by eye to make sure the weighting scheme was robust and produced the desired results.
Appendix C: List of run parameters
Table C.1 lists the initial parameters and the properties of the first and second Larson cores for all 143 runs. The columns are (from left to right):
3.: T_{0} Initial cloud temperature
4.: R_{0} Initial cloud radius
5.: Initial cloud central gas density
6.: ϵ = M_{BE}/M_{0} BonnorEbert mass over initial mass
11.: t_{sound} = R_{0}/c_{0} Sound crossing time
14.: First core mass accretion rate
15.: First core accretion luminosity
16.: First core radiative luminosity
17.: τ _{fc} First core lifetime
10.: Second core mass accretion rate
11.: Second core accretion luminosity
In the above, the quantities with a star subscript (⋆) are evaluated at the relevant core border (first or second core), and F represents the radiative flux. The simulation end time for each run is equivalent to t_{ff3} (see Sect. 3.5). (Interactive image with hyperlinks available in the PDF version.)
Simulation initial parameters and core properties for all 143 runs.
All Tables
Density ρ_{c}, temperature T_{c}, and entropy S_{c} at the centre of the system, at the end of the simulations, as a function of resolution.
All Figures
Fig. 1 Model parameter space for 6 different radii: R_{0} = 3000,5000,7500,12 000,20 000,30 000 AU. The hatched area represents clouds that are stable against gravitational collapse, while the colours code for different values of ϵ for clouds in the unstable zone. The ϵ = 1 separating line between hatched and coloured areas is defined in Eq. (5). Each run number links to the corresponding publicly available dataset. (Interactive image with hyperlinks available in the PDF version.) 

Open with DEXTER  
In the text 
Fig. 2 Top row: evolution of the central temperature (a) and entropy (b) as a function of density for all 143 runs. The lines have been colourcoded according to the cloud’s initial temperature T_{0} (see legend in panel b)). Bottom row: c) evolution of the central temperature as a function of density for runs 30, 3, and 140 (see legend). d) Optical depth of the systems as a function of radius for runs 30, 3, and 140 (same colours as in panel c)), when the central density has just reached log (ρ_{c}) = −12.75 g cm^{3} (solid), and when the central temperature has reached log (T_{c}) = 1.2 K (dashed). 

Open with DEXTER  
In the text 
Fig. 3 Evolution track in the densitytemperature plane for 10 different pieces of fluid in runs a) 69, b) 86, c) 71, and d) 115. In all panels, 10 cells were selected and followed throughout the simulation. The number of the cell in the mesh is given for each track in the bottomright inset, along with the cell’s radius at the start of the simulation (“s”) and at the end (“e”). While the cell numbers are the same in all panels, the cell radii change. The colour codes for the time in kyr; the colourbar is different in all panels, but they all display the final 700 yr of evolution. The pink dots show the position of each piece of fluid at the end of the simulation. The dashed grey line traces the final (ρ,T) profile in the entire grid, which connects the pink dots. The inset in the top left corners gives a reminder of the parameters of each run. 

Open with DEXTER  
In the text 
Fig. 4 Radial profiles of the gas density (a), velocity (b), temperature (c), and luminosity (d) for all 143 runs. The different colours represent different initial temperatures T_{0} (see legend). The vertical dotted lines mark the positions of the average first and second core radii. 

Open with DEXTER  
In the text 
Fig. 5 First and second core mass (a) and radius (b) for all 143 runs. Each mass and radius has a starting and an end point. The starting point represents the quantity when the core has just formed, while the end point is at the end of the simulation. A small arrow in the middle of each segment points in the direction from start to end. The vertical alternating grey and white bands indicate changes in initial radius (see panel b) for values of R_{0}). 

Open with DEXTER  
In the text 
Fig. 6 a) First core radius as a function of time, where t = 0 corresponds to the time when the first core accretion shock is formed. The different colours represent different initial cloud temperatures T_{0} (see legend). The thick black line is for run 69. b) First core radius as a function of central density (solid lines), alongside the (rescaled) analytical predictions of Masunaga et al. (1998, dashed lines). The colour coding is the same as in panel a). c) Difference between the downstream gas pressure and the upstream ram pressure. d) Downstream gas pressure (blue), upstream ram pressure (red) and radiation pressures (greens) just upstream (subscript 1) and downstream (subscript 2) of the first core accretion shock, in run 69. The upstream radiative fluxes, scaled with the local gas velocity for dimensional coherence, are also plotted for relative comparison in the total energy budget (black). e) First core accretion shock compression ratio for run 69 measured in the simulation (black) and predicted for a nonradiative shock (red). 

Open with DEXTER  
In the text 
Fig. 7 Ratio of the time it takes in each run until the first core is formed (t_{core1}) to various collapse time estimates, as a function of ϵ. The blue squares represent the freefall time computed using the cloud’s inital central density ρ_{c} (t_{ff1}), while the green triangles show the results from an alternate definition of the freefall time based on the cloud’s average density (t_{ff2}). The black stars represent the cloud’s initial sound crossing time (t_{sound}). Finally, the best estimate is obtained with a modified version of t_{ff1} which takes into account how unstable the cloud initially is to gravitational collapse by introducing the ϵ parameter (t_{ff3}). 

Open with DEXTER  
In the text 
Fig. 8 a) First core lifetime τ_{fc} as a function of ϵ measured in the simulations (circles). The dashed and thin solid lines represent our analytical predictions (see text). b) First core lifetime as a function of mass accretion rate for all runs. The black dashed line portrays a simple estimate from dimensional arguments (see text). c)–i) Ratio of the mass accretion rate measured in the simulation to various predictions Ṁ_{pred} from numerical and similarity solutions, as a function of ϵ. c) The hearts represent the similarity solution of Larson (1969) and Penston (1969), the stars are the solution of Shu (1977), while the triangles are from the results of Hunter (1977). In panels d) to f), the Ṁ_{pred} estimate is obtained from dividing the initial cloud mass M_{0} by the freefall times t_{ff1}, t_{ff2}, and t_{ff3}. g) The hearts represent a modified version of Penston’s formula, which is equivalent, barring a scaling factor, to the prediction using the cloud sound crossing time (diamonds). h) The circles represent the first best fitting formula Ṁ_{A} we have used to predict the lifetimes. i) The squares show the second mass accretion rate estimate Ṁ_{B} we have selected for our analysis. The colours in all panels code for different initial temperatures (see legend in panel a)). 

Open with DEXTER  
In the text 
Fig. 9 Correlation matrix for the parameters and core properties in Table C.1. The correlations of the type r_{xx} in the white diagonal have been ignored. For r_{xy}< 0.5, the values are given in the cells for indicative purposes, and have intentionally been made less visible to avoid cluttering the figure with too many insignificant values. The interested reader should be able to see the numbers by zooming in enough on the figure. Each square links to a correlation plot stored in the online database. (Interactive image with hyperlinks available in the PDF version.) 

Open with DEXTER  
In the text 
Fig. 10 a) Final density profiles for BonnorEbert (dashed) and uniform density (solid) setups with initial temperatures of 5 K (blue) and 30 K (red). b) Same as in a) but for the gas velocity. c) Ratio of the gas thermal pressure to ram pressure as a function of radius. The arrows point to the absence of a first Larson core in the most unstable uniform density run. 

Open with DEXTER  
In the text 
Fig. 11 Comparison between the radii, temperatures and masses of the prestellar cores in the Aquila molecular cloud and the parameter space covered by the initial conditions of our set of numerical simulations. In each panel, the empty circles represent observed cores with ϵ > 1, while the coloured squares are for cores with ϵ < 1. The colour of the squares code for their predicted first core lifetimes. The grey area deliniates the simulation parameter space: every corner or intersection of a dotted line with either another dotted line or the thick grey boundary represents one simulation. 

Open with DEXTER  
In the text 
Fig. A.1 Resolution study for run 69. Panels a) and c) show the temperature and entropy, respectively, of the cell at the centre of the system as a function of density using 256, 512, 1024, 2048, 4096, and 8192 cells (see legend). In panel a), the thin solid lines are for simulations with the standard timestep calculation based on relative variations of primary variables from one time step to the next, while the thick dashed lines are for the reruns when the values of dt are enforced according to the central gas density (see text). Panel b) shows the central temperature as a function of time. Panels d) and f) show profiles of the gas density and velocity, respectively, as a function of radius, at the end of the simulations. Panel e) shows the integration timesteps in the different simulations as a function of central density. The thick grey dashed line represents the forced values of dt we have chosen for our temporal resolution reruns (see text). 

Open with DEXTER  
In the text 
Fig. A.2 Evolution sequence for two simulations with different resolutions. The first and third columns are for 256 grid cells, while 4096 cells have been used in the second and fourth columns. Panels a) to n) show radial profiles of various quantities, while panels o) and p) show a small portion of the thermal evolutions from Fig. A.1a. The different colours represent a succession of epochs to illusrate the evolution of the system around the “knee” which appears in lowresolution runs. The numbers in panels o) and p) show the order in time of the colours. In panels g) and h), dashed lines represent negative values of the radiative flux (i.e. a flux pointing towards the centre of the system). 

Open with DEXTER  
In the text 
Fig. B.1 Illustration of the weighting scheme for computing correlations between two variables x and y. The top row shows the relation between L_{out} and R_{1}, the middle row is for M_{0} as a function of run number, while the bottom row displays versus T_{0}. In the left column, the light blue dots represent the different simulations (the run numbers are indicated inside each symbol), while the coloured contours show the normalised weights given to the data points. In the right column, the blue solid line represents the cumulative distances for each run, sorted in ascending order. The red dashed line is the derivative of this distance distribution. The thick vertical black solid line marks the position of the first derivative spike above 0.075 (grey dotted line), while the black arrow reads the value of the distance d_{0} on the blue curve at the location of the derivative spike. The insets show the weighting functions w(d,d_{0}). 

Open with DEXTER  
In the text 