On the stability of nonisothermal BonnorEbert spheres
II. The effect of gas temperature on the stability
^{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: 27 May 2014
Accepted: 23 August 2015
Aims. We investigate the stability of nonisothermal BonnorEbert spheres in the context of a model that includes a selfconsistent calculation of the gas temperature. In this way, we can discard the assumption of equality between the dust and gas temperatures and study the stability as the gas temperature changes with the chemical evolution of the cooling species.
Methods. We use a gasgrain chemical model to calculate the chemical evolution. The model includes a timedependent treatment of depletion onto grain surfaces, which strongly influences the gas temperature as the main coolant molecule CO depletes from the gas. The dust and gas temperatures are solved with radiative transfer calculations. For consistent comparison with previous work, we assume that the cores are deeply embedded in a larger external structure, corresponding to visual extinction A_{V}^{ext} = 10 mag at the core edge. We also study the effect of lower values of A_{V}^{ext}.
Results. We find that the critical nondimensional radius ξ_{1} derived here, which determines the maximum density contrast between the core center and the outer boundary, is similar to our previous work where we assumed T_{dust} = T_{gas}. Here, the ξ_{1} values lie below the isothermal critical value ξ_{0} ~ 6.45, but the difference is less than 10 %. We find that chemical evolution does not notably affect the stability condition of lowmass cores (< 0.75 M_{⊙}), which have high average densities and a strong gasgrain thermal coupling. In contrast, for higher masses the decrease in cooling due to CO depletion causes substantial temporal changes in the temperature and in the density profiles of the cores. In the mass range 1−2 M_{⊙}, ξ_{1} decreases with chemical evolution, whereas above 3 M_{⊙}, ξ_{1} instead increases with chemical evolution. We also find that decreasing A_{V}^{ext} strongly increases the gas temperature, especially when the gas is chemically old, and this causes ξ _{1} to increase with respect to models with higher A_{V}^{ext}. However, the derived ξ_{1} values are still close to ξ_{0}. The density contrast between the core center and edge derived here varies between 8 and 16 depending on core mass and the chemical age of the gas, compared to the constant value ~14.1 for the isothermal BES.
Key words: radiative transfer / ISM: clouds / astrochemistry / ISM: molecules
© ESO, 2015
1. Introduction
The BonnorEbert sphere (Bonnor 1956; Ebert 1955; hereafter BES), which is an isothermal gas sphere in hydrostatic equilibrium, has been used successfully to approximate the density structures of prestellar cores (Bacmann et al. 2000; Alves et al. 2001; Kandori et al. 2005; Marsh et al. 2014). However, the assumption of isothermicity is not valid generally (Zucconi et al. 2001; WardThompson et al. 2002; Pagani et al. 2004; Crapsi et al. 2007; Juvela & Ysard 2011). To accommodate a radial temperature profile, a nonisothermal version of the BES (referred here to as a modified BonnorEbert sphere; MBES) has been studied in the literature (Evans et al. 2001; Galli et al. 2002; Keto & Field 2005; Sipilä et al. 2011).
In previous studies of the MBES, it is either assumed that the dust and gas temperatures are equal or that the gas temperature has been derived based on some standard abundances for the cooling species. However, both approaches are approximations and do not hold generally. In Sipilä et al. (2011; hereafter Paper I), we studied the stability of MBESs that are deeply embedded in a larger external structure, such as a molecular cloud, corresponding to a high visual extinction mag at the edges of the studied cores. We also assumed that T_{dust} = T_{gas}, which holds well for lowmass cores with high average densities, but is not valid at lower densities where the collisional coupling between gas and dust is weak. In the present paper, we aim to generalize the analysis of Paper I by including a selfconsistent determination of the gas temperature in the stability calculations. This is accomplished by calculating chemical evolution in model cores with a comprehensive gasgrain chemical model (Sipilä 2012; Sipilä et al. 2013), followed by determining the gas temperature with a radiative transfer model at different time steps, taking advantage of the timedependent chemical abundances. In this way, we can study the stability condition not only as a function of the core mass but also as a function of “chemical time”, which is here defined as the chemical age of the core since some initial state (see Sect. 2.2). We can determine whether the results of Paper I are significantly affected when T_{dust} ≠ T_{gas}, especially in highmass cases when the cores have low average densities and consequently weak gasgrain thermal coupling.
The paper is organized as follows. In Sect. 2, we discuss our model calculations in detail. In Sect. 3, we present the results of our calculations, and discuss them in Sect. 4. We present our conclusions in Sect. 5.
2. Method
This section outlines how the calculations are carried out in practice. The chemical model is discussed here only briefly, and we refer the reader to Sipilä (2012) and Sipilä et al. (2013) for a complete description of the model.
2.1. The MBES
In what follows, we discuss the basic properties of the MBES in a rather concise form; A more detailed discussion on the MBES can be found in Paper I.
The MBES differs from the BES in that it is nonisothermal. Assuming the ideal gas equation of state and hydrostatic equilibrium, and that T = T(r), the density distribution of the MBES is given by (1)This equation can be transformed into nondimensional form by making the substitutions In the above, λ is the central density of the core; τ(ξ) = T(ξ) /T_{c}, where T_{c} is the central temperature of the core; β = kT_{c}/ 4πGm, where k and m are the Boltzmann constant and the average molecular mass of the gas (assumed here equal to 2.33 amu), respectively; ξ and ψ(ξ) are dimensionless variables (ξ represents a nondimensional radius). With these substitutions, Eq. (1) transforms to the modified LaneEmden equation (4)We impose the following boundary conditions at the core center: ψ = 0, dψ/ dξ = 0, τ = 1 and dτ/ dξ = 0. These boundary conditions ensure that T = T_{c} and ρ = λ at ξ = 0, i.e., at the center of the core. To solve Eq. (4), one also needs to supply a temperature profile, which has to be determined externally – we discuss the temperature calculations in more detail below. Finally, the density profile of the MBES is given by substituting the solution function ψ and the temperature profile into Eq. (2).
2.2. Determining the stability of the MBES
The nondimensional radius ξ is a free parameter – for a given core mass, there exist a series of core configurations corresponding to different values of ξ. In what follows, the nondimensional radius of each core configuration will be represented by ξ_{out}. It should be noted that the solution to Eq. (4) is unique to each core configuration (see Paper I), so that the solution is applicable in the interval 0 <ξ<ξ_{out} for each value of ξ_{out} (i.e., separately for each core configuration).
As discussed in Paper I, determining the stability of an MBES of given mass is analogous to finding the core configuration which is critically stable. We construct the different core configurations in the same way as in Paper I, i.e., by following an iterative process. To illustrate the process, outlined in Fig. 1, let us fix the core mass and the nondimensional radius ξ_{out}. We first construct a BES corresponding to the chosen values of core mass and ξ_{out}, and adopting T = 10 K. After this step, we determine a dust temperature profile for the BES using radiative transfer modeling (Juvela & Padoan 2003; Juvela 2005). The dust temperature profile is then used to solve Eq. (4). In Paper I, we considered two different dust models (from Ossenkopf & Henning 1994 and Li & Draine 2001), but here we only consider the former. This issue is discussed in Sect. 4.1.
Fig. 1 Steps taken in deriving the critically stable configuration of an MBES of given mass. See text for detailed explanations of each step. 

Open with DEXTER 
The solution to Eq. (4) and the dust temperature profile are used to derive the central density of the MBES from the equation (5)where τ_{out} = T(ξ_{out}) /T_{c}. After this step, the density profile of the MBES is calculated from Eq. (2). To ensure the consistency of the density profile with the temperature profile, we perform additional iterative calculations, calculating a new dust temperature profile and a new density profile in sequence. Both profiles typically reach their final solutions in a few iterations.
Next, we proceed to calculate the gas temperature. In Paper I, we made the simplifying assumption that T_{dust} = T_{gas}. However, this is a rather crude estimate for cores with total hydrogen density n_{H} ≲ 10^{5} cm^{3}, where the coupling between the dust and the gas is weak. In this paper, we carry out an updated analysis of the results of Paper I, replacing the dust temperature with the gas temperature in all of the relevant formulae. The gas temperature is solved identically to the method discussed in Sipilä (2012). That is, we divide the firstapproximation MBES density profile into concentric shells, solve the chemical evolution separately in each shell and extract chemical abundance profiles for the cooling species (see below) as functions of radial distance from the core center. The chemical model is adopted from Sipilä et al. (2013); the model includes a full description of gasgrain chemistry and the deuterated forms of chemical species with up to four atoms.
The gas temperature is calculated with a Monte Carlo radiative transfer program (Juvela 1997) which balances heating and cooling functions to solve the gas temperature. In the present model, the gas is heated externally by cosmic rays and by the photoelectric effect (Goldsmith 2001). However, the latter process is only important at A_{V} ≲ 2 mag (Juvela & Ysard 2011). Molecular line cooling is calculated for the following species: ^{12}CO, ^{13}CO, C^{18}O, C, O and O_{2}. Since the adopted chemical model does not include the isotopes of the various species, we adopt the isotopic ratios ^{12}CO /^{13}CO = 60 and ^{12}CO/C^{18}O = 500 (Wilson & Rood 1994). We also include the energy exchange (collisional coupling) between the gas and the dust grains following the description of Goldsmith (2001)^{1}. We implicitly assume that chemical timescales are longer than the radiative and dynamical timescales, so that a core is able to quickly reach a new equilibrium state as the chemistry evolves.
After the gas temperature profile has been calculated, it is used to solve Eq. (4) again, producing a new density profile which is then used to redetermine successively the dust temperature, the chemical abundances and the gas temperature as outlined above. The iteration is repeated a few times until the density and the (gas and dust) temperature profiles converge toward their respective final solutions.
Carrying out the above iterative process for different values of ξ_{out}, while keeping the core mass constant, we obtain β, λ and ψ as functions of ξ_{out} and use these to construct the pressure derivative (6)see Paper I for details on how this expression is derived. Finally, the critically stable configuration corresponds to the lowest value of ξ_{out} for which the pressure derivative is zero – all values of ξ_{out} above this value correspond to unstable configurations (Bonnor 1956; Paper I). In the next section, we present the critical ξ_{out} values for a range of core masses derived as outlined above, and we label the critical ξ_{out} value for each core mass as ξ_{1}. Similarly, we label the values of β, λ and ψ(ξ_{out}) corresponding to critical cores as β_{1}, λ_{1} and ψ_{1}, respectively.
Chemical abundances vary with time, and this influences the gas temperature as 1) atomic species are processed into molecules and 2) the main coolant species CO in its various isotopic forms freezes onto the dust grains. To investigate the effect of varying chemical abundances on the stability of the MBES, we have carried out the above analysis for four different time steps, corresponding to 5 × 10^{4}, 1 × 10^{5}, 5 × 10^{5} or 1 × 10^{6} years of chemical evolution since the initial state of the core (see below). In practice, we calculate ξ_{1} for each core mass four times, extracting the chemical abundances at either 5 × 10^{4}, 1 × 10^{5}, 5 × 10^{5} or 1 × 10^{6} years. The different gas temperature profiles at the four time steps translate to marked differences in density profiles, and also in the critical ξ_{1} values. This issue is discussed in the next section. We note that we assume external visual extinction of mag in the chemical and radiative transfer calculations, so that the model cores are assumed to exist deeply embedded inside larger structures, such as molecular clouds. This choice ensures the possibility of direct comparison of our results with those of Paper I.
In the chemical calculations, we assume that the gas is initially atomic, with the exception of hydrogen which is in molecular form (Sipilä 2012). The choice of the chemical composition of the gas in the beginning of the calculation (labeled as t_{0} in the figures) is a parameter of the model. The four timesteps defined above measure the extent of chemical development with respect to the atomic initial state. We note that an atomic chemical composition may not be consistent with an initial physical structure corresponding to a BES (as is assumed here) because chemical processing of the gas is expected to take place during the formation of the BES itself. However, here we do not attempt to impose constraints on “absolute” core ages, but to investigate if chemical evolution can influence the stability condition. Therefore, the chosen four timesteps are simply representative points along the chemical development track of a given core and its surroundings.
3. Results
Fig. 2 Critical nondimensional radius ξ_{1} (upper lefthand panel), ψ(ξ_{1}) = ψ_{1} (upper middle panel), density contrast λ_{1}/ρ_{out,1} = (T_{out,1}/T_{c,1})e^{ψ1} (upper righthand panel), central temperature T_{c,1} = 4πGmβ_{1}/k (lower left), central density λ_{1}/m (lower middle), and the outer core radius R_{out,1} = β^{1/2}λ^{− 1/2}ξ_{1} (lower right) of each critical core configuration as functions of core mass. The profiles are plotted at four different time steps, given in the legend in the bottom middle panel. The blue dotted lines in the upper panels represent the corresponding critical values for the (isothermal) BonnorEbert sphere (ξ_{0} = 6.45, ψ_{0} = 2.64, and λ_{0}/ρ_{out,0} = 14.1). 

