© ESO, 2017

1. Introduction

The Bonnor-Ebert sphere (hereafter BES; Bonnor 1956; Ebert 1955), essentially an isothermal gas sphere in hydrostatic equilbrium, has been used in many studies to represent the physical structure of starless/prestellar cores, perhaps most famously in the case of B68 by Alves et al. (2001). It is however an observationally established fact that real cores are not isothermal (Crapsi et al. 2007; Pagani et al. 2007). BES models including the possibility of non-isothermality have been studied by several authors (Evans et al. 2001; Zucconi et al. 2001; Galli et al. 2002; Keto & Field 2005; Juvela & Ysard 2011). Keto & Caselli (2008, 2010) and Keto et al. (2014, 2015) have studied the particular case of L1544 using hydrodynamics coupled with (simplified) chemistry. In the previous studies considering static cloud models, the dust or gas temperature was either calculated based on the density profile of a BES so that the model retained the physical structure of the isothermal core, or the equations determining the density profile were modified to include the possibility of a radially-varying temperature. We refer to a core model constructed using the latter approach as a modified Bonnor-Ebert sphere (MBES).

When attempting to reproduce observations using a particular core model, it is necessary to consider whether or not the model core is stable. The stability condition of the BES is well known (Bonnor 1956); cores with nondimensional radius ξ_out ≳ 6.45 represent unstable configurations. The stability of a BES with internal motions or subjected to non-ionizing radiation has been studied by Seo et al. (2013) and Seo & Youdin (2016). However, the stability condition of the MBES is a relatively little-studied subject. In Sipilä et al. (2011, hereafter Paper I), we presented an analytical expression for the stability condition of the MBES based on an extension of Bonnor’s equations to accommodate radial changes in the temperature. In that paper, we made the assumption that T_dust = T_gas, which is approximately valid only at high density. In Sipilä et al. (2015c, Paper II), we expanded the validity range of the stability analysis by switching from T_dust to T_gas in the relevant equations, made possible by including chemistry in the calculation framework. In both papers it was found that the stability condition of the MBES, parametrized by the non-dimensional radius ξ_out, is close to that of the BES.

In this paper, we investigate the potential impact that chemistry may have on the stability of an evolving MBES, that is, instead of studying the global effects on ξ_out in cores with varying masses as was done in our previous papers, we concentrate on the evolution of ξ_out for single cores. We select cores with varying initial nondimensional radii (stable and unstable configurations based on the results of Papers I and II), and study the evolution of the physical properties of the core in tandem with the chemical evolution of the gas. The particular aim of the present paper is to find out if a core can evolve from the unstable regime into the stable regime, or vice versa, owing to chemical evolution. This important question has received little attention in the literature so far.

The paper is organized as follows. In Sect. 2 we introduce our model and explain how the stability calculations are carried out in practice. The results of our analysis are presented in Sect. 3. In Sect. 4 we discuss our results and in Sect. 5 we give our conclusions.

2. Model

2.1. Basic formulae

The properties of the MBES are discussed in detail in Papers I and II. Here, we recall the relevant formulae. The density distribution of the MBES is given by $\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left(\frac{{\mathit{r}}^{\mathrm{2}}}{\mathit{\rho }}\left[\mathit{T}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{r}}\mathrm{+}\mathit{\rho }\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{r}}\right]\right)\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi Gm\rho }}{\mathit{k}}\mathit{,}$(1)which assumes the ideal gas equation of state and hydrostatic equilibrium. In the above, m and k denote the average molecular mass of the gas and the Boltzmann constant, respectively, and ρ denotes the medium density. Even though in this paper we employ a detailed chemical model that includes a wide range of chemical species, we make the approximation in the calculation of the molecular mass that the gas consists of H₂ and He only so that m = 2.33 amu. By making the substitutions $\mathit{\rho }\mathrm{=}\frac{\mathit{\lambda }}{\mathit{\tau }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}}$(2)$\mathit{r}\mathrm{=}{\mathit{\beta }}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\lambda }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{\xi ,}$(3)one can transform Eq. (1) into nondimensional form: ${\mathit{\xi }}^{-2}\frac{\mathrm{d}}{\mathrm{d}\mathit{\xi }}\left[{\mathit{\xi }}^{\mathrm{2}}\mathit{\tau }\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}\mathit{\xi }}\right]\mathrm{=}\frac{\mathrm{1}}{\mathit{\tau }}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }}\mathit{,}$(4)the modified Lane-Emden equation. In the above, ξ represents the nondimensional radius of the core (ξ and ψ(ξ) are both dimensionless); λ represents the central density; τ(ξ) = T(ξ) /T_c represents the radial (gas) temperature profile.

