Why Do Some Cores Remain Starless?

Abstract

Prestellar cores, by definition, are gravitationally bound but starless pockets of dense gas. Physical conditions that could render a core starless (in the local Universe) is the subject of investigation in this work. To this end, we studied the evolution of four starless cores, B68, L694-2, L1517B, L1689, and L1521F, a VeLLO. We demonstrate: (i) cores contracted in quasistatic manner over a timescale on the order of ~ 10⁵ yr. Those that remained starless briefly acquired a centrally concentrated density configuration that mimicked the profile of a unstable BonnorEbert sphere before rebounding, (ii) three cores viz. L694-2, L1689-SMM16, and L1521F remained starless despite becoming thermally super-critical. By contrast, B68 and L1517B remained sub-critical; L1521F collapsed to become a VeLLO only when gas-cooling was enhanced by increasing the size of dust-grains. This result is robust, for other starless cores viz. B68, L694-2, L1517B, and L1689 could also be similarly induced to collapse. The temperature-profile of starless cores and those that collapsed was found to be radically different. While in the former type, only very close to the centre of a core was there any evidence of decline in gas temperature, by contrast, a core of the latter type developed a more uniformly cold interior. Our principle conclusions are: (a) thermal super-criticality of a core is insufficient to ensure it will become protostellar, (b) potential star-forming cores (the VeLLO L1521F here), could be experiencing dust-coagulation that must enhance gasdust coupling and in turn lower gas temperature, thereby assisting collapse. This also suggests, mere gravitational/virial boundedness of a core is insufficient to ensure it will form stars.

1 INTRODUCTION

It is well-known that a number of cores do in fact remain starless. Key physical characteristics about these cores to have emerged from various detailed observational surveys are—(a) their longevity or in other words, the fact that some cores remain starless for a period significantly longer than others do (e.g. Ward-Thompson et al. 2007), and (b) that cores typically are not pockets of isothermal molecular gas, but are colder close to their respective centres. Also, gas within cores is usually subsonic or transonic at best (e.g. Evans et al. 2001; Crapsi et al. 2004, 2007, di Francesco et al. 2007; Schnee et al. 2013). Observational estimates of the typical age of cores derived from detailed studies of their chemical composition and the strength of the magnetic field threading them are on the order of ~ 10⁶ yr, which is on the order of a free-fall time for a typical prestellar core (e.g. Beichman et al. 1986; Kirk, Ward-Thompson, & André 2005; Ward-Thompson et al. 2007; see e.g. Maret, Bergin, & Tafalla 2013 for an estimate of the chemical age). Broadly speaking, models attempting to explain the evolution of prestellar cores may be classified into two paradigms: (1) those that propose quasi-static evolution of magnetised cores, and (2) those in which cores form and evolve on a relatively short timescale in turbulent clouds, the so-called paradigm of rapid star formation.

Depending on the strength of magnetic field threading a core, it may be further classified as magnetically sub-critical(mass to flux ratio less than unity), or super-critical(mass to flux ratio greater than unity). Sub-critical cores collapse on the ambipolar diffusion timescale (e.g. Shu et al. 1987; Mouschovias 1991), whilst those that are super-critical, in the absence of turbulence, collapse essentially on the free-fall timescale (e.g. Nakano 1998). This model, however, encounters a number of difficulties: (i) the sub-critical nature of cores is inconsistent with the findings of a detailed study of a handful of isolated cores that show these cores are threaded by only a relatively weak magnetic field (e.g. Kirk, Ward-Thompson, & Crutcher 2006), (ii) the ambipolar diffusion timescale, typically on the order of ~ 10⁷ yr—the lifetime of a molecular cloud (MC) (e.g. Nakano 1998; Ciolek & Basu 2001), is much higher than the estimates of typical core lifetime, and (iii) it is also inconsistent with the paradigm of rapid star formation, reinforced by the observations of several nearby star-forming regions where stellar populations are found to have ages significantly smaller than the crossing time for their respective host cloud (e.g. Hartmann, Ballesteros-Paredes, & Bergin 2001, Elmegreen 2000).

In the dynamic picture of star formation, on the other hand, supersonic turbulence, perhaps moderated by a weak magnetic field, largely regulates formation of prestellar cores as they condense out via amplification of local perturbations in the density field of MCs. Also, since turbulence decays rapidly, over less than a free-fall time, cores supported by turbulence could possibly explain various features of starless cores (e.g. Mac Low et al. 1998; Ballesteros-Paredes et al. 2003, 2006). Cores, according to this model, are in approximate (thermal) pressure-equilibrium with the ambient medium and have a density distribution that reflects the profile of a pressure-confined Bonnor–Ebert sphere (BES). One of the difficulties with this model of core formation is the transient nature of resulting cores on account of being merely in hydrostatic equilibrium (e.g. Mac Low & Klessen 2004), which is inconsistent with typical properties of starless cores. Furthermore, simulations such as those by Walch et al. (2010), for instance, demonstrate that turbulence can only delay collapse in gravitationally super-critical cores, but cannot arrest it. Evidently, a solution to the problem of starless nature of cores must be sought beyond these tested ideas such that their physical properties noted above can be reconciled.

Evolution of cores has been studied analytically by a number of authors over the last four decades. Some of the earliest models suggested by e.g. Bodenheimer & Sweigert (1968), Larson (1969), and Penston (1969) were criticised later by Shu (1977) on grounds of idealised(artificial) initial conditions. The Shu model that envisaged a test core modelled as a unstable BES which would then collapse in the inside-out manner was itself later deemed physically unrealisable, for it is difficult to reconcile the assembly of a core that is critically stable at best, and unstable at worst (Whitworth & Summers 1985; hereafter WS). WS demonstrated analytically that a pressure-confined BES could indeed be driven to collapse in outside-in fashion by a compressional wave triggered by raising sufficiently the magnitude of the externally confining pressure, $P_{\text{ext}}$ . This was later demonstrated numerically by Hennebelle et al. (2003). Physically, this would be the case in regions of active star formation where energy and momentum is dumped into the ambient medium by ongoing episodes of star formation. Similar arguments were presented by Gómez et al. (2007) who demonstrated that a sufficiently strong inwardly propagating compressional wave could induce a BES to collapse. More recently, Anathpindika and Di Francesco (2013) showed that a core modelled as a non-isothermal BES and with sufficiently cold interiors can also collapse in this manner to become protostellar. Other cores in their work where gas close to the core-centre was not sufficiently cold remained starless and in fact, oscillated radially. Similarly, Kaminski et al. (2014) demonstrated that BESs could be similarly induced to collapse by increasing the magnitude of P _ext, effected by raising the density of the confining medium. They also showed that BESs enveloped by a hot tenuous medium would likely evolve in quasi-static manner before eventually collapsing.

Models invoking detailed radiative transfer modelling of starless cores include those suggested by Keto & Caselli (2008, 2010). These authors classified starless cores into two classes viz., thermally sub-critical(central density, n_c ≲ 10⁵ cm⁻³) and super-critical (n_c ≳ 10⁵ cm⁻³) cores. They argue, gas temperature in centrally dense cores is lowered more efficiently via collisional coupling between gas and dust which renders such cores thermally super-critical and therefore, more susceptible to collapse which is likely to proceed in the outside-in manner (e.g. Williams et al. 1999; Caselli et al. 2002; Schnee et al. 2007). On the other hand, thermally sub-critical cores are likely to be relatively stable against self-gravity and may exhibit radial oscillations (e.g. Lada et al. 2003; Aguti et al. 2007). These models usually assume a density distribution to model the temperature profile of a typical starless core. Other relatively simple analytic models such as the Plummer-type suggested by Whitworth & Ward-Thompson (2001) produces rapid density enhancement on a timescale much shorter than the observational estimates of contraction timescale.

In the present work, we therefore propose to numerically investigate the evolution of a set of cores whose physical characteristics are well-known. We are particularly interested in studying the conditions under which a core is likely to be rendered unstable against self-gravity. To this end, we will self-consistently deduce the temperature profile for a core that remains starless and one that collapses. Also, we will study the temporal excursion of the virial components of our test cores to examine if gravitational and/or virial boundedness is sufficient to ensure that a prestellar core becomes protostellar. This article is organised as follows: In Section 2, we will discuss the choice of our initial conditions and the numerical scheme employed to investigate the problem. We will present our results and discuss them in respectively Sections 3 and 4 before concluding in Section 5.

2 PHYSICAL DETAILS

2.1. The model core : A brief review of the formation of cores and their density profile

Test cores in the present work were modelled by examining the density profile of a typical core that forms in a MC. We therefore refer to our recent work (Anathpindika 2015), where the formation of cores had been demonstrated via an interplay between the thermal instability (TI) and the gravitational instability (GI), in a purely hydrodynamic realisation of a self-gravitating MC. For a comprehensive discussion of the subject, the interested reader is referred to a recent review by Offner et al. (2014). Lack of adequate computational resources, however, prevented those calculations from continuing further to track the evolution of individual cores. We therefore extracted cores emblematic of those found in typical star-forming regions from an earlier simulation reported in Anathpindika (2015). Cores could be isolated (e.g. B68), or be part of an elongated filament-like natal cloud such as the countless number of them detected in the Ophiuchus MC or the Taurus MC amongst many others.

