Issue 
A&A
Volume 601, May 2017



Article Number  A113  
Number of page(s)  9  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201630146  
Published online  15 May 2017 
On the stability of nonisothermal BonnorEbert spheres
III. The role of chemistry in core stabilization
^{1} MaxPlanckInstitute for Extraterrestrial Physics (MPE), Giessenbachstr. 1, 85748 Garching, Germany
email: osipila@mpe.mpg.de
^{2} Department of Physics, PO Box 64, 00014 University of Helsinki, Finland
Received: 28 November 2016
Accepted: 21 March 2017
Aims. We investigate the effect of chemistry on the stability of starless cores against gravitational collapse.
Methods. We combined chemical and radiative transfer simulations in the context of a modified BonnorEbert sphere to model the effect of chemistry on the gas temperature, and study the effect of temperature changes on core stability.
Results. We find that chemistry has in general very little effect on the nondimensional radius ξ_{out}, which parametrizes the core stability. Cores that are initially stable or unstable tend to stay near their initial states, in terms of stability (i.e., ξ_{out} ~ const.), as the chemistry develops. This result is independent of the initial conditions. We can however find solutions where ξ_{out} decreases at late times (t ≳ 10^{6} yr), which correspond to increased stabilization caused by the chemistry. Even though the core stability is unchanged by the chemistry in most of the models considered here, we cannot rule out the possibility that a core can evolve from an unstable to a stable state owing to chemical evolution. The reverse case, where an initially stable core becomes ultimately unstable, seems highly unlikely.
Conclusions. Our results indicate that chemistry should be properly accounted for in studies of starforming regions, and that further investigations of core stability especially with hydrodynamical models are warranted.
Key words: ISM: clouds / ISM: molecules / radiative transfer
© ESO, 2017
1. Introduction
The BonnorEbert 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 nonisothermality 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 radiallyvarying temperature. We refer to a core model constructed using the latter approach as a modified BonnorEbert 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 nonionizing radiation has been studied by Seo et al. (2013) and Seo & Youdin (2016). However, the stability condition of the MBES is a relatively littlestudied 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 nondimensional 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 (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_{2} and He only so that m = 2.33 amu. By making the substitutions (2)(3)one can transform Eq. (1) into nondimensional form: (4)the modified LaneEmden 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 . 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 (5)where represents the thermal component and σ_{s,nt} = 200 m s^{1} represents the nonthermal 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 nondimensional 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 (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 timedependent 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 nonthermal 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^{5} 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 λ = 10^{5} 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 gasdust 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 gasdust 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_{2} 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_{2}. We discuss the effect of changing the initial abundances in Sect. 4.2.
In contrast to Paper II, we consider here a timedependent approach to the chemistry. Starting from the initial MBES (t = t_{0}), we solve the chemistry until some predetermined time step t_{1}, 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_{1} and continue until t_{2}, 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^{7} yr. We set the length of the initial time step to 100 yr, and considered 19 subsequent time steps logarithmically spaced between 10^{2} and 10^{7} 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 intermediatedensity 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. 
Central densities (λ) and boundary pressures (P_{ext}/k) of the core models discussed in this paper.
3. Results
Fig. 3 Lefthand 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. Righthand 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 (lefthand panel), average densities (middle panel), and average temperatures (righthand 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^{5} yr. In two of the cores (B and D) ξ_{out} begins to decrease at t> 10^{5} 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 gasdust coupling, which results in a temperature difference of ~4 K between the two components at a density of ~10^{5} 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 LaneEmden equation (Eq. (4)). We show in the lefthand 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^{6} 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 grainsurface 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 freezeout. 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 freezeout and on the other hand to stronger gasdust 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) (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 which is dominated by . 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 A_{V} 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 singlepoint chemical model assuming n(H_{2}) = 10^{4} 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 intermediatedensity cloud evolve for t = 10^{5} 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^{3} yr to a few × 10^{5} 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 latetime 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). 
Initial abundances of species relevant to the cooling calculations.
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^{2} 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^{6} yr for the molecular envelope surrounding the protostar IRAS 162932422 using the spin states of H_{2}D^{+} as a chemical clock. Other estimates for the ages of molecular clouds range from ~10^{6} 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^{6} yr have been derived also for the cores themselves (WardThompson et al. 1994; Kong et al. 2016). Given that the freefall collapse timescale is ~10^{5} 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^{6} 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 gasgrain chemistry was not found to affect starformation 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 BonnorEbert spheres with a timedependent 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 nondimensional 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 nondimensional 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_{2}). 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 (lowmass) starforming 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
 Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Black, J. H. 1994, in The First Symposium on the Infrared CASP Conf. Ser., 58, 355 [Google Scholar]
 Bonnor, W. B. 1956, MNRAS, 116, 351 [CrossRef] [Google Scholar]
 Brünken, S., Sipilä, O., Chambers, E. T., et al. 2014, Nature, 516, 219 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ebert, R. 1955, Zeitschrift für Astrophysik, 37, 217 [Google Scholar]
 Evans, II, N. J., Rawlings, J. M. C., Shirley, Y. L., & Mundy, L. G. 2001, ApJ, 557, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Galli, D., Walmsley, M., & Gonçalves, J. 2002, A&A, 394, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldsmith, P. F. 2001, ApJ, 557, 736 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hartmann, L., BallesterosParedes, J., & Bergin, E. A. 2001, ApJ, 562, 852 [NASA ADS] [CrossRef] [Google Scholar]
 Hocuk, S., Cazaux, S., Spaans, M., & Caselli, P. 2016, MNRAS, 456, 2586 [NASA ADS] [CrossRef] [Google Scholar]
 Juvela, M., & Ysard, N. 2011, ApJ, 739, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Juvela, M., Padoan, P., & Jimenez, R. 2003, ApJ, 591, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., & Caselli, P. 2008, ApJ, 683, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., & Field, G. 2005, ApJ, 635, 1151 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., Rawlings, J., & Caselli, P. 2014, MNRAS, 440, 2616 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446, 3731 [NASA ADS] [CrossRef] [Google Scholar]
 Kong, S., Tan, J. C., Caselli, P., et al. 2016, ApJ, 821, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Li, P. S., McKee, C. F., & Klein, R. I. 2015, MNRAS, 452, 2500 [NASA ADS] [CrossRef] [Google Scholar]
 Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043 [NASA ADS] [CrossRef] [Google Scholar]
 Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943 [NASA ADS] [Google Scholar]
 Pagani, L., Bacmann, A., Cabrit, S., & Vastel, C. 2007, A&A, 467, 179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pattle, K., WardThompson, D., Kirk, J. M., et al. 2015, MNRAS, 450, 1094 [NASA ADS] [CrossRef] [Google Scholar]
 Pillai, T., Kauffmann, J., Tan, J. C., et al. 2015, ApJ, 799, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Seo, Y. M., & Youdin, A. N. 2016, MNRAS, 461, 1088 [NASA ADS] [CrossRef] [Google Scholar]
 Seo, Y. M., Hong, S. S., & Shirley, Y. L. 2013, ApJ, 769, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Sipilä, O. 2012, A&A, 543, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sipilä, O., Harju, J., & Juvela, M. 2011, A&A, 535, A49 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sipilä, O., Caselli, P., & Harju, J. 2015a, A&A, 578, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sipilä, O., Harju, J., Caselli, P., & Schlemmer, S. 2015b, A&A, 581, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sipilä, O., Harju, J., & Juvela, M. 2015c, A&A, 582, A48 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 WardThompson, D., Scott, P. F., Hills, R. E., & Andre, P. 1994, MNRAS, 268, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Young, K. E., Lee, J.E., Evans, II, N. J., Goldsmith, P. F., & Doty, S. D. 2004, ApJ, 614, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Zucconi, A., Walmsley, C. M., & Galli, D. 2001, A&A, 376, 650 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
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 (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 lefthand 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 righthand 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 LaneEmden 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 latetime 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
Central densities (λ) and boundary pressures (P_{ext}/k) of the core models discussed in this paper.
All Figures
Fig. 1 Boundary pressure of a BES with central density λ = 10^{5} 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 Lefthand 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. Righthand 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 (lefthand panel), average densities (middle panel), and average temperatures (righthand 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 A_{V} 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 lefthand 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 righthand 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.