© ESO, 2011

1. Introduction

Physical conditions prevailing in prestellar cores, i.e. gravitationally bound concentrations of molecular gas, constitute the initial set-up of star formation. Prestellar cores are generally not spherical in shape, but it has been found that a Bonnor-Ebert (BE) sphere (Bonnor 1956; Ebert 1955), i.e. an isothermal gas sphere in hydrostatic equilibrium, can succesfully approximate several prestellar cores (e.g. Bacmann et al. 2000; Alves et al. 2001; Kandori et al. 2005). While the density structure predicted by a BE model can give a good fit to observations, there may be problems with e.g. the pressure required for the core to be stable (see e.g. the review by Bergin & Tafalla 2007, and the references therein). Furthermore, the isothermal nature of the observed cores cannot be ascertained, and it is possible that at least some of the “isothermal” cores possess a significant temperature gradient due to attenuation of starlight (e.g. Zucconi et al. 2001; Ward-Thompson et al. 2002; Pagani et al. 2004; Crapsi et al. 2007). A hydrostatic core with the temperature decreasing inwards is one possible configuration preceding the gravitational collapse, and can thus serve as the starting point for dynamical models.

In light of the above, it is of interest to study what kind of constraints could be imposed on the stability of a non-isothermal version of the Bonnor-Ebert sphere (BES). The non-isothermal (“modified”) Bonnor-Ebert sphere (MBES) has been discussed previously by e.g. Evans et al. (2001), Galli et al. (2002), Keto & Field (2005) and Sipilä et al. (2010). In this paper, we derive an explicit stability condition for the MBES, using a method similar to that of Bonnor (1956) in his derivation of the isothermal stability condition. The non-isothermal stability condition will allow us to calculate the critical radius for the MBES, and this in turn yields information on e.g. the possible values of the density contrast and boundary pressure for critical modified BE spheres.

The paper is structured as follows. In Sect. 2 we discuss the properties of the MBES and present a method of obtaining the non-isothermal solution (ξ, ψ). We also present an explicit stability condition and discuss how the critical radius of an MBES of given mass can be found. In Sect. 3 we present the results of our analysis, i.e. the behavior of the physical characteristics (radius, density and temperature structure) of critical cores as a function of core mass. In Sect. 4 we discuss the implications of the results presented in Sect. 3. In Sect. 5 we present our conclusions. Appendix A includes detailed derivations of some expressions presented in Sect. 2, and Appendix B discusses thermodynamics of the MBES.

2. The critical radii of non-isothermal Bonnor-Ebert spheres

2.1. The modified Bonnor-Ebert sphere

The Bonnor-Ebert sphere (Bonnor 1956; Ebert 1955) is an isothermal hydrostatic gas sphere in hydrostatic equilibrium. By assuming the ideal gas equation of state and using the hydrostatic equilibrium condition, Bonnor (1956) obtained an expression for the density distribution of such a sphere (Eq. (2.3) in Bonnor 1956). The density distribution equation can be solved by expressing the density and radius in terms of non-dimensional variables ξ and ψ (Eq. (2.6) in Bonnor 1956); this leads to the Lane-Emden equation (Eq. (2.8) in Bonnor 1956), which reduces to a pair of first-order differential equations. The solution to the Lane-Emden equation then yields the physical parameters of the sphere by back-substitution.

By considering a radial dependence of the temperature, the above analysis can be generalized to non-isothermal spheres; this leads to a modified Bonnor-Ebert sphere (Evans et al. 2001; Keto & Field 2005). In this case, the density distribution equation takes the form (1)In Sipilä et al. (2010) we discussed a modified version of Bonnor’s substitutions (Eq. (2.6) in Bonnor 1956) taking into account the radial dependence of the temperature: where τ = T/T_c, T = T(r) = T(ξ), T_c is the temperature at the core center, λ is the mass density at the core center and β = kT_c/4πGm, where k is the Boltzmann constant and m is the average molecular mass of the gas (here we assume that the gas consists of H₂ and He, so that m = 2.33 amu). We note that our substitutions are slightly different from those of Keto & Field (2005), who used the temperature at the core edge, T_out, as the reference temperature. We have chosen the central temperature T_c as the reference temperature because it depends less on the assumptions of the conditions outside the core. Now, for example, pressure and volume can be expressed in terms of ξ and ψ as where we used Eqs. (2) and (3) and assumed the ideal gas equation of state p = ρkT/m. We note that the dimension in the above equations is carried by β and λ (and of course by G), as ξ and ψ are dimensionless.

Inserting Eqs. (2) and (3) into Eq. (1), one obtains the non-isothermal (“modified”) Lane-Emden equation (6)Imposing the boundary conditions ψ = 0, dψ/dξ = 0, τ = 1 and dτ/dξ = 0 at the center one can (numerically) integrate Eq. (6) and thus obtain a density profile for the MBES, although the temperature profile must be calculated independently using, e.g., radiative transfer means. We discuss this point in greater detail in the following section.

Finally, we note that by choosing T_c = T = const. (τ = 1) one readily recovers the isothermal versions of Eqs. (1) to (6).

2.2. Solving the non-isothermal Lane-Emden equation

Because of its temperature dependence, Eq. (6) cannot be readily solved. One therefore needs to supply a temperature profile. In order to simulate the physical conditions inside prestellar cores, we have chosen to obtain the temperature profile using radiative transfer modeling of the dust component (Juvela & Padoan 2003; Juvela 2005). We assume in the following that T_gas = T_dust; this assumption is discussed below.