The images on panels of Figure 1 show the two representative cores selected from a simulation reported in Anathpindika (2015), of which, that on the left-hand panel is an isolated core whilst the one on the right-hand panel belongs to a larger natal filament-like clump. The density profile of each of the two has been shown in Figure 2. Interestingly, and equally significantly the gas distribution in either core, as is visible from the respective panels of Figure 2, is very well approximated by the density profile of a stable BES of radius, ξ _b = 3. This is also consistent with the density profiles reported by Teixeira, Lada, & Alves (2005) for their sample of cores observed in the Lupus 3 MC. We use this finding to model our test cores in this work.

Figure 1. Shown here are rendered density plots from an earlier work (Anathpindika 2015), of a core that formed in a small section of the fragmented cloud. Picture on the left-hand panel is an example of a core in relative isolation, dissociated from any large clump, whilst that on the right-hand panel shows a core that belongs to a contiguous, elongated clump in a different region of the same fragmented cloud; time in units of Myr has been marked on the top left-hand corner of each panel.

Figure 2. Distribution of gas density in the respective cores shown in Figure 1 above. Dots here represent density of individual gas particles that assembled these cores and interestingly, the density profile of a stable Bonnor– Ebert sphere having radius, $\xi _{\text{b}}$ = 3 (see text below), appears to fit the density distribution of either cores very well.

The Bonnor–Ebert sphere (BES)

The maximum stable mass of the isothermal (or the Bonnor-Ebert (BE)) sphere of a given radius, ξ _b , is

(1) $$\begin{equation} M_{\text{BE}}\equiv M(\xi =\xi _{\text{b}}) = \frac{a_{0}^{3}}{(4\pi G^{3})^{1/2}}\Big (\xi ^{2}\frac{d\psi }{d\xi }\Big )^{1/2}_{\xi _{\text{b}}} \end{equation}$$

(Bonnor 1956; Ebert 1955). Here, G is the gravitational constant, $\xi =\frac{r}{R_{0}}$ is the dimensionless radius of the BE sphere, where $R_{0}=(\frac{a_{0}^{2}}{4\pi G\rho _{c}})^{1/2}$ , is the physical scale-factor for length and, a ₀, the isothermal sound-speed. The central density, ρ _c , for the BE sphere is defined as

(2) $$\begin{equation} \rho _{\text{c}} \equiv \frac{a_{0}^{6}}{4\pi G^{3}M^{2}}\Big (\xi ^{2}\frac{d\psi }{d\xi }\Big )^{2}_{\xi _{\text{b}}}, \end{equation}$$

so that the radial distribution of density within the BES is

(3) $$\begin{equation} \rho (r) = \rho _{c}\mathrm{e}^{-\psi (\xi )}, \end{equation}$$

M ≡ M(ξ), and ψ is the well-known Lane–Emden function (Chandrasekhar 1939).

Now, Equation (1) can be used to calculate the magnitude of gas temperature required to support a core modelled as a BES of radius, ξ _b , and having mass, $M_{\text{core}}$ , against self-gravity as

(4) $$\begin{equation} T^{\text{exp}}_{\text{gas}}(\xi _{\text{b}}) = \frac{\bar{m}}{k_{\text{B}}}\Big [M_{\text{core}}(4\pi \rho _{\text{c}}G^{3})^{1/2}\Big (\xi \frac{d\psi }{d\xi }\Big )^{-1}_{\xi _{\text{b}}}\Big ]^{2/3}, \end{equation}$$

$\bar{m}$ , being the mean molecular mass of gas. Finally, the magnitude of the pressure, P _ext, confining the BES is simply

(5) $$\begin{equation} P_{\text{ext}}\ \equiv \ P(\xi =\xi _{\text{b}}). \end{equation}$$

2.2. Numerical algorithm

Simulations discussed in this work were developed using the Smoothed particle hydrodynamics (SPH) code, SEREN (Hubber et al. 2011). SEREN is a well-tested, state-of-the-art SPH code that incorporates all the basic features of this algorithm. The calculation proceeds by stacking particles representing the model core on an octal tree that is used to identify the nearest neighbours of a particle(exactly 50), and to calculate the magnitude of net force experienced by a particle. The standard cubic-spline kernel is employed to smooth out contributions from individual particles and the smoothing length of each particle is adjusted such that it has exactly 50 neighbours. The prescription of the Riemann artificial viscosity with the viscosity parameter, α = 0.5, suggested by Monaghan (1997) was preferred in this set of simulations.

2.3. Calculating dust and gas temperature

Temperature of individual gas particles in this work was calculated by explicitly solving the respective equations of thermal balance for the gas and dust. Dust in the cores is primarily heated by the interstellar radiation field (ISRF), characterised in the present work by a heating rate, $\Gamma _{{\text{ISRF}}}$ . Other functions that characterise the thermal energy budget of a dust grain are its cooling function, Λ_dust, and gas–dust coupling function, $\Lambda _{\text{g}-\text{d}}$ , that accounts for energy transfer during collisional interaction between a dust grain and a gas molecule. $\Lambda _{\text{g}-\text{d}}$ , is effectively a cooling function when dust is cooler than gas as it gains thermal energy from warmer gas molecules. On the contrary, $\Lambda _{\text{g}-\text{d}}$ , behaves like a heating function when dust is warmer than gas. In the present work, we use the standard heating/cooling rates suggested by Goldsmith (2001) for typical clouds in the Solar neighbourhood, and where applicable, those corrected elsewhere in literature. The dust heating rate due to the ISRF:

(6) $$\begin{equation} \Gamma _{{\rm ISRF}} = 3.9\times 10^{-24}\Big [\frac{n(H_{2})}{\mathrm{cm}^{-3}}\Big ]\chi \ \mathrm{erg}\mathrm{cm}^{-3}\mathrm{s}^{-1}, \end{equation}$$

where, χ ~ 10⁻⁴, is the factor by which the ISRF is attenuated in a putative star-forming clump, whilst n(H ₂) is the molecular hydrogen number density. We hold this magnitude of χ fixed over the entire duration of a simulation. This quantity is reasonably well constrained in literature; Valdivia et al. (2016) showed that χ ~ 10⁻⁴ for typical prestellar densities ( ≳ 10⁴ cm⁻³) which is also consistent with the magnitude reported for IRDCs ¹ (Goldsmith 2001).

Next,

(7) $$\begin{equation} \Lambda _{\text{g}-\text{d}} = 2\times 10^{-33}\Big [\frac{n(H_{2})}{\mathrm{cm}^{-3}}\Big ]^{2}\Big (\frac{\triangle T}{\text{K}}\Big )\Big (\frac{T_{\text{gas}}}{10\text{ K}}\Big )^{0.5} \ \mathrm{erg}\mathrm{cm}^{-3}\mathrm{s}^{-1} \end{equation}$$

for standard values of dust-to-gas ratio and dust grain size (Burke & Hollenbach 1983); $\triangle T = T_{\text{dust}} - T_{\text{gas}}$ , $T_{\text{dust}}$ , and $T_{\text{gas}}$ being respectively the temperature of dust, and gas. In one of the test cases, listed 8 in Table 1, we raised the dust grain size by about a factor of eight to mimic the effect of the so-called fluffy dust grains. Fluffy dust grains enhance the contribution due to $\Lambda _{\text{g}-\text{d}}$ and assist gas cooling when the dust is cold. This is usually true for cores in typical low-mass star-forming clouds. In fact, we then repeated calculations for all the other test cores(listed 9–12 in Table 1) to investigate if the puffed-up dust grains also lowered the gas temperature in other starless cores before eventually leading them to collapse. Finally, the dust cooling function:

(8) $$\begin{equation} \Lambda _{\text{dust}} = 4.22\times 10^{-31}\Big (\frac{T_{\text{dust}}}{\text{K}}\Big )^{6}\Big [\frac{n(H_{2})}{\mathrm{cm}^{-3}}\Big ] \ \mathrm{erg}\ \mathrm{cm}^{-3}\ \mathrm{s}^{-1} \end{equation}$$

Table 1. Physical details of test cores modelled in this work; listed at the bottom of the table are the respective references from where physical details of these cores have been obtained.

^a Q $=\frac{E_{\text{g}}}{E_{\text{p}}}$ , ratio of gravitational energy with energy due to confinement by external pressure.

^b Free-fall time for the initial configuration; $t_{\text{ff}}=(\frac{3\pi }{32G\bar{\rho }})^{1/2}, \bar{\rho }$ —average initial gas density in a test core.

^c Alves, Lada, & Lada (2001); Bergin et al. (2006).

^d Harvey et al. (2003); Chitsazzadeh (2014).

^e Tafalla et al. (2004); Kirk et al. (2006).

^f Chitsazzadeh et al. (2014).

^g Chitsazzadeh (2014).

(Bate & Keto 2015). Without an elaborate chemical network at our disposal cooling of molecular gas was implemented by adopting a parametrised cooling function, $\Lambda _{\text{gas}}$ , for undepleted molecular abundances in clumps,

(9) $$\begin{equation} \Lambda _{\text{gas}} = A\Big (\frac{T_{\text{gas}}}{10\text{ K}}\Big )^{B} \end{equation}$$