Open with DEXTER 
We present in Fig. 2 the critical nondimensional radius ξ_{1} and other related quantities (see below) as functions of core mass and chemical time, as derived from Eq. (6). The values of ξ_{1} derived here are consistently lower than in Paper I for the OH94 model (their Fig. 4). However, the ξ_{1} values predicted by the two works agree to within 10 % regardless of core mass. Hence we derive similar stability for (deeply embedded) critical MBESs regardless of whether we assume T_{dust} = T_{gas} or T_{dust} ≠ T_{gas}.
Evidently, there is marked variation in ξ_{1} both as a function of core mass and as a function of chemical time. A clear decreasing trend in ξ_{1} is seen from the lowest masses (~0.1 M_{⊙}) up to about 0.75 M_{⊙}. ξ_{1} coincides with the isothermal critical radius at ~0.25 M_{⊙}. For masses above 0.75 M_{⊙}, ξ_{1} is, at early times, constant to an accuracy of ~0.15 regardless of core mass. However, at late times the value of ξ_{1} depends on the core mass. In the range 0.75−2.5 M_{⊙}, ξ_{1} decreases as a function of chemical time while for M ≳ 2.5 M_{⊙} an increase with chemical time is evident. For an isothermal BE sphere, the translation from ξ_{out} to physical quantities is straightforward: a larger ξ_{out} translates directly to a higher density contrast between the center and the edge, while the physical radius depends on ξ_{out}, the central density, and the isothermal temperature. Here, the conversion between the various parameters is less evident because the temperature is also a function of radius, and the temperature profile affects not only the density contrast (Eq. (2)) but also the solution to the LaneEmden equation and the determination of the central density (Eqs. (4) and (5)). To understand the significance of the temporal changes in ξ_{1}, we have to consider the associated changes of λ_{1}, β_{1} (i.e., T_{c,1}), and ψ_{1}.
The central density of a critically stable nonisothermal sphere, λ_{1}, decreases monotonously with increasing mass like in the isothermal case. Chemical evolution starts to modify λ_{1} above 0.75 M_{⊙} where it increases with chemical time. The central temperature T_{c} first increases as a function of mass, reaches a maximum, and then decreases again towards the largest masses. Like λ_{1}, T_{c} shows hardly any temporal variation below 0.75 M_{⊙}. Beyond this point, however, T_{c} increases strongly with chemical time and the maximum shifts towards larger masses. Finally, ψ_{1} decreases with chemical time for masses below ~2.5 M_{⊙}, but increases for the higher masses.
The approximate constancy of ξ_{1}, λ_{1} and β_{1} over chemical time at the lowest masses (< 0.75 M_{⊙}) implies that the outer radii of critically stable lowmass cores are similar regardless of how long the gas has been chemically processed. For higher core masses, the physical size of a critical core always decreases with chemical time, but the effect is much more prominent for the highest masses. This is a consequence of the increasing average temperature caused by chemical evolution; hotter cores are better able to withstand gravity and external pressure, and thus higher densities are required for the cores to collapse.
Fig. 3 Density (black lines, scale on left yaxis) and gas temperature (blue lines, scale on right yaxis) profiles of an MBES with M = 1 M_{⊙}, ξ_{out} = 5 as functions of core radius. The profiles are plotted for four different time steps, labeled in the left panel. Left panel: external A_{V} = 10 mag. Middle panel: external A_{V} = 2 mag. Right panel: external A_{V} = 1 mag. 