The variable β in Eq. (3) represents the (scaled) velocity dispersion at the center of the core. In Papers I and II we assumed that the dispersion is due to purely thermal motions, and hence that $\mathit{\beta }\mathrm{=}{\mathit{\sigma }}_{\mathrm{s}\mathit{,}\mathrm{th}}^{\mathrm{2}}\mathit{/}\mathrm{4}\mathit{\pi G}\mathrm{=}\mathit{k}{\mathit{T}}_{\mathrm{c}}\mathit{/}\mathrm{4}\mathit{\pi Gm}$. Here, we introduce a turbulent component to the velocity dispersion with a constant value of 200 m s^-1, which is in the range of typical values observed in dense cores (e.g., Pattle et al. 2015). Therefore in this paper $\mathit{\beta }\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{s}\mathit{,}\mathrm{th}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{s}\mathit{,}\mathrm{nt}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi G}}\hspace{0.17em}\mathit{,}$(5)where ${\mathit{\sigma }}_{\mathrm{s}\mathit{,}\mathrm{th}}^{\mathrm{2}}\mathrm{=}\mathit{k}{\mathit{T}}_{\mathrm{c}}\mathit{/}\mathit{m}$ represents the thermal component and σ_s,nt = 200 m s^-1 represents the non-thermal component. We discuss the effect of assuming a purely thermal velocity distribution on our results in Sect. 4.3.

Equation (4) can be solved numerically when the temperature profile is known. Here, we determine the dust and gas temperature profiles by using radiative transfer calculations (see below). The following boundary conditions are imposed at the core center: ψ = 0, dψ/ dξ = 0, τ = 1, and dτ/ dξ = 0, which ensure that T = T_c and ρ = λ at ξ = 0, that is, at the center of the core.

2.2. Method