Goldsmith (2001); this cooling function mimics very well the observed cooling-rate in clumps within the Solar neighbourhood. The coefficients (A, B) here are the same as (α, β) used by Goldsmith (2001) and listed in his Table 1. We avoid (α, β) here for fear of confusing them with the SPH artificial viscosity coefficients. We assume that gas in a test core is heated externally only due to cosmic-rays(CRs), for it is reasonable to neglect photoelectric heating since cores are usually found in well-shielded regions of a MC. The heating rate, Γ_CR, due to CRs is

(10) $$\begin{equation} \Gamma _{{\rm CR}} = 10^{-27}\Big [\frac{n(H_{2})}{\mathrm{cm}^{-3}}\Big ] \ \mathrm{erg}\ \mathrm{cm}^{-3}\ \mathrm{s}^{-1}. \end{equation}$$

Finally, the equilibrium dust and gas temperature was calculated by solving simultaneously the respective equations of thermal balance:

(11) $$\begin{equation} \Gamma _{{\text{ISRF}}} - \Lambda _{\text{dust}} + \Lambda _{\text{g}-\text{d}} & =& 0, \\ \end{equation}$$

(12) $$\begin{equation} \Gamma _{\text{CR}} - \Lambda _{\text{gas}} - \Lambda _{\text{g}-\text{d}}& =& 0. \end{equation}$$

The equilibrium gas temperature, $T_{\text{gas}}^{\text{eq}}$ , of a particle, p, corresponds to an equilibrium energy, $u_{\text{p}}^{\text{eq}}$ , that is used to define the thermal timescale, $\text{d}t_{\text{thermal}}$ , for an SPH gas particle as

(13) $$\begin{equation} dt_{\text{thermal}} = \vert u_{\text{p}} - u_{\text{p}}^{\text{eq}}\vert \Big (\frac{\text{d}u}{\text{d}t}\Big )_{\text{net}}^{-1}, \end{equation}$$

where $(\frac{\text{d}u}{\text{d}t})_{\text{net}}$ is the net heating/cooling that the particle is subject to. It is the equivalent of (n ²Λ − nΓ); Λ and Γ being respectively the effective cooling, and heating rate for an SPH particle. The internal energy of a particle was then revised according to

(14) $$\begin{equation} u_{\text{p}}^{\prime } = u_{\text{p}}^{\text{eq}} + (u_{\text{p}} - u_{\text{p}}^{\text{eq}})\exp \Big (\frac{-\text{d}t}{\text{d}t_{\text{thermal}}}\Big ), \end{equation}$$

where $\text{d}t$ is the real time step for the particle p. In case a particle is heating/cooling in quasi-static manner, $\text{d}t_{\text{thermal}}\gtrsim \text{d}t$ , so that the above equation reduces to

(15) $$\begin{equation} u_{\text{p}}^{\prime } = u_{\text{p}} - \text{d}t\Big (\frac{\text{d}u}{\text{d}t}\Big )_{\text{net}}. \end{equation}$$

On the other hand, if it cools/heats rapidly, $\text{d}t_{\text{thermal}}\ll \text{d}t$ , so that it quickly attains its equilibrium energy and $u_{\text{p}}^{\prime } = u_{\text{p}}^{\text{eq}}$ . Temperature, T _ext (see Table 1), of the ICM particles confining a model core was held constant throughout the length of a realisation. We must note that our extant interest is confined only to the low-mass star-forming clouds in the local Universe. The question about the impact of the ambient environment on the evolution of cores, though pertinent, is beyond the scope of this work. It must also be added that the choice of the respective heating functions Γ_CR, Γ_ISRF, and the molecular-line cooling function, $\Lambda _{\text{gas}}$ , used in the calculations presented in this work is appropriate, for the external heating is significantly stronger in regions of high-mass star-formation and at high redshifts (see e.g. Redman et al. 2004). Also, the details of molecular cooling are likely to be different in these respective environs.

2.4. Initial conditions

The choice of initial conditions is relatively simple and consists of a core in approximate (thermal)pressure equilibrium with its confining medium. All test cores were initially Jeans stable. As discussed in Section 2.1, we modelled each of our test core as an isothermal pressure-confined BES( $\xi _{\text{b}}$ = 3), with the exception of the core B68 that was also modelled as a critically stable BES( $\xi _{\text{b}}$ = 6.5) and a uniform density sphere to test if the specific choice of the initial density distribution possibly bears upon the evolution of the core. The magnitude of temperature, $T_{\text{gas}}$ , initially assigned to the gas in a model core was obtained from literature for the respective core. The initial radius of the core, R _core, listed in column 4 of Table 1 was calculated using T _gas in Equation (4) above that ensures the initial setup is in pressure equilibrium. The analogous expression for a core assembled as a sphere of uniform density is

$$\begin{equation*} T_{{\rm gas}} = \frac{\bar{m}G M_{{\rm core}}}{3k_{\text{B}}R_{{\rm core}}}, \end{equation*}$$

where symbols have their usual meaning. As noted earlier, particles representing the externally confining ICM were always maintained at a fixed temperature, T _ext, which was calculated using Equation (5) for the respective BESs whilst for the Uniform density sphere, it was calculated by assuming the same centre-to-boundary density contrast as that for a pressure-confined BES( $\xi _{\text{b}}$ = 3). Listed in Column 7 of Table 1 is the ratio of the energy due to self-gravity against that due to the externally confining medium, $Q_{\text{init}}$ ; a magnitude smaller than unity implies a pressure-confined configuration. Each of the three polytropes were assembled by randomly positioning particles within them. The particle distribution in each polytrope was relaxed by allowing it to evolve isothermally for a few sound-crossing times. The resulting density distribution for the respective BESs has been shown on the two panels of Figure 3. The settled polytropes were then stretched to the desired core radius, $R_{\text{core}}$ . The model core in this work was assumed to be composed of the usual cosmic mixture.

Figure 3. Density profile of the settled stable Bonnor– Ebert sphere ( $\xi _{\text{b}}$ = 3), and the critically stable Bonnor– Ebert sphere ( $\xi _{\text{b}}$ = 6.45) has been shown respectively, on the left- and the right-hand panel. Dots represent the density of gas particles assembled in either polytrope whilst the analytically expected density distribution given by Equation (3) has been shown by a continuous line in these respective plots.

Resolution. The number of particles, $N_{\text{gas}}$ , to be assembled in a polytrope was calculated using the equation:

(16) $$\begin{equation} N_{{\rm gas}} = \Big (\frac{3}{32\pi }\Big )\cdot \Big (\frac{N_{{\rm neibs}}}{h_{{\rm avg}}^{3}}\Big )R_{{\rm core}}^{3}, \end{equation}$$

N _neibs = 50, being the number of neighbours of each particle. The average smoothing length, h _avg, that determines the spatial extent of the smallest resolvable region in a realisation was calculated such that it satisfied the Truelove criterion (Truelove et al. 1997), of resolving the Jeans instability. If $\lambda _{\text{J}}$ is the Jeans length at the minimum gas temperature to be attained ( ~ 7 K) in these realisations for the typical protostellar density, ~ 10⁷ cm⁻³, then the Truelove criterion demands, $h_{{\rm avg}}\lesssim \frac{\lambda _{\text{J}}}{4}\sim 281$ AU, which ensures there is no artificial fragmentation at the desired density. Furthermore, this choice of spatial resolution transforms into a much higher resolution for the rest of the core since $\frac{R_{{\rm core}}}{h_{{\rm avg}}}\gg$ 30 for each realisation. Consequently, these simulations not only satisfy the SPH equivalent of the Truelove criterion (Hubber, Goodwin, & Whitworth 2006), but also the more stringent criterion suggested by Federrath et al. (2014) that ensures convergence in the calculation of energy and momentum within a test core. Simulations reported here were developed with ~ 0.7 M particles out of which ~ 0.3 M represented the ICM. Each realisation required about 1 300 CPU hours with eight threads.

Apart from these, $N_{\text{icm}}$ number of particles representing the externally confining ICM and always maintained at temperature, $T_{\text{ext}}$ , were assembled in a jacket of thickness of 2 $h_{\text{avg}}$ around the core. The number of ICM particles was calculated such that the number density of particles on either side of the gas-ICM interface was identical. The core+ICM assembly was then enveloped by a layer of dead particles, the so-called boundary particles, to prevent gas/ICM particles from escaping. The initially static ICM particles exert only the hydrodynamic force on gas particles whilst the boundary particles do not contribute to the net force and always remain static in the computational domain.

Representing the protostar in case the model core became singular

The protostar in a realisation was represented by a sink particle, the density threshold, $n_{\text{sink}}$ , for which was set at $n_{\text{sink}}\sim \ 10^{7}$ cm⁻³. In other words, an SPH particle with density exceeding the threshold, $n_{\text{sink}}$ , was replaced with a sink particle if in addition it also had a negative divergence of velocity and acceleration (e.g. Bate, Bonnell, & Price 1995; Bate & Burkert 1997; Federrath et al. 2010). The radius of a sink particle was set at 2.5 times the average SPH smoothing length of the candidate sink particle and turns out to be on the order of ~ 800 AU. Gas particles that enter the sink radius were accreted by it.

3 RESULTS

3.1. Starless cores investigated in this work

