1.Introduction

It is widely thought that the gravitational collapse of dense molecular cloud cores results in the formation of stars (see, e.g., Shu 1977). Detailed calculations of the thermal structure of molecular cores show that over a large radial range they have only small temperature variations (see, e.g., Falgarone & Puget 1985), so that their hydrostatic structure can be modeled as an isothermal, self-gravitating, sphere.

In this paper, we study self-gravitating spheres with a polytropic equation of state ($P\propto {\rho }^{\gamma }$, where P is the pressure, p the density and γ the specific heat ratio of the gas) with 𝛾 values close to 1, meant as models for cloud cores with small departures from isothermality. We first compute the density profiles of these hydrostatic structures, and then evaluate their stability to radial perturbations.

We use the method of Bonnor (1956) for determining the stability of the isothermal, self-gravitating sphere. Bonnor (1956) considers a sphere that has a fixed mass M and an outer radius R_e at which the self-gravitating structure is in contact with a much hotter environment with a pressure P_e. He then calculates the variation in the outer pressure (or density) of this fixed mass sphere as a function of a changing outer radius R_e (assuming that the inner structure adjusts instantaneously to the steady, self-gravitating configuration). If the outer density p_e (or, alternatively the pressure P_e, which is a monotonically increasing function of p_e) satisfies the condition

the sphere will be stable. In the case of the isothermal sphere, the sphere is stable for small radii, and becomes unstable (i.e., $d{\rho }_e/d{R}_e$ becomes positive) at a radius R_S. Isothermal spheres with R_e > R_S are therefore unstable to radial perturbations.

Bonnor (1958) extended his analysis to the case of self-gravitating spheres of arbitrary polytropic index $n=1/(\gamma -1)$. He derived (from the numerical density stratifications of Emden, 1907) the values of the maximum stable radius R_e for the cases with n = 4, 4.5, 4.9, 5 and 6. Bonnor’s criterion has also been applied to hydrostatic structures with thermal properties calculated using an energy equation with appropriate heat source and sink terms (see Falgarone & Puget 1985).

Bonnor (1958) finds that self-gravitating, polytropic spheres with γ > 4/3 are inconditionally stable to radial perturbations. Because of this, in the present paper we limit ourselves to polytropes with $1\le \gamma <4/3$. The stability of polytropes with different 𝛾 values has been described more recently by Horedt (2013).

An interesting problem is whether or not Bonnor’s stability criterion gives “correct” results, where by “correct” we mean reproducible by time-dependent, hydrodynamic collapse simulations. Until quite recently, hydrodynamic simulations of the collapse of isothermal spheres did not reproduce (at least, not in a quantitative way) the results found from Bonnor’s criterion. Hunter (2001) speculated that this lack of agreement was probably the result of a lack of accuracy in the hydrodynamic simulations. More recently, Raga et al. (2013) computed more precise simulations which reproduced the stability found from Bonnor’s (1956) criterion for the collapse of the isothermal sphere. Even though it has not yet been proven, in this paper we will assume that Bonnor’s (1958) stability criterion for arbitrary 𝛾 also predicts successfully the stability limits found from time-dependent hydrodynamical collapse calculations.

The paper is organized as follows. In § 2 we present the Lane-Emden equation and give a summary of the analytic properties of its solutions for specific heat ratios in the ECUACION range. In § 3 we describe the density and mass distributions obtained from numerical integrations of the Lane-Emden equation. In § 4, we apply Bonnor’s (1958)) stability criterion to these solutions, and derive the properties of the maximal stable sphere as a function of γ. Finally the results are summarized in § 5.

2.The Lane-Emden Equation

2.1.General Considerations

The hydrostatic equation for a self-gravitating, non-magnetized spherical cloud is:

where R is the spherical radius, P the gas pressure, p the density and G the gravitational constant.

We now assume that we have a polytropic law

with constant C and γ. Then, equation ([hyd]) takes the form

which receives the name of “Lane-Emden equation”.

As we are thinking of modelling a gas cloud, possible values of γ are 5/3 (for an atomic/ionic, non-radiative gas), 7/5 (for diatomic molecules) and 1 (for an isothermal situation). Another interesting value is 4/3, corresponding to a gas of ultra-relativistic fermions (relevant for degenerate stellar cores).

In a stratified cloud, there is in principle no reason for C in equation (3) to be independent of radius (individual fluid parcels will evolve satisfying equation 3, but do not necessarily share the same value of C). A unique value of C is obtained for a “well mixed” gas cloud, in which mixing due to slow turbulent motions produces a stationary stratification with a constant specific entropy. The existence of such a mixing is required for modelling gas clouds with the polytropic Lane-Emden equation (equation 4).