Open with DEXTER 
For critically stable lowmass MBESs, the density constrast is always lower than that of a critical BES (~14). For M< 0.75 M_{⊙}, the density contrast increases with chemical time owing to changes in the temperature profile; the temperature at the edge of the cores increases with CO depletion (see also Sect. 4.2; note that the central temperature is not affected owing to the efficient gasgrain thermal coupling at the core center) and the consequent increase in the T_{out}/T_{c} ratio overcompensates the decrease of e^{ψ} due to the decrease of ψ as a function of chemical time. For M> 0.75 M_{⊙}, the density contrast largely follows the changes in ψ. For M> 2.5 M_{⊙}, the density contrast can slightly exceed that of the critical BES (~14).
The described tendencies are controlled by three processes: 1) thermal coupling between gas and dust; 2) molecular line cooling; and 3) depletion of molecules onto dust grains. The gasgrain thermal coupling is efficient at densities above n(H_{2}) ~ 10^{5} cm^{3}. Therefore the lowestmass stable cores which have the highest (central and average) densities are practically unaffected by changes in the chemical composition of the gas, because the temperature profile is determined by the interaction between the dust and the interstellar radiation field. The same is true for the centers of slightly more massive cores up to about 1 M_{⊙}. However, their outer parts are cooled efficiently by molecular line emission, until molecules, in particular CO, start to freeze out. With the diminished cooling, the gas temperature rises.
The most massive MBESs have low central densities and the temperature is determined by dust heating and molecular line cooling throughout. The temperature increases strongly in their central parts owing to CO depletion. For each time step, there is a local maximum in the central temperature (as a function of core mass) which shifts towards higher core masses for longer chemical times. This is because highermass cores have increasingly lower average densities and hence the gasgrain thermal coupling is progressively weaker, allowing the temperature to increase further with CO depletion. However, the depletion timescales are very long in the outer layers where the densities are low, and line cooling can operate there also at late chemical times. The steep temperature gradient results in a large density contrast between the center and edge for the most massive cores (~16 for the critical cores).
Given the relatively straightforward temporal changes in T_{c} and λ_{1}, we can deduce that the function ψ_{1} is mainly responsible for the massdependent behavior of ξ_{1} as a function of chemical time. The changes in the cooling efficiency alter the radial density distribution described by ψ, which leads either to a decrease or an increase of ξ_{1} depending on the core mass.
4. Discussion
4.1. The effect of temperature on the results
We have presented the critical radii of MBESs that are deeply embedded in a parent molecular cloud, corresponding to a high external visual extinction mag. This choice facilitates comparison with the results of Paper I where the same assumption was made. However, as pointed out in Sect. 3, the determination of the temperature may have a strong impact on our results because the solution to Eq. (4) is not unique. To study how our results might vary with different initial assumptions, we have calculated the properties of a 1.0 M_{⊙} core with two low values of (2 mag and 1 mag).
Figure 3 presents the density and temperature profiles of a 1.0 M_{⊙} core with (arbitrarily chosen) ξ_{out} = 5 at different time steps, assuming , 2 or 1 mag. The gas temperatures are similar at the core center regardless of the choice of external A_{V}, because of the high density: even for , the visual extinction at the core center is ~10 mag (we assume n(H_{2}) /A_{V} = 9.4 × 10^{20}; Bohlin et al. 1978). However, there are significant differences between the gas temperature profiles outside the core center depending on the choice of external A_{V}. This is because a decrease in visual extinction increases the heating by the photoelectric effect in the core. Although the photodissociation rate of the major coolant molecule CO is strongly increased as well, the net effect on core cooling is not large because in these circumstances atomic carbon takes over as the main coolant species. The gas temperature at the edge of the core for and t = 10^{6} yr is ~26 K.
Fig. 4 Critical nondimensional radius ξ_{1} (upper left panel) of the 1.0 M_{⊙} MBES as a function of chemical time for different values of , indicated in the figure. 