In order to justify the choice of initial conditions listed in Table 1, we begin by briefly examining the physical properties of these cores documented in contemporary literature. This endeavour will also serve us well whilst testing the validity of results deduced from our realisations.

B68. One of the most well-studied cores and is reported to have a density profile that matches the density distribution of a unstable BES of radius, $\xi _{\text{b}}$ = 6.9 and radius, $r_{\text{c}}\sim$ 0.05 pc (e.g. Alves et al. 2001), though there is a slight variation in the estimate of its radius (between ~ 0.035 and 0.05 pc), depending on the assumed distance to the core (e.g. Nielbock et al. 2012; Schnee et al. 2013).

L694-2. First classified as starless by Lee & Myers (1999) and Lee, Myers, & Tafalla (1999), using data from the SDSS survey. Estimates of its physical parameters vary and a mass between 0.5 M_⊙ (Levhakov et al. 2013, by studying NH₃ emission) and 3.1 M_⊙ (Harvey et al. 2003) has been suggested, subject to the tracer used to study it. Similarly, the core has been suggested to have a radius between 0.05 pc (Williams, Lee, & Myers 2006) and 0.15 pc (Harvey et al. 2003). In a more recent contribution, Chitsazzadeh (2014) suggested a core mass of ~ 2.5 M_⊙ and a radius ~ 0.1 pc and furthermore, their radiative transfer modelling of L694-2 showed that its cold interior was cocooned by a warmer envelope. They estimated the gas temperature to be ~ 19 K in the envelope of this core and ~ 7 K towards its centre.

L1517B. A starless core that is probably experiencing a breathing mode, i.e. exhibiting signs of in-fall close to the centre, but has an expanding envelope (Keto & Field 2005; Fu, Gao, & Lou 2011). As with other cores, estimates of its physical properties are subject to the molecular tracer. For instance, the mass and radius for this core varies between ~ 0.05 M_⊙ and ~ 0.06 pc(Benson & Myers 1989) to ~ 1.7 M_⊙(Kirk et al. 2005) and ~ 3.89 M_⊙(Fu et al. 2011). The gas temperature for this core has been reported to be ~ 9 K (Tafalla et al. 2002; Keto & Field 2005).

L1689-SMM16. A unusually massive starless core and supposedly has mass well in excess of its (thermal)Jeans mass, $M_{\text{Jeans}}$ , and therefore has been described as Super-Jeans (Sadavoy et al. 2010a, 2010b). However, signs of outwardly directed gas flow have been detected in its envelope, but is believed to be on the verge of collapse given the relatively large, $\frac{M_{\text{core}}}{M_{\text{Jeans}}}$ , ratio (Chitsazzadeh et al. 2014).

L1521F. Probably associated with a weak, poorly collimated outflow and is believed to have formed its first hydrostatic core. Such objects have a relatively small luminosity, on the order of ≲ 0.1 L_⊙, and are known as Very low luminosity objects(VeLLOs) in contemporary literature (e.g. Di Francesco et al. 2007). Some of the earliest detections of such objects were reported by Young et al. (2004) (e.g. L1014-IRS) and Kauffmann et al. (2005) (e.g. L1148-IRS). The central condensation in L1521F is sub-stellar and has a bolometric luminosity of ~ 0.05 L_⊙ (e.g. Bourke et al. 2006; Terebey et al. 2009), and estimates of mass for this core vary between 3 M_⊙ (Onishi, Mizuno, & Fukui 1999) and 5.5 M_⊙ (Crapsi et al. 2004). In the present work, we adopt the more recently deduced physical properties for this core (Chitsazzadeh 2014). Radiative transfer models developed by this author for the L1521F constrained the temperature of gas in the envelope to between 18 and 21 K.

3.2. Evolution of gas density

The polytrope in each realisation became centrally concentrated due to an inwardly propagating compressional wave triggered by the gradual cooling of gas in it. In the set of realisations discussed here, the pressure exerted by the confining ICM was initially purely thermal in nature since the particles resembling the ICM were initially static. The initial magnitudes of P _ext are therefore likely to be somewhat smaller, perhaps by a factor of a few, in comparison with the magnitude of pressure exerted by the ICM in typical star-forming clouds. Observationally, however, it is known that a non-thermal component also makes a significant contribution to the pressure, $P_{\text{ext}}$ (e.g. Bailey, Banerjee, & Caselli 2014), in these clouds. The somewhat lower initial magnitude of $P_{\text{ext}}$ is unlikely to affect our calculations here since the particles representing the ICM in these realisations are mobile during a calculation and therefore additionally contribute to $P_{\text{ext}}$ over the length of a calculation.

Starless cores : B68, L694-2, L1517B, and L1689-SMM16. The picture in Figure 4 is a representative plot showing the rendered density image of the cross-section of the mid-plane of the B68 modelled as a pressure-confined BES ( $\xi _{\text{b}}=3$ ), and listed 1 in Table 1. Having acquired its peak density, the core rebounded. This is evident from the picture on the central-panel of Figure 4 that shows the epoch at which this core acquired its peak density and evidently, the inwardly directed compressional wave appears to have given way to expansion. Shown on the top-panel of Figure 5 is the radial distribution of gas density within this core at different epochs. The model core is centrally most concentrated at t ~ 0.32 Myr when in fact, its density profile is similar, or rather steeper than that of a unstable BES having radius, $\xi _{\text{b}}$ = 7. Also, at this epoch, the radial density falls off by more than an order of magnitude over a radius, $r_{\text{c}}\sim$ 0.05 pc. We note that the core appears to be demonstrating breathing modes as it acquired a second density peak at t ~ 1 Myr when calculations for this realisation were terminated. The core B68 was also modelled as a critically stable BES ( $\xi _{\text{b}}\sim$ 6.5), and a sphere having uniform density. However, contrary to the previous realisation, the test core modelled as a critically stable BES initially showed expansionary tendency as is reflected by the initial fall in its central density. Then, as the core continued to cool, and accompanied with a compressional wave, it became centrally peaked (t ~ 1 Myr), with a density profile similar to that of a BES having radius $\xi _{\text{b}}$ = 7. Thereafter, this test core rebounded as is evident from the plots shown on the central panel of Figure 5. The uniform density sphere on the other hand, contracted rapidly and acquired a density peak at t ~ 0.16 Myr; see lower panel of Figure 5. Qualitatively, the results from these latter two realisations converge with those deduced for the core modelled as a pressure-confined BES( $\xi _{\text{b}}$ = 3). As with the core B68, the other remaining cores viz., L694-2, L1517B, and L1689-SMM16 also became centrally concentrated as an inwardly propagating compressional wave swept them up (see also Keto et al. 2015). Shown in Figure 6 are plots of the radial density distribution of gas in these latter cores at different epochs of their evolution. Observe that having acquired a centrally concentrated density distribution when the gas density in each core core mimicked the density profile of a unstable BES, the cores rebounded. The radius of each core at this epoch, i.e. the distance from the centre over which density falls over an order of magnitude, is consistent with various reported observations listed above. There is also evidence for radial oscillations in each of these realisations.

Figure 4. Rendered density plots for the model B68 core, listed 1 in Table 1, show the assembly of a centrally condensed core due to an inwardly directed compressional wave. Epochs corresponding to the images on respective panels of this figure are t = 0.001 Myr, 0.32 Myr, and 1 Myr; arrows overlaid on these plots represent the direction of the underlying gas flow.

Figure 5. Shown on the upper, central, and the lower panel are plots of the radial distribution of gas density at different epochs of the core B68 modelled as respectively, the stable Bonnor– Ebert sphere( $\xi _{\text{b}}$ = 3), the critically stable Bonnor– Ebert sphere( $\xi _{\text{b}}$ = 6.45) and a uniform density sphere. In each of these three realisations, the model core acquires a centrally condensed form when its density profile mimics that of a unstable Bonnor– Ebert sphere plotted with a dashed-black line, though the epoch at which this happens depends on the initial model of the core (see text).

Figure 6. Shown on the upper, central, and the lower panel are the plots showing the radial distribution of gas density at different epochs for respectively the cores L694-2, L1517B, and L16189-SMM16. Note that each one them acquired a centrally condensed form whence their respective density profile mimicked that of a unstable BES. The epoch when each core acquires its peak density and the radius, $\xi _{\text{b}}$ , of the BES that fits its density distribution is of course different.

The core L1521F. The realisation listed 7 in Table 1 for this core was developed by adopting the standard dust grain size in the gas–dust coupling function, $\Lambda _{\text{g}-\text{d}}$ , given by Equation (7). The calculation was then repeated with dust grain size an order of magnitude higher (listed 8 in Table 1). This effectively raised the contribution to gas cooling due to the collisional coupling between gas and dust. The core in the former realisation did not collapse and like the starless cores discussed previously, simply acquired a centrally peaked distribution before rebounding. In the latter case, however, the core did in fact become singular and formed a sub-stellar object represented by a sink particle. The temporal variation of the radial distribution of gas density in either realisation of this core is shown on respectively the upper and central panel of Figure 7. Evidently, the higher dust grain size proved efficient towards lowering the temperature of gas within the core which in turn assisted in-fall. Finally, shown on the lower panel of this figure is the radial density profile for other test cores (realisations 9–12 in Table 1), but now with gas temperature calculated with puffed-up dust grains. The resulting density distribution for these respective cores is similar to that of the VeLLO L1521F. It is therefore clear that irrespective of the physical parameters of an individual core, a larger size of dust grains leading to a higher contribution due to gas–dust coupling assists core collapse. We will expand upon this further in Section 3.4 below where the temperature profile deduced for these cores will be presented.