Before a temperature profile can be calculated, however, a first-approximation density profile for the model core is needed. A natural starting point for the calculations is the BES because the solution (ξ, ψ) to the Lane-Emden equation is known.

We first fix the mass of the core. The mass of an isothermal core can be written, using the variables introduced previously, as (7)where β_iso = kT/4πGm. In the above and from this point on, the subscript “out” refers to a value taken at the edge of the core. Next, a value for the non-dimensional radius ξ_out has to be chosen. The choice of ξ_out determines the stability of the BES. The critical value is ξ₀ ~ 6.45; cores with ξ_out < ξ₀ are stable, whereas cores with ξ_out > ξ₀ are unstable (see e.g. Bonnor 1956). Also, low values of ξ_out represent cores with a low central density and fairly shallow density gradients, whereas cores with high ξ_out are centrally dense and present a steep density gradient.

Having fixed M and ξ_out, we proceed to solve the central density λ using Eq. (7). A density profile is then constructed using ρ = λexp(− ψ), the isothermal counterpart of Eq. (2). A dust temperature profile for the core can now be calculated. In the temperature calculation, we assume that the core is embedded in a larger molecular cloud; we set A_V = 10 at the edge of the core. The spectrum of the unattenuated interstellar radiation field is taken from Black (1994). The final calculated temperature depends on the properties of the dust component. We carried out separate temperature calculations assuming two different types of grains, using grain opacity data from Ossenkopf & Henning (1994, henceforth OH94) and from Li & Draine (2001, henceforth LD01); the LD01 extinction curve was slightly modified as described in Sipilä et al. (2010). In the latter case, we consider both silicate and graphite grains, with grain material densities 3.5 g cm^-3 and 2.5 g cm^-3, respectively. From now on, we will call core models with OH94 grains “type 1” models, and core models with LD01 grains “type 2” models. The OH94 and LD01 grain models were chosen for this study because these models describe grains in the dense medium inside prestellar cores, where grain properties and possibly their size distribution are thought to differ from those in the diffuse ISM (e.g. Steinacker et al. 2010; Pagani et al. 2010; Juvela et al. 2011).

In all calculations in this paper, we assume that T_gas = T_dust. This assumption should hold well for low-mass cores and in the central parts of more massive cores (M ~ 4−5   M_⊙; hereafter these will be referred to as “high-mass cores”), but generally in the outer parts of high-mass cores T_gas ≠ T_dust (e.g. Galli et al. 2002; Keto & Field 2005). However, high-mass cores should be approximately isothermal owing to the low average density (see e.g. Keto & Field 2005), and hence we can (qualitatively) expect the high-mass MBES to behave like the BES.

After the temperature calculation, the modified Lane-Emden equation (Eq. (6)) is solved using the determined temperature profile. This yields the function ψ as a function of ξ. The function ψ is not yet self-consistent, however, because the solution was obtained starting from the isothermal core. We proceed by updating the central density of the core using the new function ψ. The mass of a non-isothermal core can be expressed as (8)where now β_noniso = kT_c/4πGm. The mass of the core and the non-dimensional radius ξ_out are kept constant during the iteration, so that the above equation yields a new estimate of the central density λ. A new density profile is then calculated according to Eq. (2). A new temperature profile is resolved, yielding a new solution to Eq. (6). The iteration is continued until the function ψ converges; this happens typically after 3–4 iterations. When the iteration is complete, the physical parameters of the core can be derived using Eqs. (2) to (5).

2.3. The critical radius of the MBES

As discussed in Sect. 2.2, the initial choice of ξ_out determines the density structure of the isothermal core, which one uses as a starting point for constructing the non-isothermal core. Because in the isothermal case the Lane-Emden equation does not depend on temperature, the solution function ψ is always the same regardless of the choice of ξ_out; the solution is, in this sense, universal. In the non-isothermal case the situation is different – each value of ξ_out represents a unique core configuration with a unique temperature structure. This means in particular that the solution function ψ depends on ξ_out.

The stability condition of both the BES and the MBES can be found by constructing the derivative of the boundary pressure in terms of core volume. For a stable core, we expect pressure to increase in a contraction of the core, corresponding to a negative sign of the derivative. The critical point is the lowest value of ξ_out for which the pressure derivative turns positive (δp/δV > 0). As shown by Bonnor (1956), all values of ξ_out beyond this point represent unstable cores. We thus need to solve Eq. (6) for multiple values of ξ_out, and look for the value of ξ_out for which the pressure derivative changes sign. In practice, when the iteration is complete, we extract the values of the parameters β, λ and ψ¹ corresponding to each ξ_out and combine them to form “global” functions β(ξ_out), λ(ξ_out) and ψ(ξ_out). These will be needed to calculate the pressure derivative.

The pressure derivative can be constructed by using variational calculus; our method is analogous to that of Bonnor (1956), with the exception that we have extended the discussion to non-isothermal models. In practice, the extension is carried out by considering variations of three parameters (β, λ and ξ), instead of just the two latter ones needed in the isothermal analysis. In the case of the MBES, the pressure derivative takes the form (see also Appendix A) (9)We insert into the above equation the global functions β(ξ_out), λ(ξ_out) and ψ(ξ_out) and look for the value of ξ_out for which the pressure derivative changes sign; this value is the non-dimensional critical radius ξ₁. We calculated this for a range of core masses; the results of the analysis are presented in Sect. 3.