In Paper II, we examined the stability of the MBES under the assumption that the non-dimensional radius ξ_out and the total mass M of a model core stay constant as the gas temperature changes with the chemical evolution of the gas. As discussed in Paper II, these assumptions may not hold in reality, which means that the results presented in that paper represent the possible configurations that MBESs of various masses can exist in when embedded in regions of different chemical ages. Here, we change our approach and keep the external pressure p_ext constant so that the model represents a core evolving inside a fixed environment. The boundary pressure can be written as ${\mathit{p}}_{\mathrm{ext}}\mathrm{=}{\mathit{\rho }}_{\mathrm{out}}\mathrm{(}{\mathit{\sigma }}_{\mathrm{s}\mathit{,}\mathrm{th}}^{\mathrm{2}}\mathrm{+}{{\mathit{\sigma }}_{\mathrm{s}\mathit{,}\mathrm{nt}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi G\beta \lambda }{\mathrm{e}}^{\mathrm{-}{\mathit{\psi }}_{\mathrm{out}}}}{{\mathit{\tau }}_{\mathrm{out}}}\mathit{,}$(6)where τ(ξ) = T(ξ) /T_c (see Paper I) and the subscript “out” refers to a value taken at the nondimensional outer radius of the core, ξ_out, which is here a time-dependent variable contrary to the scenario discussed in Papers I and II. Note that the expression of p_ext in the above is different than the one given in Paper I (Eq. (4) therein) because here the β parameter also includes a non-thermal contribution.

In addition to the external pressure, we are free to fix one of two properties of the core: the central density or the mass (while ξ is not fixed). Here we choose to fix the central density. This choice allows us to determine the nondimensional radius ξ_out using Eq. (6). By tracking the changes in ξ_out as the gas temperature changes with the chemical evolution of the core, we can study the potential impact of the chemistry on the core stability. We explain below how this is achieved in practice.

Setting the core mass and the external pressure as the problem constants leads sometimes to situations where a solution cannot be found, that is, the evolution of ξ_out cannot be determined. We discuss models with fixed mass in more detail in Appendix A.

2.3. Determining the core evolution

We constructed the initial MBES by first constructing a BES with predetermined values of central density, (gas and dust) temperature, and external pressure. The choice of the external pressure is, however, not completely free. This is demonstrated in Fig. 1, where we plot the boundary pressure of a BES as a function of ξ_out. Here the central density is set to λ = 10⁵ cm^-3 and the temperature (constant throughout the core) to T = 10 K; we also assumed a turbulent component corresponding to a velocity dispersion of σ_s,nt = 200 m s^-1. Different values of ξ_out correspond to different configurations that a core with the chosen values of central density and temperature can exist in. The steepness of the density gradient within a core increases with ξ_out. Also, the core mass, for example, is different for each value of ξ_out. The value of the boundary pressure for the configuration with ξ_out = 0 represents the maximum value of p_ext possible for a BES with the chosen parameters. The lower the boundary pressure, the more unstable the core configuration is, since solutions with ξ_out ≳ 6.5 represent unstable cores (Bonnor 1956).

Fig. 1 Boundary pressure of a BES with central density λ = 105 cm-3, T = 10 K as a function of nondimensional radius ξout.

Here, when constructing the initial MBES, we chose the boundary pressure of the BES to be lower than the maximum value allowed by the choices of central density and initial central temperature. We then calculated the dust temperature of the BES using the Ossenkopf & Henning (1994) dust model (see also Papers I and II). The initial MBES was constructed in the following steps: 1) solve Eq. (4); 2) find the value of ξ that satisfies Eq. (6); 3) calculate the density profile using Eq. (2). After the initial MBES was set up, we proceeded to calculate the chemical evolution in the core following the methods laid out in Paper II; we separated the model core into spherical shells and calculated the chemical evolution separately in each shell. The initial conditions of the model are discussed in Sect. 4.2. The abundances of various cooling species were extracted from the chemical modeling results and used as input for the radiative transfer code LOC (Juvela, in prep.) to determine the gas temperature profile, which was in turn used to update the MBES following the steps described above. The gas temperature is determined by the relative strengths of heating and cooling processes. The gas is assumed to be heated by cosmic rays, and cools by line radiation and loss of energy in collisions with dust grains. In this paper we used the gas-dust coupling scheme of Goldsmith (2001). We note that a stronger coupling between gas and dust (e.g., Young et al. 2004) would reduce the role of chemistry in determining the gas temperature because it increases the range of physical conditions where the gas-dust coupling is important.

We assumed in our models that the visual extinction at the edge of the core is A_V = 2 mag. However, the radiative transfer calculations also include a layer outside the core itself that absorbs line radiation, and can emit it back into the core, preventing the gas temperature from dropping to very low values (6−7 K) at the edge of the core at early times when the abundances of the main coolants are high. The external molecular layer is assumed to have the same density as the edge of the core, and its thickness is set so to correspond to one magnitude of visual extinction. The layer is not extended all the way to A_V = 0 mag because there, H₂ should dissociate, and in this work we only desire to model what happens inside the core. However, the inclusion or exclusion of the external molecular layer, or the assumed amount of extinction at the core edge, have only a very limited effect on our results (see Sects. 3 and 4.1). We note that in this paper we do not model possible temporal changes in dust properties. The unattenuated ISRF spectrum is taken from Black (1994).

We adopt the chemical model from Sipilä et al. (2015a,b), although in this paper we do not include deuterium in the chemical calculations in order to speed up the process. We adopt the initial chemical abundances from Table 3 in Sipilä et al. (2015a), where the gas is initially atomic with the exception of hydrogen, which is in H₂. We discuss the effect of changing the initial abundances in Sect. 4.2.

In contrast to Paper II, we consider here a time-dependent approach to the chemistry. Starting from the initial MBES (t = t₀), we solve the chemistry until some predetermined time step t₁, and after solving the gas temperature and updating the density profile of the core following the steps 1 to 3 discussed above, we pick up the chemistry at t₁ and continue until t₂, and so on. In this way we can track how the nondimensional radius ξ_out evolves as a function of time, and study whether the stability of the core is affected, using also the results from Papers I and II. We studied four combinations of central densities and boundary pressures (see Table 1), allowing chemical development up to 10⁷ yr. We set the length of the initial time step to 100 yr, and considered 19 subsequent time steps logarithmically spaced between 10² and 10⁷ yr. Tests varying the amount of time steps were carried out and it was found that our main results, described below, remained virtually unaffected by changes in the time resolution.