Open with DEXTER 
Fig. 5 Thermal pressure (lefthand panel), density (middle panel) and gas temperature (righthand panel) of an M = 1.0 M_{⊙} MBES with ξ_{out} = 6 as functions of distance from the core center at four different time steps, labeled in the lefthand panel. The vertical lines in each panel represent the outer radius of the core at the different time steps. 

Open with DEXTER 
Evidently, the central densities are similar in all cases, even though the determination of the central density also depends on the temperature and ψ profiles, and not only on the central temperature (Eq. (5)). The density profiles are also similar for chemical times up to t = t_{0} + 5 × 10^{5} yr. For t = t_{0} + 1 × 10^{6} yr at low , the strongly increased gas temperature changes the mass distribution in the core so that the slope of the density profile becomes somewhat steeper.
The differences in temperature depending on the choice of translate to very different solutions to Eq. (4), and consequently the nondimensional radius of the critical configuration is different in all cases. Figure 4 shows the change of the critical nondimensional radius ξ_{1} of the 1.0 M_{⊙} core as a function of chemical time, assuming , 2 or 1 mag. Evidently, the critical radius increases as decreases. Also, for low , the critical radius tends to increase as a function of chemical time, whereas for mag, there is a slight decreasing trend.
We find that the ξ_{1} values are within ~15% of each other regardless of the choice of , when M = 1.0 M_{⊙}. For lower core masses, we expect the results to be closer together because of the high average densities, which lead to significant temperature differences only at the very edge of the cores. However, for M> 1.0 M_{⊙}, the average density is low, so the A_{V} gradient through the core is shallow, and we expect temperature effects to be more prominent. The increase of ξ_{1} with decreasing A_{V} implies that stable isolated cores can have larger density constrasts (between the core center and the edge) than cores embedded in molecular clouds.
We note that in Paper I, we found similar values for ξ_{1} regardless of the dust model, although the dust temperature given by the model of Li & Draine (2001) is lower than given by the model of Ossenkopf & Henning (1994) (see Paper I). In the present paper, we have not considered the Li & Draine (2001) dust model, as we do not expect our results to be strongly dependent on the dust model. This is motivated, on the one hand, by the results of Paper I and, on the other hand, by the fact that ξ_{1} does not change significantly even if the gas temperature profile is radically different (Fig. 4).
4.2. Core evolution
In the stability calculations presented in this paper, we considered a series of model cores defined by the mass and the nondimensional radius, so that we could derive values for the critical nondimensional radius according to an analytical formula (Eq. (6)). However, the choice of the parameters used to define an MBES is free and when constructing a model core one could, instead of the nondimensional radius, set for example the external pressure or the (dimensional) outer radius. In a realistic scenario, one can expect all of these quantities to change as the core and its surroundings evolve; regardless of how the two required parameters are chosen, implicit assumptions are made simultaneously on the properties of the medium outside the core.
If we consider a fixed nondimensional radius like earlier in this paper, the properties of the core are limited by the external thermal pressure exerted on it. This can be understood by studying the evolution of the medium at the core boundary. In Fig. 5, we plot the thermal pressure, density and gas temperature profiles of an MBES with M = 1.0 M_{⊙} and ξ_{out} = 6. The ξ grid was converted to physical radius R using Eq. (3), and vertical lines marking ξ_{out} = 6 were inserted at each timestep to help estimate the change in R. The gas temperature increases as a function of time owing to CO depletion. For late chemical times with increasing thermal pressures, the configuration defined by ξ_{out} = 6 decreases in size (and increases in average density). We note that the hydrodynamical models of Keto & Field (2005) predict that MBESs can exhibit oscillatory behavior. Expansion as a function of time is not achieved in our models because the gas temperature is nearly always increasing, leading to more efficient compression of the cores at long chemical timescales.
The changes in physical radius, density profile etc. at the different time steps considered here do not necessarily represent evolutionary tracks for the cores, because in a realistic scenario neither the nondimensional radius nor the mass of a core are constant. However, the present study demonstrates that the conditions for the core stability change with chemical evolution. For example, in the case of the 1 M_{⊙} MBES, we expect a core configuration with ξ_{out} ~ 6 to be critically stable if the medium is chemically young, while the same configuration would be unstable if associated with chemically old gas. Nevertheless, when a core has been found, e.g., based on the column density and temperature distributions, to agree with a MBES model, the diagrams presented in Fig. 2 of this paper can be used to estimate its stability.
5. Conclusions
We analyzed the stability of nonisothermal (modified) BonnorEbert spheres with a new model that includes a selfconsistent determination of the gas temperature. We compared our results with those of Sipilä et al. (2011), where it was assumed that T_{dust} = T_{gas}. We found that the critical nondimensional radius ξ_{1} changes with the chemical evolution in the core and its surroundings, especially for cores with M> 1 M_{⊙}. This is because of the depletion of the coolant species (mainly CO) onto grain surfaces, which raises the gas temperature. Therefore, cores that exist in a chemically young environment are expected to have a different stability condition than those that exist in regions of chemically old gas.
The ξ_{1} values derived in the present work are slightly lower than in the models of Sipilä et al. (2011). However, the results of the two works agree to within 10 %. The bulk of our analysis was carried out for deeply embedded cores with external visual extinction mag; test calculations for lower values of yield higher critical radii. In summary, our results indicate that the stability of the modified BonnorEbert sphere is similar regardless of whether one assumes T_{dust} = T_{gas} or T_{dust} ≠ T_{gas} and, by extension, also similar to the stability of the classic isothermal BonnorEbert sphere (see Sipilä et al. 2011) – when one assumes that the core is deeply embedded in a larger structure.
We note that Young et al. (2004) have considered a stronger gasgrain coupling than in the model of Goldsmith (2001), which might somewhat compensate for the depletion of coolant molecules at intermediate densities if adopted in our model.
Acknowledgments
We thank the anonymous referees for helpful comments that improved the paper. O.S. acknowledges financial support from the European Research Council (ERC; project PALs 320620), from the Academy of Finland grant 250741, and from the Department of Physics of the University of Helsinki.
References
 Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bacmann, A., André, P., Puget, J.L., et al. 2000, A&A, 361, 555 [Google Scholar]
 Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Bonnor, W. B. 1956, MNRAS, 116, 351 [NASA ADS] [CrossRef] [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, Z. Astrophysik, 37, 217 [NASA ADS] [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]
 Juvela, M. 1997, A&A, 322, 943 [NASA ADS] [Google Scholar]
 Juvela, M. 2005, A&A, 440, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Juvela, M., & Padoan, P. 2003, A&A, 397, 201 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Juvela, M., & Ysard, N. 2011, ApJ, 739, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Kandori, R., Nakajima, Y., Tamura, M., et al. 2005, AJ, 130, 2166 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., & Field, G. 2005, ApJ, 635, 1151 [NASA ADS] [CrossRef] [Google Scholar]
 Li, A., & Draine, B. T. 2001, ApJ, 554, 778 [NASA ADS] [CrossRef] [Google Scholar]
 Marsh, K. A., Griffin, M. J., Palmeirim, P., et al. 2014, MNRAS, 439, 3683 [NASA ADS] [CrossRef] [Google Scholar]
 Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943 [NASA ADS] [Google Scholar]
 Pagani, L., Bacmann, A., Motte, F., et al. 2004, A&A, 417, 605 [NASA ADS] [CrossRef] [EDP Sciences] [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. 2013, A&A, 554, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 WardThompson, D., André, P., & Kirk, J. M. 2002, MNRAS, 329, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [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]
All Figures
Fig. 1 Steps taken in deriving the critically stable configuration of an MBES of given mass. See text for detailed explanations of each step. 

Open with DEXTER  
In the text 
Fig. 2 Critical nondimensional radius ξ_{1} (upper lefthand panel), ψ(ξ_{1}) = ψ_{1} (upper middle panel), density contrast λ_{1}/ρ_{out,1} = (T_{out,1}/T_{c,1})e^{ψ1} (upper righthand panel), central temperature T_{c,1} = 4πGmβ_{1}/k (lower left), central density λ_{1}/m (lower middle), and the outer core radius R_{out,1} = β^{1/2}λ^{− 1/2}ξ_{1} (lower right) of each critical core configuration as functions of core mass. The profiles are plotted at four different time steps, given in the legend in the bottom middle panel. The blue dotted lines in the upper panels represent the corresponding critical values for the (isothermal) BonnorEbert sphere (ξ_{0} = 6.45, ψ_{0} = 2.64, and λ_{0}/ρ_{out,0} = 14.1). 

Open with DEXTER  
In the text 
Fig. 3 Density (black lines, scale on left yaxis) and gas temperature (blue lines, scale on right yaxis) profiles of an MBES with M = 1 M_{⊙}, ξ_{out} = 5 as functions of core radius. The profiles are plotted for four different time steps, labeled in the left panel. Left panel: external A_{V} = 10 mag. Middle panel: external A_{V} = 2 mag. Right panel: external A_{V} = 1 mag. 

Open with DEXTER  
In the text 
Fig. 4 Critical nondimensional radius ξ_{1} (upper left panel) of the 1.0 M_{⊙} MBES as a function of chemical time for different values of , indicated in the figure. 

Open with DEXTER  
In the text 
Fig. 5 Thermal pressure (lefthand panel), density (middle panel) and gas temperature (righthand panel) of an M = 1.0 M_{⊙} MBES with ξ_{out} = 6 as functions of distance from the core center at four different time steps, labeled in the lefthand panel. The vertical lines in each panel represent the outer radius of the core at the different time steps. 

Open with DEXTER  
In the text 