The core mass and hence the total number of particles is conserved in a contraction of the core, i.e. δN = 0. This leads to a condition between δβ, δλ and δξ_out (see Appendix A), which can be used to eliminate for example δβ in Eq. (9). If this is done, it is straightforward to verify that in the isothermal limit, Eq. (9) reduces to the isothermal pressure derivative (Eq. (2.16) in Bonnor 1956).

3. Results

We have derived the physical parameters (radius, density, temperature) of critical MBESs for a range of core masses using the iterative approach described in Sect. 2. In this section we first compare the radial density distributions within a critical MBES and a critical BES, and then describe how the “global” properties, i.e. the central density and the temperature profile of an MBES depend on the core mass.

Fig. 1 H2 number density (solid lines, left y-axis) of critical 0.25   M⊙ type 1 (OH94, blue) and type 2 (LD01, red) MBESs as a function of radial distance from core center. Also plotted is the density profile of the critical 0.25   M⊙ BES corresponding to the mean temperature of the type 1 MBES (green solid line; see text). The temperature profiles of the respective cores are plotted as dashed lines (right y-axis).

3.1. Non-isothermal vs. isothermal: internal density structure

The density distributions of an MBES and a BES for a given mass are compared in Fig. 1. Here we plot the gas number density (n(H₂), solid lines, left y-axis) as a function of radial distance from the core center for a critical MBES of mass 0.25 M_⊙. The type 1 core is plotted in blue and the type 2 in red. Also plotted (green solid line) is the density profile of a critical BES of the same mass and a temperature that is equal to the average temperature of the type 1 core (cf. Eq. (B.4)). The temperature profiles are plotted with dashed lines, keeping the same color codes as for the densities (the scale is on the right).

Evidently, the higher central temperature of a BES as compared with an MBES also permits a higher central density, and the BES is more compact. A similar effect is seen between type 1 and type 2 cores. The type 1 core is slightly warmer and more centrally concentrated than the type 2 core. The temperature difference between the center and the outer boundary is clearly larger for the type 1 core (about 2.5 K) than for type 2 (about 1 K). This is caused by the different dust models.

3.2. The central densities and temperature profiles of critical non-isothermal cores

Fig. 2 Central H2 number densities of type 1 (OH94, blue solid line) and type 2 (LD01, red solid line) critical MBESs as a function of core mass. The green solid line represents central densities of critical BESs corresponding to the mean temperatures of the type 1 MBESs. Also plotted are the outer number densities of the respective core types (dashed lines).

Figure 2 plots the central number densities of critical MBESs and BESs as functions of the core mass. Also plotted are the MBES number densities at the outer boundary. The overall tendency seen in this figure that the central density decreases steeply with the core mass can be understood from Eqs. (7) and (8), according to which the central density (represented by the parameter λ) is approximately proportional to M^-2 when the non-dimensional radius ξ_out does not change. The slight deviations from the M^-2 power-law are caused by temperature changes and small variations of ξ_out discussed below. The central densities of BESs are higher than those of MBESs in the entire mass range, and the type 1 cores are always denser than the type 2 cores. The ratio of central and outer number densities of the MBESs remains approximately constant, regardless of core mass. There is a small increase in this ratio toward low core masses however; we will discuss this in Sect. 3.3.

Figure 3 plots the temperature structures of type 1 (blue) and type 2 (red) critical MBESs as a function of core mass. The upper bounds of the shaded areas represent the outer temperatures of the cores, the lower bounds represent central temperatures. Also plotted are the mean temperatures of the cores as defined by Eq. (B.4) (solid lines). The temperature difference between the core center and the edge increases toward lower masses. This is caused by the increase in the central density (see Fig. 2), which leads to greater attenuation of external radiation. Low-mass MBESs have noticeable temperature gradients, while cores with higher masses are almost isothermal (type 1 cores still present a temperature gradient of  ~1 K at M = 5   M_⊙). Looking at the mean temperatures, we see that for high-mass cores the mean temperature approaches the outer temperature, indicating that most of the core mass lies in the outer parts of the core. For low core masses the mean temperature approaches the central temperature, and the medium is more centrally concentrated.

Fig. 3 Temperature structures of type 1 (blue) and type 2 (red) critical MBESs as a function of core mass. The upper bounds of the shaded areas represent the outer temperatures of the cores, the lower bounds represent central temperatures. The solid lines indicate mean temperatures (see text).

3.3. The critical non-dimensional radius ξ₁

The non-dimensional outer radius ξ_out characterizes the shape of the density profile, which can be compared with high-precision observations as in the famous case of B68 (Alves et al. 2001). The critical value of ξ_out where an MBES becomes unstable against an increase in the external pressure is denoted here by ξ₁ (see Sect. 2.3). Figure 4 plots the value of ξ₁ calculated according to Eq. (9) as a function of core mass. The blue and red lines correspond to type 1 and type 2 models, respectively. Also plotted in the figure is the isothermal critical value ξ₀ ~ 6.45 (green dashed line), which does not depend on core mass.

Fig. 4 Non-dimensional critical radius ξ1 for type 1 and type 2 models (blue and red lines, respectively) as a function of core mass. Also plotted is the isothermal critical radius ξ0 (green dashed line).