The values of the core parameters (λ, p_ext) were chosen in such a way that we obtained a range of core masses and nondimensional radii (see below). We concentrated on intermediate-density cores to bring out the effect of the chemistry; at high density the dust and gas are efficiently coupled through collisions and the role of the chemistry in determining the gas temperature diminishes.

Fig. 2 Nondimensional radii ξout of the four cores studied in this paper as functions of time. The core models (see Table 1) are indicated in the figure. The black lines represent our fiducial model described in the text, while the blue lines represent a model where the external molecular layer outside the core is not included.

3. Results

Fig. 3 Left-hand panel: function ψ for core A at selected time steps, indicated in the middle panel of the figure. The green horizontal lines mark the value of ψ that satisfies Eq. (6) at each time step, and the red circles are a guide for the eye to locate the corresponding value of ξout (see text). Middle panel: dust (red) and gas (blue) temperature profiles for core A as functions of distance from the core center at different times. Right-hand panel: density profile for core A as a function of distance from the core center at different times. The quoted times are approximate because the time steps considered in the model, save for the first and last, are not exact multiples of ten.

Fig. 4 Outer radii (left-hand panel), average densities (middle panel), and average temperatures (right-hand panel) of the four cores studied in this paper (indicated in the figure) as functions of time.

We plot in Fig. 2 the evolution of ξ_out as a function of time for the four core models tabulated in Table 1. A clear trend is seen in all of the cores; ξ_out tends to stay approximately constant until t ~ 10⁵ yr. In two of the cores (B and D) ξ_out begins to decrease at t> 10⁵ yr, while cores A and C present somewhat more complex behavior. The evolution of ξ_out can be understood by studying the evolution of the function ψ (which is related to the density profile; see Eq. (2)). Any changes in ψ are tied to the evolution of the temperature profile. As an example, we show in Fig. 3 the function ψ, the temperature profile, and the density profile for core A for selected time steps in the model. At early times, the high abundance of CO helps to maintain a gas temperature of less than 10 K at the center of the core, but as CO depletes, the gas temperature is determined solely by the gas-dust coupling, which results in a temperature difference of ~4 K between the two components at a density of ~10⁵ cm^-3 (see also Goldsmith 2001). At late times, the central temperature of the core changes very little as a function of time (i.e., β~ const.), and so the value of ξ_out is controlled mainly by ψ_out and τ_out (see Eq. (6)), with the former being determined by the temperature profile through the Lane-Emden equation (Eq. (4)). We show in the left-hand panel of Fig. 3 the solution of Eq. (6), i.e., ψ_out = −log (p_extτ_out/ 4πGβλ) for different time steps; the value of ξ_out at each time step is found at the intersection of a horizontal line with a ψ curve of the same linestyle.

Note that the gas temperature rises sharply toward the core edge at t ≳ 10⁶ yr. This effect is a result of the depletion of the cooling species onto grain surfaces, which takes a long time at the low density of the core edge. Then, most of the carbon is locked in grain-surface methane, which is not readily broken into C or CO by photodissociation, depriving the gas of the main coolant species.

To complement the results discussed above, we show in Fig. 4 the outer radii, average densities, and average temperatures (as calculated with Eq. (B.4) in Paper I) of the four cores studied here. The evolution of the outer radius follows that of ξ_out. Cores B and D, which show an eventual decrease of ξ_out, display contraction and an associated increase in average density at late times. The evolution of the average temperature is similar for all cores, showing an increase of 2−3 K from early times to late times owing to molecular freeze-out. However, cores B and D are warmer than cores A and C because they start at lower values of ξ_out and hence their density contrast (central density vs. edge density) is lower. The higher densities toward the edge translate on the one hand to stronger molecular freeze-out and on the other hand to stronger gas-dust coupling. Both effects lead to increased average gas temperatures with respect to the cores with the same central density but higher density contrast.