Figure 7. Shown on the upper panel of this plot is the radial distribution of gas density at different epochs of the core L1521F modelled as a stable Bonnor– Ebert sphere( $\xi _{\text{b}}$ = 3), and listed 7 in Table 1. The density distribution for the core that did in fact collapse to become a VeLLO, and listed 8 in Table 1, is shown on the central panel. Shown on the lower panel is the radial density distribution for other cores listed 9–12 in Table 1 at the time respective calculations were terminated; see text for details.

Let us now examine the accretion history and the mass of the VeLLO L1521F. We remind, the central object in this test core was represented by a sink particle with a density threshold, ~ 10⁷ cm⁻³. This density threshold is at least an order of magnitude higher than the central density ( ~ 10⁶ cm⁻³) reported by Onishi et al. (1999) and Crapsi et al. (2004) for this core. Shown on the left-hand panel of Figure 8 is the rate, $\dot{M}_{{\rm acc}}\equiv \frac{\text{d}M_{*}}{\text{d}t}$ , at which the sink particle in this core accretes the in-falling gas. After an initially sluggish rate of accretion, it steadily gains momentum reaching a peak rate of ~ 10⁻⁶ M_⊙ yr⁻¹, before eventually petering off over a period of ~ 2 × 10⁵ yr. This observed evolution of $\dot{M}_{{\rm acc}}$ is far from uniform and in fact, is qualitatively consistent with the one derived by Whitworth & Ward-Thompson (2001) for a core modelled with Plummer-like density profile. Importantly, the magnitude of $\dot{M}_{\text{acc}}$ observed here is consistent with the observationally deduced magnitude by Takahashi, Ohashi, & Bourke (2013). Using this observed rate of accretion, we can calculate the approximate mass of the central accreting object, M _*, with the expression for Bondi–Hoyle accretion:

(17) $$\begin{equation} \dot{M}_{{\rm acc}}=\frac{\bar{\rho }G^{2}M_{*}^{2}}{v_{\text{r}}^{3}}, \end{equation}$$

(Bondi 1952), where $\bar{\rho }\sim 10^{-17}$ g cm⁻³, is the average density of in-falling gas, $v_{\text{r}}\sim$ 0.1 km s⁻¹ and, G, the gravitational constant; other symbols have their usual meanings. The accretion-rate defined by Equation (17) gives a upper limit to the magnitude of $\dot{M}_{{\rm acc}}$ ; by contrast, the one in the Shu model of the collapse of a singular isothermal sphere, is constant (Anathpindika 2011). For an average accretion rate, $\dot{M}_{{\rm acc}}\sim \ 2\times 10^{-7}$ M_⊙ yr⁻¹, Equation (17) yields, M _* ~ 0.04 M_⊙, a sub-stellar object, which is roughly consistent with the mass accreted by the sink-particle at the centre of this core as is visible from the plot shown on the right-hand panel of Figure 8. Now, although the sink particle is likely to continue acquiring mass via accretion, the observed relatively low rate of accretion makes it likely that the object in this case will remain sub-stellar which is consistent with the fact that it is a VeLLO.

Figure 8. The plot on the left-hand panel shows the rate of mass-accretion of the sink particle that forms at the centre of the core whose density profile was shown on the central panel of Figure 7. The accretion rate peaks at ~ 10⁻⁶ M_⊙ yr⁻¹ before petering off. Shown on the right-hand panel of this figure is the total mass accreted by this sink particle.

3.3. Velocity field within contracting cores

In the interest of brevity, we will restrict the present discussion only to the two contrasting cases viz., the starless core B68 and the VeLLO, L1521F; the behaviour of the remaining starless cores was similar to that of the core B68. Shown on the upper panel of Figure 9 is the radial component of gas velocity in this core. The negative magnitude of velocity on this plot as usual denotes inwardly directed gas flow. Evidently, the onset of a strong compressional wave pushed gas inward that enhanced the density near the centre of the model core. Although this compressional wave can be seen to be the strongest just before the core achieved its central density peaks (t ~ 0.20 Myr and t ~ 0.76 Myr), the gas flow within the core remained subsonic at all times. The magnitude of inward/outward gas velocity is between ~ 0.02 and 0.15 km s⁻¹ which is consistent with that inferred by for instance, Maret, Bergin, & Lada (2007) through observations of a number of emission lines towards B68 and that suggested by radiative-transfer modelling of the core (e.g. Keto et al. 2006). At early times, t ~ 0.2 Myr, gas within the interiors of the core appears to be outwardly directed before being swept-up by the compressional wave.

Figure 9. Upper-panel: As in the upper-panel of Figure 5, but now the variation of the radial component of gas velocity in the core B68. Central-panel: Radial component of the velocity of in-falling gas at different epochs of the evolution of the core L1521F(VeLLO), listed 8 in Table 1. Lower-panel: Radial variation of gas-velocity dispersion at different epochs of the core B68 listed 1 in Table 1; see text for details.

Next, shown on the central panel of this figure is the radial velocity of the in falling gas in the VeLLO L1521F, the magnitude of which is consistent with that reported by Onishi et al. (1999). We note that the velocity of the in falling gas in the envelope of the core in this case is transonic or even mildly supersonic in contrast to the gas-velocity observed for starless cores such as the B68. This inward/outward gas motion in the core generates a velocity dispersion, the magnitude of which at different radii within the core B68 has been shown on the lower-panel of Figure 9. Observe that the magnitude of velocity dispersion is relatively large close to the centre of the test core and its magnitude steadily increases as the core continues to evolve. Nevertheless, the velocity field generated appears sufficient to arrest the collapse of the core. We will revisit this point in Section 3.5 with reference to the BE stability criterion.

3.4. Temperature profile within test cores

Shown on the left-hand panel of Figure 10 is the radial distribution of gas temperature calculated using the standard grain size in Equation (7) within the test cores at the epoch when the respective cores acquired their first density peak. The radial temperature plots shown here were constructed by dividing the core into concentric shells followed by the calculation of the density averaged temperature for each shell. From this plot, it is readily visible that although gas in the outer envelope of each of these cores remained relatively warm and approximately isothermal, gas temperature close to their respective centre exhibited a sharp fall. It is only in this very small radius about the centre of the core that there is any evidence of significant cooling. Nevertheless, it proved insufficient to sustain collapse. However, the temperature profile and the magnitude of temperature deduced here for cores known to be starless viz., B68, L694, L1517B, and L1689-SMM2, is consistent with that reported for these cores following several detailed observational studies discussed in Section 3.1 above. The core L152F, known to be a VeLLO, remained starless in our realisation (listed 7 in Table 1); where the standard dust grain size was used in these calculations.

Figure 10. Left-panel: As in Figure 9, but now the radial variation of gas temperature within respective test cores that remained starless at the epoch when they acquired their peak density; see text in Section 3.4 . Right-panel: Radial distribution of gas temperature within respective cores with puffed-up dust grains. Attention is especially drawn to the fact that in these latter realisations the cores developed a more extended cold region close to their centre fairly early in their evolutionary cycle.

Shown on the right-hand panel of Figure 10 is the radial distribution of gas temperature within these cores early in their respective evolutionary cycle, but with the temperature now calculated using the puffed-up dust grains; see Section 2.3. In this case, not only could we induce the core L1521F to collapse (see central-panel of Figure 7), but all the other cores that previously remained starless also collapsed (see lower panel of Figure 7). The remarkable difference between the plots on the respective panels of Figure 10 probably hold the key to the question about the stability of a core, and its propensity to become protostellar. In this case, gas closer to the centre of a core and therefore more dense than that in its envelope, appears to have cooled much more uniformly. This likely assisted the in-fall.

Our simulations have shown that not all starless cores remain thermally sub-critical, i.e. n_c ≲ 10⁵ cm⁻³. It appears that whilst a thermally super-critical core could be more likely to become protostellar, though it need not necessarily become protostellar. Consequently, a competing cooling mechanism must assist the gas to cool as a core continues to contract. This can be easily demonstrated with a simple deduction that follows. A core will become protostellar only if gas in it satisfies the Jeans criterion at all radii (Anathpindika & Di Francesco 2013). Thus, $M(r) > M_{\text{Jeans}}(r)$ , i.e.

(18) $$\begin{equation} \frac{4}{3}\Big (\frac{G^{3}}{\pi }\Big )^{1/2}r^{3}\rho ^{3/2}(r) > a^{3}(r), \end{equation}$$

where symbols have their usual meanings. Rearranging this equation yields

$$\begin{equation*} \rho ^{3/2}(r) > \frac{3}{4}\Big (\frac{\pi }{G^{3}}\Big )^{1/2}\frac{a^{3}}{r^{3}}. \end{equation*}$$

This equation yields the condition that gas density in a core must satisfy in order to become protostellar—