One can notice here that for type 1 cores ξ₁ clearly increases toward lower masses, whereas for type 2 cores ξ₁ is nearly constant and close to the isothermal value ξ₀, except for the wiggle at the low-mass end. The different behavior of ξ₁ in the two cases is caused by the different temperature structures of the two core types (Fig. 3). We first consider type 1 cores. At high core masses the temperature gradient of the critical core is  ~1 K, whereas for low core masses the gradient is  ~3.5 K. For type 2 cores the high-mass temperature gradient is  ~0.3 K, whereas the low-mass gradient is  ~1.0 K. Comparing these results with Fig. 4 we see that the critical radius correlates fairly well with the behavior of the temperature gradients; for type 1 models ξ₁ grows as the temperature gradient grows, and for type 2 models ξ₁ is approximately constant until the temperature gradient starts to grow toward low core masses. A more general conclusion is that low-mass cores with a high central density and steep temperature gradient are able to remain stable for slightly higher values of ξ₁.

The increasing temperature gradient toward lower masses is also reflected in the density contrast, (10)For MBESs the critical density contrast is generally close to that of a critical BES (~14), but exceeds the isothermal value for the lowest masses.

We note that the values of ξ₁ plotted in Fig. 4 are subject to numerical uncertainty, mostly because of Monte Carlo fluctuation in the temperature calculations, which can be up to 0.1 K. This fluctuation is reflected in the pressure derivative (Eq. (9)) through its dependence on β, λ and ψ and their numerical derivatives. To analyze the effect of numerical fluctuation, we performed the iterative calculations presented in Sect. 2 with multiple ξ_out grids of varying accuracy (meaning specifically the spacing between neighboring values of ξ_out). We found that the values of ξ₁ for a given core mass can change by up to 0.1 depending on the grid used, but the general shapes of the lines plotted in Fig. 4 remain more or less the same in all cases. The numerical fluctuation is probably responsible for the small dip in the type 2 data present in Fig. 4. Furthermore, we would expect ξ₁ to approach ξ₀ in the high core mass limit (nearly isothermal cores) – the apparent lack of convergence is here attributed to numerical uncertainty.

Galli et al. (2002) studied the structure and stability of MBESs by varying the external pressure and the strength of the ISRF (assuming outside shielding corresponding to A_V = 1). They report different values for the density contrast than the ones presented here. In the particular case of M = 5   M_⊙, they derive ρ_c/ρ_out ~ 20. This is an interesting result, because our models predict nearly isothermal values for the density contrast for high core masses. We attempted to reproduce their results by considering the density profile for the M = 5   M_⊙ sphere as calculated here, but calculating the dust temperature assuming A_V = 1. We calculated the gas temperature in the same way as in Galli et al. (2002), i.e. by balancing heating and cooling functions using the parametrization of Goldsmith (2001). We calculate T_gas ~ 12 K at the core center, decreasing to  ~9 K at about half radius. At the core edge, the temperature rises again to above  ~9.5 K. Using the gas temperature profile to solve the modified Lane-Emden equation (Eq. (6)) and solving the density contrast from Eq. (10) yields ρ_c/ρ_out ~ 20.

The large density contrast is probably caused by the increased temperature at the core edge, where the increased pressure stabilizes the core. The parametrization of Goldsmith (2001) does not take into account the geometry at the edge of the core; in a more “realistic” scenario, increased photon escape probability at the core edge may allow gas temperatures to keep decreasing toward the core edge unless the cloud is so weakly shielded that the external UV field is able to heat its surface layers (through photoelectric effect). This scenario has recently been studied by Juvela & Ysard (2011). In the case of lower gas temperatures at the core edge, the predicted density contrast is probably closer to the isothermal value. However, proving this statement would certainly require a quantitative study.

Finally we note that while there is a dependence on core mass, our models predict that the values of ξ₁ are equal to the isothermal critical value ξ₀ ~ 6.45 up to an accuracy of  ~5%. The density contrasts of the critical MBESs are also very similar to the density contrast of the critical BES,  ~14.

Fig. 5 Physical radii of type 1 (blue solid line) and type 2 (red solid line) critical cores as a function of core mass. The black curve represents homogeneous spheres in virial equilibrium (zero outside pressure) with temperature corresponding to the mean temperature of the 1 M⊙ type 1 MBES. Also plotted are the physical radii of critical isothermal cores corresponding to the mean temperature of each non-isothermal core (dashed lines; see text).

3.4. Physical radii and boundary pressures

The physical radius, R_out, of an MBES of given mass depends on the central temperature and the non-dimensional radius ξ_out according to (11)where c_s,c = kT_c/m is the sound speed in the core center. For critical spheres, ξ_out = ξ₁ is roughly constant (Fig. 4), and changes in the central temperatures are modest, so the radius is roughly proportional to the mass.

We plot in Fig. 5 the physical radii of critical type 1 (blue solid line) and type 2 (red solid line) cores as a function of core mass. The critical radii of the “corresponding” BESs with the temperature , the average temperature of the MBES of the same mass, are plotted with dashed lines. The radii represent the minima of the static equilibrium below which compression leads to gravitational collapse. Also plotted is the physical radius of a virialized homogeneous, isothermal sphere with zero outer pressure, for which . We used the assumption here as well. Evidently the physical radii of MBESs are only slightly larger than those of the corresponding BESs. For the latter the critical radius can be approximated by .