Finally, we should note that many times the solutions to Lane-Emden’s equations are identified with the “polytropic index” $n=1/(\gamma -1)$ (which appears naturally under some traditional changes of variables used for this equation). However, we will use γ (and not n) throughout the present paper.

2.2.The “Small R”, Non-Singular Solution

We assume that we have a non-singular, “flat-topped” solution, with $(d\rho /dR{)}_0=0$ in the cloud center. Proposing a second-order Taylor series for the p as a function of R, inserting it in equation (4) and equating terms with the same powers of R one straightforwardly obtains:

where p₀ is the central density and

is the “core radius” of the density stratification.

For the variables $r=R/R_c$ and $\rho \text{'}=\rho /{\rho }_0$ the Lane-Emden equation takes the dimensionless form:

with second-order solution $\rho \text{'}=1-r^2$.

It is possible to construct higher order solutions by proposing appropriate Taylor expansions for the density as a function of radius. However, it is a well known problem that for γ ≈ 1 the resulting series solutions converge to the exact solution only for (see, e.g., Chandrasekhar 1967).

We should note that the core radius R_c used for the adimensionalization (see equations 6 and 7) has a clear dependence on γ. In order to compare solutions with different γ values, it is possible to choose solutions with the same central isothermal sound speed:

where the first equality is the definition of the isothermal sound speed and the second equality is obtained using equation (3). Substituting into equation (6) we then obtain

Where R_c,1 is the core radius of the isothermal (γ = 1) solution with the same central sound speed c₀.

Therefore, for comparing solutions with different γ values but with the same central temperature (corresponding to a sound speed c₀), we plot the solutions as a function of

see equations (7) and (9).

The mass of the self-gravitating spheres as a function of the spherical radius is given by:

where

is the dimensionless mass. It is straightforward to see that ${\gamma }^{3/2}m\left(r\right)$ is the mass normalized to $4\pi {\rho }_0{R}_{c,1}^2$ (where R_c,1 is the core radius of the γ =1 solution, see above).

2.3.The Singular Solutions

In order to obtain the singular solutions, one inserts a proposed solution ${\rho }_{sing}=A/R^p$ into the Lane-Emden equation (4). Equating the resulting powers of R and multiplying constants one obtains:

Introducing a similar power law into the dimensionless Lane-Emden equation (7) one obtains:

Clearly, these singular solutions exist only for γ < 4/3.

From numerical solutions of the Lane-Emden equation, it is found that for γ < 1.2 the non-singular solution (with the small r solution $\rho \text{'}=1-r^2$ and large r regime obtained by numerically integrating equation 7) converges to the singular solution for r >> 1. For γ = 6/5 = 1.2 (corresponding to a polytropic index n = 5), the Lane-Emden equation has the analytic, non-singular solution

This solution has a large r behaviour of the form $\rho \text{'}\propto r^{-5}$, clearly not matching the singular solution (which for γ = 6/5 has a shallower $\rho \text{'}\propto r^{-5/2}$ dependence, see equation 15). The transition between the large r, singular solution convergence obtained for γ = 6/5 and the γ = 6/5 solution (given by equation 16) is explored in § 3.

For γ > 1.2, the numerical non-singular solution of the Lane-Emden equation has a density that goes to zero at a finite radius. Therefore, the non-singular solution has a $r\gg 1$ convergence to the singular solution only for γ = 6/5. The nature of this convergence is described in the following subsection.

2.4.The γ = 6/5, r >>1 Convergence to the Singular Solution

It is customary to write

and use equation (7) to obtain a differential equation for the fractional deviation q(r) from the singular solution. Assuming that q(r), $\dot{q}=dq/dr$ and $\ddot{q}=d^{2}q/d{r}^2\ll 1$, this differential equation can be linearized to obtain:

which for $a>3-2\sqrt{2}$ has the general solution

where B and φ are arbitrary constants,

and

It is straightforward to see that the $a>3-2\sqrt{2}$ condition (necessary for equation 19 to be valid, see above) is satisfied for the $\gamma =1\to 1.2$ range.

Therefore (for all values of γ in the $1\to 1.2$ range) the large r convergence to the singular solution is oscillatory (going as $\mathrm{c}\mathrm{o}\mathrm{s}\Gamma \mathrm{l}\mathrm{n}r$). This oscillatory behaviour has an amplitude that decreases as r^-a (see equations 19-21).

The linear problem which we have presented here does not allow us to fix the values of the B and φ constants (see equation 19). These constants have to be determined by fitting equation (17), calculated with the linear q solution (see equation 19), to the non-singular solution calculated with a numerical integration of the Lane-Emden equation (see equation 7). This fitting can be done to a quite high level of accuracy (see, e. g., the work of Ito et al. 2018 for the γ =1 solution).

3. Numerical integrations

We have numerically integrated the dimensionless Lane-Emden equation (equation 7) starting from the quadratic solution $\rho \text{'}=1-r^2$ at small r. The results obtained for several values of γ in the $1\to 1.2$ range are shown in log-log plots in Figure 1.

Fig 1
The solid lines show the dimensionless nonsingular density vs. radius solutions for the isothermal, γ =1 case (top), for γ =1.5 center, and for γ = 1.190, 1.199 and 1.200 (bottom). The dashed curves show the γ =1, 1.15 and 1.2 singular solutions (in the top, center and bottom frames, respectively). The color figure can be viewed online.

For r < 1, all of the logarithmic density distributions are quite flat, with $\rho \text{'}\approx 1$. For r > 1, the γ = 1 and 1.5 distributions start a steeper decrease, and approach the corresponding singular solutions (shown with the straight, dashed lines) with the oscillatory convergence described in § 2.3. These oscillations have larger amplitudes for solutions of increasing γ.

The bottom panel of Figure 1 shows solutions with γ ≈ 1.2 values. The singular γ =1.2 solution (not distinguishable from the γ = 1.190 and 1.199 singular solutions at the resolution of the plot) is shown with a dashed line, and the non-singular γ = 1.2 solution (obtained either from the numerical integration or from equation 16) is shown with the red line. It is clear that the γ = 1.190 and γ = 1.199 solutions follow the (analytic) γ = 1.2 solution out to $r_t\approx 10$ and 100 (respectively), and then have an oscillatory transition to the singular solution. We find that the transition radius r_t becomes progressively larger as $\gamma \to 1.2$. Solutions with $\gamma \approx 1.2$ have been previously studied by Horedt (2013).

Figure 2 shows logarithmic depictions of the $\rho \text{'}\left(r\right)$ solutions obtained for γ values in the $6/5\to 4/3$ range. For all of the γ > 1.2 solutions we obtain densities that go to zero at a finite dimensionless radius $r_{max}=R_{max}/R_c$ (which can be straightforwardly determined from the numerical integrations).

Fig. 2
Dimensionless nonsingular density (p/p₀) vs. radius (r=R/R_c , where R_c is the core radius defined in equation 6) solutions for γ =1.20, 1.21, 1.25 and $1.\hat{3}$. While the γ = 1.20 solution goes out to infinite r, the solutions for larger γ values reach zero densities at at an outer radius r_f (which decreases for increasing γ).

Finally, in Figure 3 we show the inner regions of density p/p₀ and mass ${\gamma }^{3/2}m$ (see equation 12) stratifications with y in the $1\to 4/3$ range, as a function of the radius ${\gamma }^{1/2}r=R/R_{c,1}$, normalized to the core radius of the γ = 1 solution (see equations 9-10). The results shown illustrate the differences between self-gravitating spheres with different γ, but with the same central density p₀ and isothermal sound speed c₀. We see that the lower γ solutions have lower densities and masses within the core radius, and higher densities and masses at large radii.

Fig 3
Bottom frame: inner region of the nonsingular density vs. radius solutions for γ = 1.0, 1.1, 1.2. 1.3 and $1.\hat{3}$. The solutions for higher y values have steeper density decreases at larger radii. Top frame: dimensionless mass ${\gamma }^{3/2}m=M/\left({\rho }_0{R}_{c,1}^3\right)$, where p₀ is the central density) as a function of radius. The mass distribution has lower growths at large radii for higher y values. Both the mass and the density profiles are shown as a function of ${\gamma }^{1/2}r=R/R_{c,1}$ (where R_c,1 is the core radius of the γ = 1 solution with the same central sound speed, see the text).

4. The Radial Stability of the Self-Gravitating Spheres

4.1.Bonnor’s Stability Criterion

For a polytropic, self-gravitating sphere of specific heat ratio γ, Bonnor’s (1956, 1958) stability criterion can be written as:

where $r_e=R_e/R_c$ is the dimensionless outer radius of the sphere,

and

The derivation of equation (22) is completely analogous to the derivation of equation (17) of Raga et al. (2013) for the γ = 1 case.

We then take the numerically derived (p´(r) density stratifications described in § 3, calculate f(r), $\dot{f}\left(r\right)$ and F(r) (equations 23-24), and then evaluate the term on the left-hand-side of the inequality of equation (22) to find the smallest radius r_s at which the inequality is violated. In this way we derive the radius r_s = R_s/R_c of the maximal stable sphere for a given value of γ.

4.2.The Maximal Stable Solutions

After determining the maximum radius for stability r_s = R_s/R_c from Bonnor’s criterion (equations 22-24), we calculate the dimensionless mass m_s = m(r_s) and center-to-edge density ratio ${\rho }_0/{\rho }_s=1/\rho \text{'}\left(r_s\right)$ of the maximal stable solution. This exercise is repeated for several values of the specific heat ratio in the $\gamma =1\to 4/3$ range. For γ > 4/3 , Bonnor’s stability criterion is satisfied for all radii of the density stratifications (see also Bonnor 1958).

In Figure 4, we show the dimensionless mass ${\gamma }^{3/2}{m}_s=3M_s/\left(4\pi {\rho }_0{R}_{c,1}^3\right)$ (see equations [Mr]-[mr]), outer radius ${\gamma }^{1/2}{r}_s=R_s/R_{c,1}$, and center to edge density ratio p₀/p_s of the maximal stable solution as a function of γ. We see that:

• p₀/p_s grows as a function of γ, and reaches high values for γ ≈ 1.3, finally diverging at γ = 4/3,

• ${\gamma }^{1/2}{r}_s$ monotonically grows as a function of γ, and reaches a ${\gamma }^{1/2}{r}_s\approx 5.4$ for γ = 4/3 = 4/3,

• the dimensionless mass ${\gamma }^{3/2}{m}_s$ shows very small variations for γ = 4/3 in the $1\to 4/3$ range.

Fig. 4
Properties of the maximal stable solution as a function of γ. Bottom frame: center to edge density ratio p₀/p_s. Center: dimensionless outer radius ${\gamma }^{1/2}{r}_s=R_s/R_{c,1}$. Top: dimensionless mass ${\gamma }^{3/2}{m}_s$ (see equation 12).

For the $6/5<\gamma \le 4/3$ range (with density stratifications that end at an outer radius R_max, see § 2 and Figure 3), we obtain numerically the dimensionless outer radius $r_{max}=R_{max}/R_c$ and the mass m_max = m(r_max) (see equation 12) of the full stratification. In Figure 5, we plot the outer radius ${\gamma }^{1/2}{r}_{max}=R_{max}/R_{c,1}$ and the ratios $M_s/M_{max}$ and $R_s/R_{max}$ between the mass/radius of the maximal stable solution and the corresponding values of the full stratification.

From Figure 5we see that:

• the outer radius ${\gamma }^{1/2}{r}_{max}$ of the full stratification decreases monotonically from its infinite γ = 1.2 value to a value of ≈ 5.63 at γ = 4/3,

• the $R_s/R_{max}$ ratio increases from 0 (for γ = 1.2) to 1 (for γ = 4/3), with a faster increase as 𝛾=4/3 is approached,

• the $M_s/M_{max}$ ratio gradually increases from (for γ =1.2) to 1 (for γ = 4/3).

Fig. 5
Bottom frame: dimensionless maximum radius ${\gamma }^{1/2}{r}_{max}=R_{max}/R_{c,1}$ as a function of γ. Top frame: ratios M_s/M_max and R_S/R_max of the mass M_s and radius R_s of the maximal stable solution to the maximum mass M_max and radius R_max of the polytropic spheres as a function of γ. We show the results obtained for $\gamma =(6/5{)}^{+}\to 4/3$, which have density distributions which go to zero at a maximum radius R_max.

Therefore, we see that the fraction of the total radius and mass of the finite, γ > 1.2 self-gravitating spheres that is stable increases with γ. The (radial or mass) fraction that is stable tends to 1 as γ =4/3 is approached.

5. Summary

We have studied the non-singular density stratifications that are solutions to the polytropic Lane-Emden equation for self-gravitating spheres. We have first presented general analytical considerations about the proper adimensionalisation of the density and radius (using the second-order, small radius approximation to the non-singular density stratification), the derivation of the singular solutions, and the convergence at large radius of the non-singular to the singular solutions (which occurs only for γ < 6/5). The results that have been presented (see § 2) can be found in the previous literature (notably, in Chapter IV of the book of Chandrasekhar 1967 and in the book of Horedt 2004).

We have then computed numerical integrations of the Lane-Emden equation to obtain non-singular density stratifications for γ in the $1\to 4/3$ range (see § 3). Solutions around γ = 1.2 show the curious transition between the γ < 1.2 density stratifications (which at large radii converge to the singular solutions) and the γ < 1.2 stratification (which is analytic, and differs quite dramatically from the singular solution at large radii). It is not clear whether or not this rather curious mathematical feature has any relevance for astronomical applications.

Finally, we have applied Bonnor’s stability criterion (for self-gravitating solutions of arbitrary γ, see equation 22) to our numerically obtained density stratifications. We find that for spheres with $\gamma =1\to 4/3$ but with the same central temperature (or isothermal sound speed), the maximal stable sphere has (see Figure 4):

• a center to edge density ratio that grows quite substantially with increasing γ, diverging for $\gamma \to 4/3$,

• an outer radius that slowly grows with γ in the $1\to 1.2$ range, and has a steeper growth for $\gamma \to 4/3$,

• an almost constant mass as a function of γ.

Also, for the finite, γ > 6/5 spheres, the maximal stable sphere has a mass and radius that grow with γ. For γ > 4/3, the whole finite sphere is stable.

Applying these results to molecular clouds cores, we see that relatively small deviations from isothermality (which could be modelled with γ values somewhat larger than 1) are not likely to lead to the formation of substantially different self-gravitating structures, in the sense that the maximal stable solutions will have almost the same mass (regardless of the value of γ in the $1\to 4/3$ range). The outer radius of the maximal stable sphere grows for γ > 1, reaching a factor 1.5 larger value for γ = 1.2. The largest effect is found for the center to edge density ratio, which grows by a factor of ≈ 2 when γ is changed from 1 to 1.2 and by a factor of ≈10 when γ =1.3.

Therefore, we conclude that small deviations from isothermality will result in the production of maximal stable spheres with approximately the same mass, with somewhat larger outer radii, and with substantially higher center to edge density ratios. Clearly, it will be interesting to see whether the hydrodynamical collapse of the $\gamma =1\to 4/3$ density stratifications resembles the well studied collapse of the isothermal (γ = 1) sphere.

We end by discussing the work of Kandori et al. (2005), who use near-infrared images to derive the column density distribution for a set of ≈10 Bok globules. They carry out fits of the non-singular solution of the isothermal Lane-Emden equation to the derived structures. From these fits, they find the center to edge density ratios and the dimensionless outer radius ${\xi }_{max}$ (which, due to a different adimensionalisation is a factor of $\sqrt{6}$ larger than the ${\gamma }^{1/2}{r}_s$ outer radius shown in our Figure 4) of the observed structures. Most interestingly, they find that while some of the globules of their sample are close to the stability limit, about half of them have center to edge density ratios and outer radii larger than the stability limit of the isothermal sphere. This result leads them to conclude that these globules above the isothermal stability limit are probably already collapsing.

Our models provide a possible alternative explanation for these results. As one can see from Figure 4, self-gravitating spheres with departures from isothermality (modeled as spheres with γ > 1) have maximal stable solutions with center to edge density ratios and dimensionless outer radii larger than the ones of the maximal isothermal (γ =1) solution. The higher center to edge density contrast globules observed by Kandori et al. (2005) could therefore correspond to stable self-gravitating spheres with such departures from isothermality.

Acknowledgments

AR acknowledges support from the DGAPA-UNAM grant IG100218. JCRR acknowledges the Brazilian agency FAPESP grant 2017/12188-5. We thank an anonymous referee for pointing out the work of Kandori et al. (2005), which provides a most interesting connection between our models and reality.

References

Bonnor, W. B. 1956, MNRAS, 116, 351

______, 1958, MNRAS, 118, 523

Chandrasekhar, S. 1967, An introduction to the study of stellar structure (Dover, New York)

Falgarone, E., Puget, J. L. 1985, A&A, 142, 157

Horedt, G. P. 2004, Polytropes - applicaations in astrophysics and related fields (Kluwer: Dordrecht)

______, 2013, ApJ, 773, 131

Hunter, C. 2011, MNRAS, 328, 839

Ito, Y., Poje, A., Lancelotti, C. 2018, NewA, 58, 15

Kandori, R., Nakajima, Y., Tamura, M. et al. 2005, AJ, 130, 2166

Raga, A. C., Rodríguez-Ramírez, J. C., Rodríguez-González, A., Lora, V., & Esquivel, A. 2013, MNRAS, 49, 127

Shu, F. H. 1977, ApJ, 214, 488

Notes