$$\begin{equation*} \frac{\text{d}\rho }{\text{d}r} > \frac{2\rho }{a}\Big (\frac{\text{d}a}{r} - \frac{3a}{r}\Big ), \end{equation*}$$

or equivalently, the sound-speed in the core must vary as

$$\begin{equation*} \frac{1}{a}\frac{\text{d}a}{\text{d}r} < \frac{3}{r} + \frac{1}{2\rho }\frac{\text{d}\rho }{\text{d}r}. \end{equation*}$$

A simple integration of the above expression gives us the expression for sound speed in a core that is likely to collapse—

(19) $$\begin{equation} a(r) < \mathcal {C}r^{3}\rho ^{1/2}, \end{equation}$$

$\mathcal {C}$ , being the constant of integration. If we denote the subscripts ‘1’ and ‘2’ to denote the gas temperature on respectively the inward(towards core centre) and outward (away from core centre) side of r, then, a(r) ≡ a ₁ − a ₂, moving out from the centre of the core. From Equation (19), we have

$$\begin{equation*} a_{1} < a_{2} + \mathcal {C}r^{2}\rho ^{1/2}, \end{equation*}$$

or after plugging the constant of integration

(20) $$\begin{equation} a_{1} < a_{2} + G^{1/2}\frac{r^{3}}{R_{\text{core}}^{2}}\rho ^{1/2}(r). \end{equation}$$

In other words, the sound speed(or equivalently the temperature), looking inward to the centre of a core should decrease as defined by Equation (20) for collapse to ensue.

3.5. Bonnor–Ebert(BE) stability criterion

Next, we undertake a comparison of the mass of the core, $M_{\text{core}}$ , with its critical BE mass, M _BE, given by Equation (21) below, to ascertain if the BE mass is a reliable indicator of the virial boundedness of a core:

(21) $$\begin{equation} M_{{\rm BE}} = \frac{a_{0}^{4}}{(4\pi G^{3})^{1/2}P_{{\rm ext}}^{1/2}}\cdot \mu (\xi _{\text{b}})\mathrm{e}^{-\psi (\xi _{\text{b}})}, \end{equation}$$

$\mu (\xi )=\xi ^{2}\frac{\text{d}\psi }{\text{d}\xi }$ . For each of our test cores that remained starless (see Figures 5–7), we used the respective magnitude of $\xi _{\text{b}}$ , the radius of the BES whose density profile fitted the observed density distribution of our model cores, to calculate M _BE. For instance, in the specific case of the core B68, the above equation simply becomes

(22) $$\begin{equation} M_{{\rm BE}} = 1.125\frac{a_{0}^{4}}{G^{3/2}P_{{\rm ext}}^{1/2}}. \end{equation}$$

The upshot of plugging in the appropriate values of gas temperature, ~ 15 K for the core B68 (see left-hand panel of Figure 10; t ~ 0.32 Myr), and $\frac{P_{{\rm ext}}}{k_{\text{B}}}\sim 4.5\times 10^{4}$ K cm⁻³, after including velocity dispersion of the ICM, is, M _BE ~ 2.66 M_⊙; $\frac{M_{{\rm core}}}{M_{{\rm BE}}} <$ 1. In other words, the core cannot possibly become singular, which in fact, is consistent with the observed fate of the model core.

Similarly, plugging-in the sound speed corresponding to the gas temperature ( ~ 16, ~ 13, and ~ 16 K for respectively the cores L694-2, L1517B, and L1689-SMM2), along with the radius of the BES whose profile fitted the extant density distribution for the individual core (see plots in Figure 6), we find that $\frac{M_{{\rm core}}}{M_{{\rm BE}}} >$ 1 for each of the three cores. Evidently, the comparison between the BE mass and the core mass yields contradicting inferences as regards the virial stability of a core, for the unstable BES, by definition, is susceptible to collapse. However, as discussed above, such a collapse did not materialise in any of these cores leading us to conclude, the BE stability criterion is a unreliable predictor of the future evolution of a test core. The problem with this analysis is alleviated when the sound speed in Equation (21) is corrected to include the velocity dispersion (see e.g. Johnstone et al. 2000). Our inference is consistent with the recent findings reported by Pattle et al. (2015) where a number of cores surveyed in the Ophiuchus as part of the JCMT Gould Belt survey ² , though unstable according to the BE stability criterion, did not show any credible signs of star formation.

3.6. The virial chart for cores

Finally, we discuss the virial state of our test cores over the course of their evolution. For a purely hydrodynamic calculation, the virial equation is of the form—

(23) $$\begin{equation} \frac{1}{2}\ddot{\mathrm{I}} = 2E_{\text{k}} + E_{\text{g}} + E_{\text{p}}, \end{equation}$$

$E_{\text{k}}, E_{\text{g}}$ , and $E_{\text{p}}$ being respectively the kinetic energy, gravitational energy, and the energy due to the pressure confining the core. The second moment of inertia, $\ddot{\mathrm{I}}$ , is greater than zero for a dispersing core, less than zero for one that is in-falling, and equal to zero for a core in virial equilibrium. The virial components for the cores tested in this work were calculated as follows:

(24) $$\begin{equation} E_{\text{k}} = \frac{3}{2}\Sigma _{\text{i}}m_{\text{i}}\sigma _{\text{i}}^{2}, \end{equation}$$

where the summation extends over all the gas particles and $\sigma _{\text{i}} = (a_{\text{i}} + (\sigma _{\text{v}})_{\text{i}})$ ; $a_{\text{i}}$ , being the sound speed for a particle, $(\sigma _{\text{v}})_{\text{i}}$ , its velocity dispersion and $m_{\text{i}}$ , the mass of individual particle. The gravitational energy, $E_{\text{g}}$ , as

(25) $$\begin{equation} E_{\text{g}} = -GM_{{\rm core}}\Sigma _{i}\frac{m_{i}}{\vert \mathbf {r}_{\text{c}} - \mathbf {r}_{\text{i}}\vert }. \end{equation}$$

As in Equation (24), the summation here also runs over all gas particles; $\mathbf {r}_{\text{c}}$ , is coordinate of the centre of the core, and $\mathbf {r}_{\text{i}}$ , the position of an individual particle in the core. Finally, the energy corresponding to the magnitude of external pressure acting on the core:

(26) $$\begin{equation} E_{\text{p}} = P_{{\rm ext}}V_{{\rm core}}\equiv k_{\text{B}}V_{{\rm core}}\bar{n}_{\text{icm}}T_{\text{icm}}, \end{equation}$$

where V _core is the volume of the test core and $\bar{n}_{\text{icm}}$ , the mean particle density of the externally confining medium. We now define two quantities, $Q\equiv \frac{E_{\text{g}}}{E_{\text{p}}}$ , that quantifies the strength of self-gravity relative to the external pressure, and $X\equiv \frac{-(E_{\text{g}} + E_{\text{p}})}{2E_{\text{k}}}$ , the virial ratio. Note that the energy, $E_{\text{p}}$ , exerted by the externally confining medium, like the gravitational energy, $E_{\text{g}}$ is also taken as negative since it is inwardly directed.

A core is defined to be virially bound when X ⩾ 1, and pressure-confined, when Q < 1 (e.g. Pattle et al. 2015). Shown in Figure 11 is the variation in the magnitude of these quantities at different epochs of our test cores; note that the virial chart for B68 shown here is that for the core modelled as a stable BES( $\xi _{\text{b}}$ = 3). As is visible from the respective charts for our model cores, each of these cores started as a virially bound and pressure-confined entity(t ~ 0.001 Myr for each). As the cores continued to contract and acquire a centrally concentrated density profile, then on the one hand, the extent of gravitational boundedness in each of these cores can be seen to be increasing steadily, but on the other, they tended to become less virially bound. In the latter stages of their evolution, however, cores that remained starless slided down the chart and eventually acquired a configuration that is both, gravitationally bound (Q > 1), as well as virially bound. It appears, under normal circumstances, these cores will remain starless, albeit bounded. The virial chart on the lower right-hand panel for the core L1521F, the VeLLO, where the gas temperature was calculated with puffed-up dust grains is significantly different from that for the starless cores; again in the interest of brevity, we have shown here the virial chart for only this core although calculations were repeated with dust grains of enhanced cross-section for all the cores. Indeed, the extent of gravitational boundedness is much higher for the VeLLO as reflected by a significantly larger magnitude of Q, in comparison with that for the L1521F that rebounded. Furthermore, after its rebound, this latter core acquired a state of approximate equilibrium where the self-gravity was relatively stronger than the external pressure. On the other hand, self-gravity continued to dominate the VeLLO. In the next section, we will discuss the implications of this observation on the ability of a core to form stars.

Figure 11. Plots showing the temporal excursion of test cores on their respective Viral chart. The region rightward of X > 1 corresponds to the virially bound state whereas to the left, virially unbound. Similarly, the region upward of Q > 1, corresponds to the state when a core is dominated by self-gravity, whilst that below, to one that is confined by external pressure.

4 DISCUSSION

4.1. Evolution of prestellar cores