Wo conclude this section with a note about the boundary pressures of the critical MBESs. The boundary pressures p/k of type 1 and type 2 cores as a function of core mass are plotted in Fig. 6 (blue and red lines, respectively). The boundary pressures of the corresponding isothermal spheres are plotted with dashed lines. These diagrams represent the maximum pressures above which an increase in the outer pressure leads to a gravitational collapse. For sub-critical cores, i.e. those with ξ_out < ξ₁, the boundary pressure is lower and the physical radius is larger. The boundary pressures of critical type 2 cores are lower than those of type 1 cores. Low-mass critical BESs have higher boundary pressures than the corresponding MBESs, but the difference decreases toward higher masses and is practically zero for M = 5   M_⊙.

Fig. 6 Boundary pressures p/k of critical type 1 (blue line) and type 2 (red line) MBESs as a function of core mass. Also plotted are the boundary pressures of BESs with temperatures corresponding to the mean temperatures of the MBESs (dashed lines).

We compared our calculations with those of Keto & Field (2005), who presented in their Fig. 1 a p–V curve for an MBES of M = 5   M_⊙. For the thermally supported sphere, Keto & Field (2005) calculate for the critical sphere values of the gas density and boundary pressure of  ~10^4.3 cm^-3 and  ~5.5 × 10³ cm^-3 K, respectively. As also indicated in Figs. 2 and 6, we predict²  ~ 10^4.1 cm^-3 and  ~5 × 10³ cm^-3 K for the central density and boundary pressure, respectively. Our results are thus comparable to those of Keto & Field (2005).

4. Discussion

In this section we discuss the implications of the results presented in the previous section. We make some general remarks about the nature of the stability and discuss the implications of the physical properties of the critical cores on the core chemistry. We conclude the section with a discussion on very small, roughly Jupiter mass cores.

4.1. The nature of the stability

Bonnor (1956) showed by studying the second pressure derivative that the minimum volume for a stable core corresponds to the configuration for which the pressure derivative with respect to volume dp/dV at the core boundary becomes zero. The result applies here as well, although we consider non-isothermal cores. It follows (see also Sect. 3.4) that each physical radius shown in Fig. (5) represents the smallest stable core configuration for a given core mass and model type. That is, for a given core mass, configurations with R > R₁ are stable, whereas all configurations with R < R₁ are unstable.

At this point, however, it is prudent to consider what the stability actually means. The stability condition (Eq. (9)) deals with small first-order fluctuations of the relevant parameters (central density, central temperature and non-dimensional radius). Thus in this context, the statement that a core is stable implies that the core is stable against small fluctuations of these parameters. The possibility that sufficiently large fluctuations could induce collapse even in a “stable” core cannot be ruled out based on the analysis presented here.

Also, our analysis does not give insight into what happens after a core begins to collapse. It may well be (as has been demonstrated in the literature and by our preliminary MHD collapse models) that “unstable” core configurations settle into a stable state relatively soon after collapse begins. This is caused by the increase of the temperature inside the core: as the core collapses, the thermal pressure inside the core increases and may be able to support the slightly contracted core, depending on the core-cooling mechanisms. Indeed, the stability condition presented here provides an initial condition for the collapse of non-isothermal cores, rather than trying to predict the “absolute” stability of the cores.

4.2. Physical structures

Figures 1 to 3 demonstrate that a type 1 MBES has a different physical structure than a type 2 MBES of the same mass. In general, type 1 cores are denser, warmer and smaller (Fig. 5) than type 2 cores, i.e. type 1 cores are more compact. Toward low core masses, the differences in radius and density diminish, but type 1 cores always present a steeper temperature gradient.

If one computes the mean temperature of a given MBES, a “corresponding” BES of the same mass can be constructed. Figure 2 shows that the central density of a BES constructed in this way is higher than that of the MBES – this is also apparent in Fig. 1 in the special case of M = 0.25   M_⊙. The BES is, however, smaller. This result is general for low core masses. However, when using the gas temperature the situation changes to the opposite for higher masses (see e.g. Keto & Field 2005).

One can, of course, construct a BES with outer radius R exactly equal to that of a given MBES. In this case, the temperature of the BES is clearly lower than the mean temperature of the MBES, but the density structures turn out to be nearly identical, with very small differences caused by the different temperature structures.

4.3. Chemistry

We now briefly discuss the impact of the physical properties of the cores on their chemical evolution.

Let us first concentrate on the differences between the type 1 and the type 2 MBES. As discussed above, a type 2 MBES is larger than a type 1 MBES of the same mass, although the former is less dense and colder. The typical density difference in the center is less than 15% between the core types, however, so in this case the major factor in the possible different chemical evolution of the cores is the temperature gradient.

Even small temperature changes can have significant consequences for the chemistry and the interpretation of observations. Aikawa et al. (2005) studied the chemical evolution in collapsing clouds with initial conditions close to critical Bonnor-Ebert spheres. Comparing models with a central density of 3 × 10⁶ cm^-3 and temperatures in the range 10 K–15 K, the molecular depletion in the central parts of the core was found to be very sensitive to the temperature. In particular, the abundances of NH₃ and N₂H⁺ could differ by an order of magnitude or more. However, because the central abundances are always low compared to the outer regions, it may be impossible to observationally confirm this temperature dependence.

Let us then turn our attention to the differences between the MBES and the BES. As discussed in Sect. 4.2, the BES is always more compact than the corresponding MBES. The temperature dependence is more complex than when comparing two MBESs – the BES is warmer in the center, but colder toward the edge. These differences may again produce notable differences in chemical evolution.