Figure 2 also shows the results of calculations where the external molecular layer is not included in the modeling. This change decreases the gas temperature at the core edge because of the increased photon escape probability. However, the two approaches to the details of the radiative transfer modeling produce the same general conclusions on the evolution of the core stability as a function of time. In particular, this shows that the temporal variation in ξ_out at late times is a real effect of the model and is not induced by our underlying assumptions.

Fig. 5 Masses of the four cores studied in this paper (indicated in the figure) as functions of time.

Because we fixed the boundary pressure and central density, the core mass is not a constant. The mass is given by (see Papers I and II) $\mathit{M}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{\beta }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\lambda }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{\xi }}_{\mathrm{out}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\tau }}_{\mathrm{out}}\hspace{0.17em}(\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}\mathit{\xi }}{)}_{\mathrm{out}}\mathrm{·}$(7)We plot in Fig. 5 the masses of the four cores considered in this paper as functions of time. Evidently, the evolution of the core mass is similar to that of ξ_out. This is expected based on Eq. (7): the central density is a constant of the model and the central temperature varies generally very little as a function of time (at least at late times, see Fig. 3), so that the change in mass is governed mainly by ${\mathit{\xi }}_{\mathrm{out}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\tau }}_{\mathrm{out}}\hspace{0.17em}{\left(\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}\mathit{\xi }}\right)}_{\mathrm{out}}$ which is dominated by ${\mathit{\xi }}_{\mathrm{out}}^{\mathrm{2}}$. The changes in core mass imply that the core exchanges mass with its environment in order to meet the condition of fixed central density while the outer radius evolves.

The results presented here show that the nondimensional radius ξ_out tends to stay approximately constant with time – or to decrease at late times – in spite of chemical evolution within the core. Thus, chemical evolution may make no difference to the stability of the core, or may in some cases make it less likely for a core to collapse gravitationally. We discuss this issue further in Sect. 4.4.

4. Discussion

4.1. Effect of the external medium

The present analysis is based on the assumption that the central density and the boundary pressure remain constant during core evolution. From a physical standpoint, considering the central density as a constant is reasonable if we aim to model a cloud core in quiescent conditions where little dynamical evolution should take place. However, setting the boundary pressure to a constant value means that we also constrain the properties of the medium outside the core. As mentioned in Sect. 2.3, we assumed in the above analysis that the (gas) temperature of the external medium is equal to whatever value applies at the core boundary at any given time step. This implies also that the medium density outside the core is equal to that at the core edge, assuming that the turbulent velocity component has the same constant value in the core and in the surrounding cloud. Here we aim to study the temporal evolution of ξ_out in a general setting and we do not attempt to fit our results to any observed core. We simply assume that the properties of the external medium are at all times such that our model assumptions are fulfilled.

However, for the sake of consistency with Papers I and II, we studied the effect of varying the assumed visual extinction at the core edge on our results. The extinction (assumed equal to 2 mag in the results presented above) is expected to influence our results, because it determines the strength of the external heating of the core and thus contributes to the gas temperature profile. The external heating in the present model is provided by the ISRF (Black 1994), and the photoelectric effect (see Juvela et al. 2003; Juvela & Ysard 2011). Increasing the extinction will decrease the gas temperature especially near the core edge, affecting the solution to Eq. (4) and thus the derived value of ξ_out. To quantify the possible effect on our results, we ran the core A model with three different external A_V values (2, 5 and 10 mag), corresponding to situations where the model core lies progressively deeper inside a parent cloud. In these calculations the external molecular layer was not included. The results of this test are shown in Fig. 6. Evidently, increasing the extinction increases ξ_out, but the difference between the test cases is very small, in fact clearly smaller than the difference arising from including or excluding the external molecular layer (Fig. 2).

Fig. 6 Nondimensional radius ξout of core A as a function of time, assuming different values of AV at the core edge.

4.2. Initial conditions

The analysis presented above is based on the assumption that the density and temperature structures of the cores correspond to modified BESs, that is, that a (dust) temperature gradient exists already in the beginning of the core evolution. Furthermore, we assumed that the gas is initially atomic with the exception of hydrogen, which is entirely in molecular form. In this section we explore the effect of changing these initial assumptions on our results.