A unified picture of core formation and evolution in which a putative star-forming core is assembled in the outside-in fashion by converging gas flows, has recently been suggested by Gong & Ostriker (2009, 2011). The assembled core in this case remains stable as long as its mass is less than the critical mass, i.e. the critical mass, M _BE, for the marginally stable BES. Once the core acquires mass greater than M _BE, it begins to collapse in the inside-out fashion. There are a few difficulties with this proposition: (i) recent observations of several starless cores have shown that despite individual cores having a density profile resembling that of a unstable BES, they do not show any discernible signs of collapse (e.g. B68, Alves et al. 2001), (ii) it is difficult to perceive how nature could possibly assemble a core that is marginally stable. Furthermore, recent observations of a number of starless cores have also shown that some cores with masses exceeding their thermal Jeans mass or indeed, their critical BE mass do in fact, remain starless, the so-called Super Jeans cores (e.g. Sadavoy et al. 2010a, 2010b), (iii) if nature does indeed prefer to assemble cores via converging supersonic gas flows, in that case, the broadening of spectral emission lines from cores should be significantly broader than observed to reflect the agitated nature of gas within a core, and finally (iv) the inside-out collapse is inconsistent with observations of a number of starless cores that exhibit inwardly directed gas flow (e.g. Tafalla et al. 1998, 2004; Lee et al. 2007, Schnee et al. 2007), with velocity magnitude on the order of ≲ 0.1 km s⁻¹ which is significantly smaller than the typical free-fall velocity.

In the present work, we modelled each of our test cores as a pressure-confined BES ( $\xi _{\text{b}}$ = 3), which was then allowed to evolve by itself in the presence of self-gravity. The gas temperature was calculated in a self-consistent manner by explicitly solving the respective equations of thermal balance for dust and gas. We demonstrated that the eventual fate befalling a core bore largely on the temperature profile it developed rather than on its initial density distribution. In the specific case of the core B68, we showed that merely assembling the test core as a unstable BES was insufficient to induce it to collapse to singularity. We observed, irrespective of the choice of the polytrope used to model this core, with the standard dust grain size, the test core did not collapse even after acquiring a density distribution that mimicked the profile of a unstable BES. In this instance, although our test cores developed a cold central region with gas temperature ≲ 10 K, they remained starless. In fact, we argue that a core can possibly become singular only if gas in its interiors can cool efficiently during its contraction so as to be able to sustain its collapse.

To establish our argument about the importance of gas temperature in the evolution of a prestellar core, we performed two realisations for the test core L1521F. This core is believed to have formed its first hydrostatic core and has been classified as a VeLLO in literature (e.g. Takahashi et al. 2013). In the realisation listed 7 in Table 1, performed using the standard grain size for dust in the interstellar medium, the core was seen to rebound. The core in the realisation listed 8 in Table 1 where a higher dust grain size was adopted, collapsed after the gas close to its centre acquired a temperature between ~ 7 and 10 K. We repeated calculations for all the other starless cores viz., B68, L694-2, L1517B, and L1689-SMM2 with this higher grain size. It was observed that the test cores in all these latter realisations developed cold interiors and collapsed. Evidently, the temperature gradient in a core or in other words, the extent of cooling within a core must somehow determine if whether it will remain starless or become protostellar. This inference begs the question about the possible mechanism/s that might assist gas cooling. It is well-known that gas in a core cools primarily via ro-vibrational transitions of excited molecules and through the interaction between gas and dust, provided of course that the dust is cooler than gas. Detailed modelling of the cooling function for molecular gas in dark clouds by Goldsmith (2001) has shown that the molecular cooling becomes significantly stronger for gas densities upward of ~ 10⁵ cm⁻³, a phase described as thermally super-critical in literature (e.g. Keto & Caselli 2008). This remains true even after accounting for the depletion of molecular species at higher densities.

Next, a cursory glance at Column 2 in Table 2 shows that some of our test cores remained starless despite becoming thermally super-critical. Earlier in Section 3, we have already ruled out the critically stable BES and the uniform density sphere as a model core for the B68 so the corresponding entries need not be considered here. It is however, quite clear from this table that acquiring the thermally super-critical state for a core is not sufficient to ensure its collapse.

Table 2. The observed lifetime of model cores.

^a Free-fall time corresponding to the peak density acquired by the core.

^b Contraction timescale for the core (see text).

^c Listed 7 in Table 1.

4.2. Prestellar core lifetime

Listed in column 3 of Table 2 is the contraction timescale( $t_{\text{c}}$ ), i.e. the timescale over which a core acquires its peak density, for our test cores. We note, this timescale is on the order of ~ 10⁵ yr for each test core. Thus, the model cores in this work appear to have contracted in a quasi-static manner over a few free-fall times and easily survived much longer (the simulations were terminated much later after the cores had contracted). Let us now examine the observational evidence on this count. In one of the early surveys of prestellar cores, Lee et al.(1999) found 220 cores of which 17( ~ 8%) showed signs of in-fall; specifically, seven cores definitely showed evidence of in-fall whilst 10 others showed tentative evidence of in-fall. Later, Gregersen & Evans (2000) surveyed a subset of about 50 cores originally identified by Beichman et al. (1986). This survey identified 6 cores ( ~ 10%) with signs of in-fall. Similarly, Enoch et al. (2007) identified 92 protostellar objects out of 201 cores surveyed collectively in the Ophiuchus, Perseus and the Serpens MCs. In other words, ~ 46% of the cores in this sample had collapsed. More recently, Schnee et al. (2013) reported observations of 26 starless cores in the Perseus MC, 11( ~ 42%) of which showed signs of in-fall.

Now, observationally the total lifetime of starless cores is estimated from the number ratio of cores with and without a protostar. This estimate is on the order of, $t_{\text{SL}}\sim \,$ 0.3–1.6 Myr (Lee & Myers 1999), which is consistent with earlier estimates suggested by Beichman et al. (1986) and Ward-Thompson et al.(1994); see also Ward-Thompson et al. (2007). Using the fraction, $f_{\text{inf}}$ , of cores showing signs of in-fall from literature the typical contraction timescale, $t_{\text{c}}$ , of a prestellar core can be approximated as, $t_{\text{c}}\sim t_{\text{SL}}*f_{\text{inf}}$ , estimates of which vary from 8 × 10⁴ yr for data from Lee et al. (1999), ~ 10⁵ yr for that from Gregersen & Evans (2000), ~ 5 × 10⁵ yr for data from Enoch et al. (2007), and ~ 4 × 10⁵ yr for data from Schnee et al. (2013). Evidently, the estimate of $t_{\text{c}}$ using the data from Lee et al. (1999) is somewhat asynchronous with the magnitude calculated from the rest which are mutually consistent.

Variation in the number of cores showing signs of in-fall is therefore not totally surprising. On the other hand, the contraction timescale, $t_{\text{c}}$ , observed for the test cores in this work is also consistent with the chemical age reported for some starless cores. Maret et al. (2013), for instance, estimated the chemical age for cores L1517B and L1498 on the order of a few 10⁵ yr by radiative transfer modelling of these respective cores, coupled with a detailed chemical network. Nevertheless, the timescale over which the model cores in this work contracted to acquire their respective peak density is consistent with the timescale, $t_{\text{c}}$ , deduced from the surveys reported above. Equally significant is the fact that we could reconcile the observed longevity of starless cores with purely hydrodynamic models without including any additional support from the magnetic field and/or turbulent velocity field.

4.3. Virial state of cores

In a more recent contribution, Pattle et al. (2015) discussed the stability of a sample of cores detected in the Ophiuchus MC with the SCUBA-2 instrument on the JCMT. They found that a majority of the cores in their sample were either gravitationally or virially bound. We introduce here the virial chart, a plot showing the temporal evolution of the virial components of the test core. The virial chart for our test cores show that over the course of its evolution, a core can acquire a configuration that is gravitationally bound, pressure-confined, or virially bound. The rate at which a core evolves of course, depends on its temperature profile. Furthermore, as has already been pointed out by Pattle et al. (2015), the BE mass, M _BE, is not a particularly good indicator of the propensity of a core to become protostellar.

Simulations discussed in this work also demonstrate that the mere fulfilment of the criterion, $\frac{M_{{\rm core}}}{M_{{\rm BE}}}>$ 1, is not sufficient for a core to become singular and in fact, correction of the sound speed in the expression for M _BE to include the velocity dispersion appears to be a better test to determine the stability of a core. The second and third columns in Table 3 call for attention in this regard. Listed in these columns are the virial coefficients for the respective cores at the epoch when they were centrally most concentrated. A cursory glance shows that the respective cores were both gravitationally as well as virially bound with L1517B being the only exception as it remained pressure-confined, as well as virially unbound. Consequently, all these starless cores(except perhaps L1517B), were good enough to become protostellar and yet, remained starless. This leads to the next conclusion: Being gravitationally and/or virially bound, though necessary, is not a sufficient condition for a core to become protostellar. Realisations of the core L1521F, listed 7 and 8 in Table 1, are good examples to illustrate this point.

Table 3. Virial coefficients for test cores(modelled as Bonnor–Ebert sphere with $\xi _{\text{b}}$ = 3) at the epoch when they reach their maximum density.

4.4. The case for enhanced dust grain size