It is difficult to determine core temperature profiles from observations and, similarly, theoretical predictions contain significant sources of uncertainty. Goldsmith (2001) showed that molecular depletion can significantly increase the gas temperature predicted for cores. The effect can be several degrees even in dense clouds although, once the density is  ~10⁵ cm^-3 or above, T_gas tends toward the dust temperature irrespective of the decreasing line cooling power. Juvela & Ysard (2011) demonstrated further how, in addition to abundance variations, the radial temperature profiles can be modified by the velocity field (affecting the efficiency of line cooling) and the grain size distribution (affecting the gas-dust coupling). Therefore, each hydrostatic object is likely to exhibit slightly different density profiles.

The above implies that chemistry plays a role in the dynamical evolution of a core as well. A study of the chemistry of the MBES including line cooling is planned – this will allow a more realistic estimate of the gas temperature, also serving as a starting point for dynamical studies. In this context, one can also quantitatively study how the gas temperature affects the stability analysis presented here. Finally we note that the stability analysis presented in this paper indicates that an MBES can be just as stable as a BES of the same mass and (roughly) the same size, making it possible for advanced chemistries to develop in both cases in a typical core lifetime.

4.4. Jupiter mass non-isothermal cores

In the above, we presented results of a stability analysis for core masses ranging from 0.1   M_⊙ to 5.0   M_⊙. It has been suggested that small dark clouds of roughly Jupiter mass could make up for part of the submm sources observed in e.g. SCUBA maps, or even account for part of the dark matter (Lawrence 2001). Although this scheme has later been found to be unlikely (see e.g. Drake & Cook 2003; Almaini et al. 2005), the existence of such objects is not entirely ruled out. For this reason, we thought it interesting to model these roughly Jupiter mass clouds (JMC) as either a BES or an MBES to find out what kind of constraints would be imposed on the stability of these objects.

In our calculations thus far we assumed A_V = 10 at the core edge, so that the objects are embedded in a bigger construct, e.g. a molecular cloud. However, if the JMCs are thought to be isolated, then the assumption of A_V = 10 at the core edge is no longer justified. Indeed, if a JMC is isolated, the outermost parts are most likely ionized and heated to high temperatures by the ISRF, and molecular gas can be thought to be found at a minimum visual extinction of A_V ~ 2. We have thus carried out the stability analysis for a core mass of 0.01   M_⊙ (roughly ten Jupiter masses), assuming outside A_V = 2.

As a point of reference, let us first look at the BES. The critical radius of a 10 M_Jup BES at 7 K is 140 AU, and the corresponding boundary pressure is p/k = 4 × 10⁸ cm^-3 K. As discussed in Sect. 3, these values represent the minimum physical radius and maximum boundary pressure for a core in stable equilibrium. The critical radius is directly proportional to the mass and the critical boundary pressure is proportional M^-2. So for example, a 1 M_Jup BES at 7 K would have a critical radius of 14 AU and a boundary pressure of 4 × 10¹⁰ cm^-3 K. According to Lawrence (2001), the blank-field submm sources have angular radii of about 1″ or less, and if interpreted as cold dark clouds, they should be nearby objects with a characteristic distance of 100 pc. Hence, a typical radius for this kind of cloud should be about 100 AU. This radius roughly equals the critical radius of a 10 M_Jup BES. For a 1 M_Jup BES at 7 K, the radius R = 100 AU corresponds to a sub-critical dimensionless radius of ξ_out ~ 1. The boundary pressure for this configuration is still very high, p/k ~ 2 × 10⁸ cm^-3 K.

Let us then look at the MBES. For a critical type 1 core, the physical radius is  ~200 AU. The temperature gradient is fairly steep (from 5 K at the center to 10.5 K at the edge), and the boundary pressure is p/k ~ 10⁸ cm^-3 K. The temperature gradient of the type 2 core is shallower (from  ~4.8 K at the center to  ~8.5 K at the edge), but R and p/k are nearly identical to those of the type 1 core. This analysis suggests that the minimum size of a stable 10 M_jup MBES is  ~200 AU, irrespective of the dust model used. However, the boundary pressure required to maintain the critical configuration in equilibrium is unlikely to be present in interstellar conditions, at least if the core is isolated. For clearly subcritical configurations (ξ < ξ₁, R ~ 1000 AU), the boundary pressure is still  >10⁷ cm^-3 K.

We conclude that owing to the high boundary pressures required, roughly Jupiter mass critical BESs or MBESs are unlikely to be able to exist as isolated objects in interstellar space where the pressure is typically  ~10⁴−10⁵ cm^-3 K (e.g. McKee & Ostriker 2007). Subcritical configurations can have radii of  ~1000 AU, but the boundary pressure required for stable equilibrium is very high for these as well.

5. Conclusions

We studied the stability of non-isothermal Bonnor-Ebert spheres. The physical parameters of critical cores were derived for a range of core masses. Two different types of dust grains, corresponding to optical data from Ossenkopf & Henning (1994) and Li & Draine (2001), were considered in the modeling. The analysis was accordingly separated into two distinct core types (type 1 and type 2, respectively), each containing one type of dust grains.