Although CO is usually the main gas coolant (Goldsmith 2001), cooling by atomic C can also be important at early times (Sipilä 2012). Therefore, it is not immediately clear if our results would be greatly affected if we assumed that (some) CO exists already in the initial state of the core, which is equivalent to the reasonable assumption that CO production has taken place preceding the formation of the core itself. To test this hypothesis we repeated the stability calculations presented above, now assuming initial abundances which were determined by a single-point chemical model assuming n(H₂) = 10⁴ cm^-3, T = 12 K, and A_V = 2 mag. This model roughly represents the conditions in a molecular cloud and is here referred to as the “molecular model”. We started from an initially atomic chemical state and let the intermediate-density cloud evolve for t = 10⁵ yr, and used the resulting abundances as initial abundances in the stability calculation.

The initial abundances for the species relevant to the cooling calculations are given in Table 2. Because of the low visual extinction (A_V = 2 mag), CO is efficiently destroyed by photons and its abundance stays at a constant level from ~10³ yr to a few × 10⁵ yr, never exceeding that of atomic carbon, which also stays more or less constant in abundance during this time interval. Therefore, the time at which the abundances are extracted is of little consequence here. For higher values of external A_V, translating to less efficient photoprocesses, CO would be the main carbon reservoir at this density.

Figure 7 shows the nondimensional radii ξ_out as functions of time in the molecular model. The results are very close to those presented in Fig. 2; the ξ_out values of each core hardly change when molecular initial abundances are considered. The cooling power of C (in our fiducial model) and CO (in the initially molecular model) is similar at early times so that the model results do not greatly deviate from one another; all of the other relevant quantities (ψ, mass etc.) are also only marginally affected by the choice of the initial abundances. We emphasize that changes in the initial abundances do not cause a late-time increase of ξ_out, which either stays approximately constant or decreases eventually.

Fig. 7 Nondimensional radii ξout of the four cores studied in this paper as functions of time. Linestyles are the same as those in Fig. 2. Red lines reproduce the data from Fig. 2, while blue lines correspond to calculations in the molecular model (see text).

We also investigated the possibility that the initial core is isothermal, that is, corresponds to a BES. We found that this assumption also makes very little difference to our results; as long as we choose the same boundary pressure and central density in both cases, starting from a BES or an MBES yields the same ξ_out values within ~1% after t ~ 10² yr, that is, the first time step in the model. Already at the first time step the difference between the models is less than ~10%. Our main results thus seem robust given their weak dependence on the initial conditions.

4.3. Models without turbulent pressure

For completeness, we ran another set of calculations eliminating the nonthermal component from Eqs. (5) and (6) to study whether or not this has a fundamental effect on our results. Figure 8 shows the results of these calculations. It is observed that the ξ_out values are generally scaled down with respect to the model that includes turbulent support. This result is expected because removing the turbulent component reduces the maximum boundary pressure allowed for a given combination of central density and temperature, and thus the boundary pressures used (Table 1) correspond to lower values of ξ_out in a model where turbulent pressure is neglected (see also the discussion in Sect. 2.3). Our general conclusions are however the same in both types of model; ξ_out tends to stay constant or decrease at late times, and does not increase ultimately.

Fig. 8 Nondimensional radii ξout of the four cores studied in this paper as functions of time, calculated with a model including thermal velocity dispersion only. Linestyles are the same as those in Fig. 2. Red lines reproduce the data from Fig. 2, while blue lines correspond to calculations in the molecular model (see text).

4.4. Possible core evolution towards the stable regime

The lifetimes of molecular clouds and cloud cores are not strictly constrained. Recently, Brünken et al. (2014) derived an age of at least 10⁶ yr for the molecular envelope surrounding the protostar IRAS 16293-2422 using the spin states of H₂D⁺ as a chemical clock. Other estimates for the ages of molecular clouds range from ~10⁶ yr (Hartmann et al. 2001) to up to ~10 Myr (Mouschovias et al. 2006). The age of a molecular cloud sets a natural upper limit to the ages of cloud cores residing within it. Lifetimes of up to ~10⁶ yr have been derived also for the cores themselves (Ward-Thompson et al. 1994; Kong et al. 2016). Given that the free-fall collapse timescale is ~10⁵ yr for typical starless core densities, it seems clear that stabilizing forces must be at play for the cores/clouds to live up to Myr timescales. Magnetic fields may also play a significant role in the dynamics of molecular clouds by providing pressure to act as a means of support (for observational and theoretical discussions see, e.g., Mouschovias et al. 2006; Pillai et al. 2015; Li et al. 2015).