The temperature of gas in a core is sensitive to the strength of the CR ionisation rate. Keto & Caselli (2008) demonstrated that a CR ionisation rate stronger than average raised the gas temperature in the outer regions of a core irrespective of whether it was thermally sub-critical or super-critical. Furthermore, these authors also demonstrated that the temperature of gas close to the centre of a thermally super-critical core was lowered significantly to ~ 6 K when the dust opacity was raised by a factor of four. This, however, did not appear to lower the gas temperature below ~ 10 K close to the centre of a sub-critical core. Presently, there is some tentative evidence to support the idea of enhanced opacity due to relatively larger, i.e. fluffy dust grains. One of the first suggestions in this regard was made by Krügel & Siebenmorgen (1994) after comparing CO and dust observations for a few typical dark clouds. They argued, dust grains in dark clouds likely augmented their size as depleted molecular species gradually adsorbed on the grain surface. These dust grains could possibly coagulate to form the so-called puffed-up, fluffy dust grains with enhanced opacity probably then followed (e.g. Ossenkopf & Henning 1994). Indirect evidence in favour of higher dust opacity has also been suggested on the basis of estimates of core masses. Evans et al. (2001), for instance, derived core masses using the standard dust opacity and found them well in excess of the maximum stable mass for a BES. On the contrary, raising the dust opacity lowered the estimates of core mass such that they were then consistent with the BE stability criterion. Similarly, Keto et al. (2004), who calculated core masses using the N₂H⁺ lines found that opacities of dust grains would have to be significantly higher than that for the ISM in order to be consistent with the core masses calculated in the survey.

The possibility of grain growth can also be studied from the inferred values of the dust opacity spectral index that is known to be sensitive to the maximum grain-size, but insensitive to fine grains smaller than ~ 1 μm. Detailed models that study the growth of dust grains have shown that the dust opacity spectral index lies between ~ 1.5–2 for fine granular dust, then increases to between 2 and 3 for grains typically on the order of ~ 10²– ~ 10³μm and then falls to less than unity for even bigger dust particles. By comparison, the opacity spectral index for dust in the interstellar medium is ~ 2 (see e.g. Testi et al. 2014). Interestingly, in a more recent work, Chen et al. (2015) have reported dust spectral indices between 2 and 3 and between 0 and 1.5 towards a few clumps in the Perseus MC. These magnitudes suggest that dust coagulation leading to the formation of larger grains is perhaps underway in these clumps. In a recent work that lends more support to this conclusion, Juvela et al. (2015) have reported a median dust spectral index of ~ 1.84 for cores detected in the Galactic Cold cores survey carried out with data from the Herschel and Planck.

Another point that favours puffed-up dust grains is the fact that the latter tends to lower gas temperature in a core more uniformly as is visible in the plot on the right-hand panel of Figure 10 for L1521F and indeed for all the other test cores. The plot on the right-hand panel of this figure is more consistent with the temperature( ~ 6–7 K), found close to the centre of typical cores (e.g. Crapsi et al. 2007), and indeed, towards the centre of the VeLLO L1521F (Chitsazzadeh 2014). Puffed-up dust grains have a higher opacity and enhance gas cooling by efficiently coupling with the gas that effectively lowers the gas temperature within a core.

4.5. Other similar work and limitations of this work

The proposition of modelling a prestellar core with a pressure-confined BES is qualitatively similar to a number of models suggested in the past. Like those models, our model also has a central region that has approximately uniform density, beyond which it peters off. Such models have been discussed by Henriksen, André, & Bontemps (1997) and Whitworth & Ward-Thompson (2001). Whilst Henriksen et al. (1997) assumed a constant central density, Whitworth & Ward-Thompson (2001) assumed the central density to be rising according to a power law. Models using ambipolar diffusion have predicted similar density profiles (e.g. Safier, McKee, & Stahler 1997). The model core in this latter case contracted on a timescale longer than popular estimates of contraction time, $t_{\text{c}}$ , and could not reproduce the observed relatively large in-fall velocities. On the other hand, we present a more general model for the initial density distribution of a typical prestellar core. Our model is inspired by the profile of a core extracted from a previously developed numerical simulation aimed at studying the formation of prestellar cores. It does not have any free parameters. Importantly, it is not our endeavour to fit either the density distribution or the distribution of gas temperature within a core at any particular epoch of its evolution. Instead, we follow the temporal evolution of a test core in the presence of self-gravity. By doing so, we derived the density distribution, the temperature profile and the velocity of in-falling gas at different epochs of evolution of a test core. In general, results from our realisations are consistent with the findings reported by various authors in literature for the respective cores.

Limitations. Without a full-fledged chemical network at our disposal, gas temperature in realisations discussed here was calculated by employing various cooling functions to account for heating/cooling of gas and dust. Consequently, the calculations presented here do not account for depletion of various molecular species at higher gas densities. As a result, the cooling rate at higher densities is likely to have been somewhat over-estimated in these calculations. That, however, is unlikely to significantly affect our conclusions, for the cooling function due to depleted molecular species behaves approximately similarly to the function for undepleted species, though the cooling rate is somewhat lower leading to slightly higher(by only a fraction of a Kelvin as a result of CO depletion) gas temperature (Goldsmith 2001; Keto & Caselli 2008). This is simply because the equilibration timescale of the CO is only a fraction of the dynamical timescale of a core; in other words, the depleted CO is rapidly replenished via the inverse process of desorption (Keto & Caselli 2010). Furthermore, employment of a formal radiative-transfer scheme to calculate gas and dust temperature would have permitted us to assess the extent to which the ISRF is attenuated as a function of the column density of the test core. In the present work though, we have assumed constant heating due to the ISRF with a fixed attenuation within a test core and characterised by the parameter, χ ~ 10⁻⁴. This is a reasonable assumption, consistent with magnitudes reported for typical IRDCs (e.g. Goldsmith 2001; Valdivia et al. 2016). We also note that conclusions drawn in this work are for cores that evolve in quasi-static manner so that they are probably only true for cores in relatively quiescent clouds where rapid increase in the magnitude of external pressure is unlikely.

5 CONCLUSIONS

In light of the results obtained from numerical simulations discussed in this work we argue, the fate of a prestellar core bears largely upon its temperature profile. Our realisation of the core B68 with different choices of the initial distribution of gas shows that the likely fate of a test core bears more strongly on the temperature profile that it develops over the course of its evolution. For even the test core B68 modelled as a unstable BES, like other test cores that remained starless in this numerical exercise, rebounded soon after acquiring a centrally peaked density profile. Furthermore, other cores that remained starless were seen to mimic the profile of a unstable BES at the epoch when they were centrally most concentrated and yet, remained starless. We note, some of our test cores viz., L694-2, L1689-SMM2, and L1521F(with standard dust grain size) did in fact become thermally super-critical before rebounding, whilst the remaining two viz., B68 and L1517B remained sub-critical with a peak density less than 10⁵ cm⁻³. It must be noted here that the core L1521F became super-critical and eventually a VeLLO only with a larger dust grain size. In fact, with a larger grain size all other starless cores also developed approximately uniform cold interiors before eventually collapsing. These observations from the present work lead to the conclusion that thermal super-criticality is probably only a necessary condition, though not sufficient, to ensure that a core becomes protostellar.

We suggest, the strength of gas–dust coupling in a core likely holds the key to the eventual fate that befalls the respective core. In this case, we demonstrate, gas temperature within a core is lowered more uniformly. With the choice of grain size for dust typically found in the ISM, we observed that the core L1521F developed a temperature profile similar to that observed in a typical starless core such as the B68. In this latter case, gas temperature in the core was approximately uniform, on the order of 12–14 K and showed a gradual reduction to ≲ 10 K only very close to its centre (see plots in Figure 10). The conclusion about fluffy dust grains is consistent with the earlier suggestion by Keto & Caselli (2008) in this regard. This idea is presently supported by at least tentative evidence from a few detailed studies of dust properties towards some prestellar cores and clumps. Our result here therefore calls for more studies about the nature of dust grains in prestellar cores and in cores that remain starless.

Furthermore, realisations of cores that remained starless in this work also demonstrate that these cores can indeed survive over at least a few free-fall times (see Table 2). In fact, realisations discussed here reproduce the typical core contraction time, the in-fall velocities, as well as the rate of accretion for a typical VeLLO. These simulations also show that the virial state of a core changes over the course of its evolution. We tracked this change on the so-called virial chart that marks the temporal variation of virial components of a core. Cores were seen to evolve from a state that was pressure-confined to one that was gravitationally bound. Similarly, we have also shown that cores make a transition from a state that is less virially bound to one that is more strongly virially bound. Listed in Table 3 are the magnitudes of virial coefficients Q, and X, for respective cores modelled in this work. The magnitudes of these coefficients for cores that remained starless in this work call for attention; not only are these cores virially bound (X > 1), but also relatively strongly bound (Q > 1), except L1517B. This suggests, gravitational and/or virial boundedness of a core is not a sufficient condition for it to collapse, though it is probably necessary. A prospective protostellar core must continue to cool efficiently as it continues to contract in quasi-static manner. In the foreseeable future, we propose to expand the remit of this work by coupling the hydrodynamic calculations discussed here by including the effects of depletion of various molecular species as a core continues to acquire higher densities during its collapse. Also, we should like to examine the impact of ambient environment on the evolution of cores and furthermore, include the magnetic field in these calculations to study how it affects the dynamic evolution of a core.