As a general trend, the central density of a core increases as core mass decreases. The increase in density also increases the temperature gradient because of more efficient attenuation of the external radiation field (in this paper, we studied the dust temperature only). However, the absolute values of the temperature depend on the adopted dust model. Type 1 cores present steeper temperature gradients than type 2 cores, but type 2 cores are colder. The difference also translates to core size: owing to the higher temperature, the thermal pressure in the centers of type 1 cores is higher and hence a type 1 core can assume a smaller size than a type 2 core of the same mass. Thus, type 1 cores are more compact. Isothermal cores corresponding to non-isothermal cores (i.e. with same mass and same average temperature) are more centrally dense and slightly smaller than their non-isothermal counterparts. Considering gas temperature instead of dust temperature may affect core stability (see e.g. Keto & Field 2005). Our results should hold well at least for cores with masses up to  ~0.5   M_⊙, however, for which the average density is  ~10⁵ cm^-3.

It was found that the non-isothermal critical radius ξ₁ increases above the isothermal critical value ξ₀ ~ 6.45 toward low core masses. The effect is particularly pronounced for type 1 cores, which present steep temperature gradients. Nevertheless, the change in ξ₁ with core mass is fairly small, and in the mass range studied here, ξ₁ is equal to ξ₀ to within 5%. Furthermore, the ratio of central and outer densities for critical MBESs was found to present a similar increase toward low core masses, but in this case as well the ratio is close to the isothermal value (~14).

Although the differences in physical parameters between the two MBES types are not great, the different temperature gradients in particular may affect chemical evolution in these objects, which might be observable through line emission radiation from these objects. This may also be the case when comparing isothermal and non-isothermal cores. Line radiation cooling can also affect the stability of the cores. A quantitative study of chemical evolution in the different types of cores could consequently be justified.

We studied the physical parameters of MBESs of roughly ten Jupiter masses (~0.001−0.01   M_⊙) to validate or disqualify their possible existence as isolated objects. We found that the boundary pressure required to maintain critical cores equilibrium is unlikely to be found in interstellar conditions. Even subcritical configurations require boundary pressures of p/k ~ 10⁷⁻⁸ cm^-3 K. We conclude that very low mass MBESs are unlikely to be able to exist in the interstellar medium as isolated objects – this applies to very low mass BESs as well.

Finally, we note that the stability analysis presented in this paper considers small, first-order variations of the relevant parameters (central density and temperature, core radius). Thus, stable cores as defined by the stability condition derived here are stable against linear perturbations, and sufficiently large perturbations could induce collapse even in these “stable” cores. Furthermore, what happens after collapse begins is not predicted by the equations. However, the near-critical gas spheres studied here represent plausible, albeit idealistic models for cores at the very beginning of collapse.

Acknowledgments

O.S. acknowledges support from the Väisälä Foundation of the Finnish Academy of Science and Letters. The study has also been funded by the Academy of Finland through grants 132291 and 127015. The authors thank the referee Dr. Daniele Galli for helpful comments, which improved the paper.

References

Appendix A: The stability condition

In this appendix we present a short derivation of the general form of the derivative of the boundary pressure with respect to core volume (Eq. (9)). We also present an expression for the total number of particles in terms of β, λ and ξ_out. Throughout the calculations, we make use of Eqs. (2) to (6) in the main text.

A.1. The pressure derivative

Calculating the variation of the pressure (Eq. (4)) with respect to the three non-dimensional variables β, λ and ξ_out yields (A.1)The variation of the cloud volume is (A.2)The above equations then yield the pressure derivative in general form: (A.3)Identifying leads finally to Eq. (9).

A.2. Mass conservation

We now present a non-isothermal version of the mass conservation equation discussed by Bonnor (1956). Using Eqs. (2) and (3), we can write the total number of particles as (A.4)the mass of the core is M = Nm. Because N is calculated for a single core (that is, for a given value of ξ_out), one can move β and λ out of the integral; for any single core, these are constant. In the integral, the variable ξ represents the internal structure of the core and therefore we make use of the solution to the Lane-Emden equation (Eq. (6)) that is unique to each core. The above formula simplifies to (A.5)where τ_out represents the temperature contrast of the core, that is the ratio of outer and central temperatures. During the collapse process we expect the total number of particles (i.e. the total mass of the core) to be conserved, i.e. that δN = 0. Calculation of this condition leads to the equation (A.6)It is readily verified that in the isothermal limit (δβ = 0, τ = 1), Eq. (A.6) reduces to its isothermal analogue (Eq. (2.11) in Bonnor 1956).

Appendix B: Thermodynamics of the MBES

Bonnor (1956) discussed the form of Boyle’s Law for a massive isothermal gas sphere. The radial change of the temperature causes a slight modification to Bonnor’s equation of state, and we think there is a reason to make a brief excursion to the thermodynamics of an MBES. Besides writing down the equation of state, we also present the first law of thermodynamics for an MBES, which can be useful for dynamical calculations or when discussing thermal instability (see Keto & Field 2005).

It turns out that the internal energy of an MBES can be described by the simple ideal gas expression of internal energy when the temperature T is replaced by the mass-averaged temperature  defined below. In general terms, an MBES can be considered as a thermodynamic system containing N molecules in a volume V at temperature . The thermodynamic potentials and the equation of state can be derived much in the same way as for an ideal gas (see e.g. Landau & Lifshitz 1969; Mandl 1988) from a partition function Z, which is the product of the perfect gas partition function, (which may also contain rotational and vibrational parts), and the configurational partition function , where Ω is the gravitational potential energy. For example, the Helmholtz free energy is then , and the entropy can be obtained from . Below, we use a more straighforward method in the derivation of the internal energy. In any case, Ω has to be determined from the density distribution, and in this prescription it is a function of N (mass), V (outer radius), and an additional variable ξ_out that describes the distribution of mass within the volume V.

B.1. Internal energy U