Our results indicate that significant stabilization as a result of chemical evolution is not expected. Of the four cores studied here, cores B and D show decreasing values of ξ_out at late times, but these cores are stable to begin with. Cores A and C, which are initially unstable, show oscillating behavior at late times, but do not deviate significantly from their initial states in terms of stability. Based on our results one cannot rule out the possibility of a core configuration that starts out unstable and evolves into a stable state. However, our results do imply that such a transition would likely occur over a timescale of ~10⁶ yr, which agrees with the upper end of the range of expected cloud lifetimes estimated with theoretical and observational arguments. The fact that in general the chemistry has little to no effect on the stability of the cores is in line with recent calculations by Hocuk et al. (2016), where the consideration of gas-grain chemistry was not found to affect star-formation rates, although their study was limited to the cloud scale and did not probe core evolution.

Finally we note that our models represent equilibrium solutions, and it is implicitly assumed that the dynamical timescales are shorter than the chemical timescales so that the cores are able to physically adapt to the changes in the chemistry. This type of an idealized situation may not occur in real cores. To investigate the role of chemistry on core stability in a more realistic setting, the chemical and radiative transfer models should be combined with a (magneto)hydrodynamical model of core evolution. A similar study has been undertaken previously by, for example, Keto & Field (2005). However, this type of investigation is beyond the scope of the present paper, and is reserved for future work.

5. Conclusions

We studied the stability of modified Bonnor-Ebert spheres with a time-dependent model that combines chemical and radiative transfer calculations for a consistent determination of the gas temperature as a function of time. The stability of the model cores was determined based on the results of our previous papers (Sipilä et al. 2011, 2015c). Unlike in our previous works, which were based on the conservation of core mass and non-dimensional radius ξ_out, we assumed here that the central density and external pressure are constant, which describes the evolution of a core in a fixed quiescent environment. We assumed that the internal velocity dispersion is composed of a thermal and a nonthermal component.

We found that the non-dimensional radius ξ_out tends to remain approximately constant as a function of time, or in some cases to decrease late into the evolution of the core. An increase in instability, at least in a reasonable timescale, seems highly unlikely based on our models. However, the most likely scenario is that the chemistry does not make a significant difference to the core stability. Our fiducial model starts with an initially atomic gas (with the exception of hydrogen which is in H₂). Changing the initial conditions was found to influence the values of ξ_out that we derive, but does not change our conclusions of the evolution of core stability.

We cannot completely rule out the possibility that chemistry can in certain conditions play a role in determining the lifetime of a starless core, particularly in cases where the core is on the border of being gravitationally stable or unstable. Further studies of the interplay between chemistry and core stability are called for, and may shed more light on the (low-mass) star-forming process in general. Here we considered equilibrium solutions with the assumption that the dynamical timescales are shorter than the chemical timescales. More comprehensive studies of the connection between chemistry and core stability should incorporate a description of dynamics in the form of (magneto)hydrodynamic modeling.

Acknowledgments

We thank the anonymous referee for useful comments and suggestions that helped to improve the manuscript. P.C. acknowledges financial support of the European Research Council (ERC; project PALs 320620). M.J. acknowledges the support of the Academy of Finland Grant No. 285769.

References

Appendix A: Models with fixed mass and external pressure