Assuming that the medium consists of ideal gas with effective degrees of freedom³f = 3, we can equate internal energy with total thermal energy (see e.g. Chandrasekhar 1957; Kippenhahn & Weigert 1994) and write (B.1)Using Eqs. (2) to (4), we can write the integral in terms of β, λ, and ξ, and obtain (B.2)On the other hand, the assumption of ideal gas means that the total internal energy of an isothermal core is (B.3)It can be shown that Eq. (B.3) is valid also for a non-isothermal sphere if T is replaced with the average temperature defined as (B.4)When substituting the expressions of and N from Eqs. (B.4) and (A.5) to Eq. (B.3) one obtains the formula for U given in Eq. (B.2).

B.2. Gravitational potential energy Ω and the virial theorem

The total potential energy of a gas sphere of radius R can be written as (B.5)where ρ = ρ(r) and we have written the mass of a spherical shell as dM = 4πr²ρdr. Using the hydrostatic equilibrium condition (B.6)we can rewrite Eq. (B.5) into an easily integrable form: (B.7)The latter term in Eq. (B.7) can be identified as twice the thermal energy (see Eq. (B.1)). Rearranging terms yields the virial theorem (B.8)Writing Ω in terms of β, λ and ξ_out and using Eqs. (3) and (8), it can be shown that (B.9)where f is a function of ξ_out only. To illustrate, we plot in Fig. (B.1) the ratio f(ξ_out) = Ω/(− GM²/R) for a BES and an MBES of M = 1   M_⊙. Also marked is the isothermal critical value ξ₀ ~ 6.45. The figure shows that the gravitational potential energies of the BES and the MBES differ from each other somewhat. For low values of ξ_out the cloud is nearly homogenous and f(ξ_out) approaches the value 0.6 = 3/5.

Fig. B.1 Function f(ξout) for a BES and an MBES of M = 1   M⊙.

B.3. Equation of state

Substituting the expression of the internal energy in terms of the average temperature  to the virial theorem (Eq. (B.8)), one obtains the equation of state at the core boundary (B.10)Using Eq. (B.9), this can be written in the form suggested by Terletsky (1952, Bonnor’s Eq. (1.2)), criticized by Bonnor, if the numeric factor α is taken to be a function of ξ_out.

B.4. First law of thermodynamics

Considering the internal energy only, the first law of thermodynamics states that (B.11)where the last term represents the work done by the gravitational potential energy. For a reversible process, the first term can replaced by heat dQ supplied to the system. From Eq. (B.9) it can be seen that (B.12)Furthermore, using the equation of state, the change of the internal energy becomes (B.13)The validity of this formula can be verified by using on the right-hand side one of the Maxwell relations (B.14)On the other hand, from the definition of heat capacity: (B.15)where we used the result that for a monatomic ideal gas . Equations (B.14) and (B.15) together yield (B.16)Substituting this into the right-hand side of Eq. (B.13), one obtains the identity (B.17)in accordance with Eq. (B.3).

The total energy of the core equals to the sum of internal energy and gravitational potential energy, E = U + Ω, from which follows that (B.18)According to Eq. (B.9), we can write Ω = Ω(V,ξ_out) when the mass is kept constant. Then (B.19)Inserting Eqs. (B.11) and (B.19) into Eq. (B.18) yields (B.20)\newpage\noindentUsing Eq. (B.16) above the change of the total energy takes the form (B.21)The above equation is another expression of the first law of thermodynamics.

All Figures

Fig. 1 H2 number density (solid lines, left y-axis) of critical 0.25   M⊙ type 1 (OH94, blue) and type 2 (LD01, red) MBESs as a function of radial distance from core center. Also plotted is the density profile of the critical 0.25   M⊙ BES corresponding to the mean temperature of the type 1 MBES (green solid line; see text). The temperature profiles of the respective cores are plotted as dashed lines (right y-axis).

In the text

	Fig. 2 Central H2 number densities of type 1 (OH94, blue solid line) and type 2 (LD01, red solid line) critical MBESs as a function of core mass. The green solid line represents central densities of critical BESs corresponding to the mean temperatures of the type 1 MBESs. Also plotted are the outer number densities of the respective core types (dashed lines).
In the text

	Fig. 3 Temperature structures of type 1 (blue) and type 2 (red) critical MBESs as a function of core mass. The upper bounds of the shaded areas represent the outer temperatures of the cores, the lower bounds represent central temperatures. The solid lines indicate mean temperatures (see text).
In the text

	Fig. 4 Non-dimensional critical radius ξ1 for type 1 and type 2 models (blue and red lines, respectively) as a function of core mass. Also plotted is the isothermal critical radius ξ0 (green dashed line).
In the text

Fig. 5 Physical radii of type 1 (blue solid line) and type 2 (red solid line) critical cores as a function of core mass. The black curve represents homogeneous spheres in virial equilibrium (zero outside pressure) with temperature corresponding to the mean temperature of the 1 M⊙ type 1 MBES. Also plotted are the physical radii of critical isothermal cores corresponding to the mean temperature of each non-isothermal core (dashed lines; see text).

In the text

	Fig. 6 Boundary pressures p/k of critical type 1 (blue line) and type 2 (red line) MBESs as a function of core mass. Also plotted are the boundary pressures of BESs with temperatures corresponding to the mean temperatures of the MBESs (dashed lines).
In the text

	Fig. B.1 Function f(ξout) for a BES and an MBES of M = 1   M⊙.
In the text