As an alternative to the method discussed in the main part of this paper, where we fix the central density and boundary pressure of the core, we could instead fix the core mass and external pressure. In this case we eliminate λ from Eqs. (6) and (7), which yields the relation ${\mathit{\xi }}^{\mathrm{2}}\hspace{0.17em}\left(\frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}\mathit{\xi }}\right)\hspace{0.17em}{\mathit{\tau }}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{\psi }\mathit{/}\mathrm{2}}\mathrm{=}\frac{\mathit{M}}{\mathrm{(}\mathrm{4}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\hspace{0.17em}{\mathit{\beta }}^{-2}\hspace{0.17em}{\mathit{G}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{p}}_{\mathrm{ext}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(A.1)The new value of ξ_out is determined by the function on the left, while the right side of the equation is simply a number. The most meaningful comparison between the two approaches (fixed central density vs. fixed mass) is achieved by starting from the same initial condition and finding out if the evolution of ξ_out proceeds differently in the two cases. We tested this alternative method for the four cores studied in the paper (Table 1).

Fig. A.1 The left-hand side (LHS) of Eq. (A.1) as a function of ξ for the first time step for the core B model assuming fixed core mass. The dashed line marks the value of the right-hand side of Eq. (A.1).

The initial state of core B corresponds to M ~ 1.3 M_⊙ and ξ_out ~ 3.4. We proceeded with the same workflow as for the fixed central density case, that is, we calculated the chemistry and the gas temperature followed by a solution of the Lane-Emden equation. This gives the shape for the left side of Eq. (A.1) shown in Fig. A.1. One solution exists at ξ_out ~ 2.5, which is now the new nondimensional radius of the core. Continuing the iterative process leads to the evolution of ξ_out shown in Fig. A.2, where we also show the evolution of ξ_out for core D. The late-time fluctuation in ξ_out for both cores corresponds to an oscillation of the central density by a factor of ~3. The timescale of the oscillation is affected by the chosen time resolution.

The fixed mass method does not work in all cases. For core A the situation at the first time step is such that no solution to Eq. (A.1) can be found. For core C the same situation occurs at the second time step. The fixed mass method is more sensitive to the shape of the temperature profile, and to how the temperature profile continues beyond the core boundary, than the fixed central density method. Therefore the former is more ambiguous. It is conceivable that solutions could be found also for cores A and C if different model assumptions were adopted.

Fig. A.2 Nondimensional radii ξout of cores B and D (indicated in the figure) as functions of time assuming fixed core mass.

All Tables

All Figures

	Fig. 1 Boundary pressure of a BES with central density λ = 105 cm-3, T = 10 K as a function of nondimensional radius ξout.
In the text

	Fig. 2 Nondimensional radii ξout of the four cores studied in this paper as functions of time. The core models (see Table 1) are indicated in the figure. The black lines represent our fiducial model described in the text, while the blue lines represent a model where the external molecular layer outside the core is not included.
In the text

Fig. 3 Left-hand panel: function ψ for core A at selected time steps, indicated in the middle panel of the figure. The green horizontal lines mark the value of ψ that satisfies Eq. (6) at each time step, and the red circles are a guide for the eye to locate the corresponding value of ξout (see text). Middle panel: dust (red) and gas (blue) temperature profiles for core A as functions of distance from the core center at different times. Right-hand panel: density profile for core A as a function of distance from the core center at different times. The quoted times are approximate because the time steps considered in the model, save for the first and last, are not exact multiples of ten.

In the text

	Fig. 4 Outer radii (left-hand panel), average densities (middle panel), and average temperatures (right-hand panel) of the four cores studied in this paper (indicated in the figure) as functions of time.
In the text

	Fig. 5 Masses of the four cores studied in this paper (indicated in the figure) as functions of time.
In the text

	Fig. 6 Nondimensional radius ξout of core A as a function of time, assuming different values of AV at the core edge.
In the text

	Fig. 7 Nondimensional radii ξout of the four cores studied in this paper as functions of time. Linestyles are the same as those in Fig. 2. Red lines reproduce the data from Fig. 2, while blue lines correspond to calculations in the molecular model (see text).
In the text

	Fig. 8 Nondimensional radii ξout of the four cores studied in this paper as functions of time, calculated with a model including thermal velocity dispersion only. Linestyles are the same as those in Fig. 2. Red lines reproduce the data from Fig. 2, while blue lines correspond to calculations in the molecular model (see text).
In the text

	Fig. A.1 The left-hand side (LHS) of Eq. (A.1) as a function of ξ for the first time step for the core B model assuming fixed core mass. The dashed line marks the value of the right-hand side of Eq. (A.1).
In the text

	Fig. A.2 Nondimensional radii ξout of cores B and D (indicated in the figure) as functions of time assuming fixed core mass.
In the text