Gravitational instabilities of isothermal spheres in the presence of a cosmological constant

Gravitational instabilities of isothermal spheres are studied in the presence of a positive or negative cosmological constant, in the Newtonian limit. In gravity, the statistical ensembles are not equivalent. We perform the analysis both in the microcanonical and the canonical ensembles, for which the corresponding instabilities are known as `gravothermal catastrophe' and `isothermal collapse', respectively. In the microcanonical ensemble, no equilibria can be found for radii larger than a critical value, which is increasing with increasing cosmological constant. In contrast, in the canonical ensemble, no equilibria can be found for radii smaller than a critical value, which is decreasing with increasing cosmological constant. For a positive cosmological constant, characteristic reentrant behaviour is observed.


Introduction
In a seminal work [1], Antonov described a thermodynamic instability of self-gravitating systems in the microcanonical ensemble, that later became known as 'gravothermal catastrophe' [2]. A classic review on thermodynamics and statistical mechanics of selfgravitating systems is the one of Padmanabhan [3] and a more recent one is written by Katz [4]. An extension of Antonov's instability to the canonical ensemble, named 'isothermal collapse', was given by Chavanis [5]. Extended reviews on the self-gravitating gas at thermal equilibrium, from the statistical mechanics point of view, are given by de Vega & Sanchez [6,7] and Destri & de Vega [8]. Thermodynamics of self-gravitating systems can be realized as the pioneering part of a, nowadays, more general, rapidly developing, new field of research, that is the thermodynamics of systems with long-range interactions [9,10,11].
In a recent letter [12] we reported on the effect of the cosmological constant to the Antonov's gravito-thermal instability in the microcanonical ensemble. In the present of gravitational systems. In this perspective, we study the simpler possible setup, i.e. a spherically symmetric, Newtonian, bounded system, to gain a basic understanding of the effect of the cosmological constant to the stability of self-gravitating systems.
It has been proved by de Vega and Siebert [28,29] that a thermodynamic limit, different from the usual one, does exist for a self-gravitating gas in the presence of a cosmological constant and that the mean field approximation correctly describes this limit. In the presence of the cosmological constant we find that, in the mean field approximation and the Newtonian limit, the negative cosmological constant (AdS case) tends to destabilize the system, while the positive cosmological constant (dS case) tends to stabilize it. This result further supports recent investigations of AdS instabilities [30,31]. In dS case many novel features arise. The system presents a reentrant behavior [12]. A second critical radius, above which metastable states are restored, emerges, and in the canonical ensemble the system undergoes reentrant phase transition, since there appear two critical temperatures. Reentrant phase transitions were known to occur for statistical systems with long-range interactions [32,33,34,35] but not for gravitating systems. In addition, the homogeneous solution of dS has a turning point of stability, which we calculate analytically, and there exist infinite non-uniform solutions in the homogeneous radius as well as multiple series of equilibria for any radius.
The paper is organized as follows. In section 2 we calculate the entropy extrema, in section 3 we present the criteria for stability that we used in our analysis, in section 4 the temperature and energy of the system are calculated and in section 5 the way to numerically generate the series of equilibria is presented and their asymptotic behavior is analytically studied. In section 6 is studied the stability of the homogeneous solution in dS in both ensembles. The main results of the microcanonical ensemble are presented in section 7, where the correspondence of our 'dS case' with the Schwartzschild-dS space is discussed, as well. The results in the canonical ensemble are given in section 8.

Entropy extremization and free energy
Consider N >> 1 identical particles (stars), bounded inside a spherical shell with insulating and perfectly reflecting walls. In order to calculate the entropy of the system, one should calculate the N -body distribution function f N ( r 1 , . . . , r N , p 1 , . . . , p N ). This seems an impossible task. However, if the correlations between the particles are not significant and an intermediate scale where the granularity of the system can be ignored exists, one can work in the mean field approximation [3] using the 1-body distribution function f ( r, p, t). This can be defined by the N -body distribution function as f ( r 1 , p 1 , t) = f N ( r 1 , . . . , r N , p 1 , . . . , p N )d 3 r 2 . . . d 3 r N d 3 p 2 . . . d 3 p N (1) or one can think as f giving the mass dm inside a volume d 3 rd 3 p: We assume that all particles have massm = 1, so that we work with the velocity υ instead of the momentum p. Once f has been determined, the density ρ( r) can be found by and the number of particles by The total mass is M = Nm. The question posed is which f extremizes the Boltzmann entropy with constant energy E and constant number of particles N , where d 6 τ = d 3 rd 3 υ. For the case without a cosmological constant, it has been proved [1,2] that only spherical configurations maximize the entropy. For a discussion on spherical configurations in the presence of a cosmological constant see Ref. [36]. We will consider only spherical, static distributions. In the derivation of f we follow [2]. Using the Lagrange's multipliers β = 1/kT , µ the variation condition with respect to f is The energy is E = K + U , where K is the kinetic energy Regarding the gravitational potential energy U , we have to work on the gravitational potential φ( r). In the Newtonian limit (see [37] for dynamical effects of the cosmological constant in the Newtonian limit), the Poisson equation in the presence of a cosmological constant Λ [38] is where ρ Λ = Λc 2 8πG . For an analytical derivation see Appendix Appendix B. The validity of the Newtonian approximation is analytically discussed in [12]. In the Newtonian limit it should hold ΛR 2 ≪ 1 while the cosmological constant is negligible if where R S = 2GM/c 2 is the Schwartzschild radius of the system. Since for Newtonian systems it normally is R/R S ≫ 1, the cosmological constant is not in principle negligible. For spherically symmetric configurations, bounded in r ∈ [0, R], the potential can be written as Therefore the potential energy can be written as Using equations (7), (12) the variation of energy is Using equations (4), (5) and (13) we get So that, in order for equation (6) to hold for all δf , we get for A = e µ−1 . We get the Maxwell-Boltzmann distribution likewise 'flat case'. The cosmological constant enters to the equation implicitly through the potential. Now, the density distribution can easily be calculated Absorbing the constants to the initial value φ(0) of the field and the central density ρ 0 , we finally get In spherical coordinates, the equation (8), substituting equation (17), gives This equation holds for r ≤ R, where R is the radius of the bounding wall. For r > R, it is of course 1 r 2 d dr r 2 d dr φ(r) = 0 and φ, φ ′ should be continuous at R. Let introduce the dimensionless variables Then, equation (18) becomes which we call the Emden-Λ equation. The initial conditions are The first initial condition comes from equation (19), while the second from spherical symmetry (gravitational field at the center is zero). We call z the value of x at R: We want to generate series of equilibria, that is to find y(z) for various z , λ (various isothermal spheres) and not just find y(x) for some λ. This should be done by solving (20) for various z, λ keeping N , i.e. M , constant at each case. For the 'flat' λ = 0 case this was very easy to perform. One could just solve Lane-Emden equation for one z and interpret the result y(x) as y(z) with no inconsistency, since the mass could be considered fixed: there is no mass scale for this system. Another difficulty that enters, is that varying λ cannot be realized as varying ρ Λ , since λ contains ρ 0 , as well, that is different for each equilibrium configuration. Therefore, the situation in general becomes rather complex. We will see in section 4, how we resolved the problems by introducing a new parameter and constructing an appropriate computer code. The canonical ensemble can be studied by defining the Helmholtz free energy F = E − T S with use of the Boltzmann entropy (5). It is equivalent to work with the Massieu function [40,5] It is shown by Chavanis [5] that the maximization of J with constant T is equivalent to the maximization of S with constant E to first order in variations δρ. This means that the two ensembles have the same equilibria defined by the distribution function It is easy to check, performing for J the previous calculations for S, that this holds true in the presence of a cosmological constant, too. What is different in the two ensembles is the stability analysis, i.e. the second order variation of entropy and free energy as we will see in section 3.

Criteria for stability
We want to calculate the second order variation of entropy w.r.t. perturbations δρ. We follow closely Padmanabhan [3]. Maximizing the entropy for a given distribution ρ(r) we get the Maxwell-Boltzmann distribution The entropy can therefore be written as We vary ρ(r) and therefore U , K through φ, T respectively keeping E and M fixed. Keeping E fixed gives the constraint δK + δU = 0 From this equation we derive in Appendix C that So that, we get with the second order variation being If at an equilibrium, for any perturbation δρ it is δ (2) S| equil < 0 then the entropy is a (local) maximum and the equilibrium is (locally) stable. If there exists one or more perturbations for which δ (2) S| equil > 0, then the equilibrium is unstable. Let us see how the sign of δ (2) S| equil can be deduced by an eigenvalue equation [3] generated by (27). Since the total mass is constant, for δρ should hold: Let us concentrate on spherical symmetric perturbation δρ = δρ(r) and introduce the mass perturbation q(r) = δM (r) Then The force due to perturbed distribution (δφ) ′ = G q r 2 has to be finite everywhere and therefore q should go like q → r 3 for r → 0. This means that q(0) = 0. Then equation (28) gives that q(R) = 0. Thus, the boundary conditions are Substituting (29), (30) into equation (27) and performing several integrations by part, we get The last expression can be written as The sign of δ (2) S is therefore determined by the eigenvalues of the 'matrix' K(r, r ′ ) At this equilibrium, where there is a transition from stability to instability, it should be ξ = 0. If for an equilibrium is found one perturbation F ξ for which ξ > 0 then this equilibrium is unstable. For an equilibrium to be (locally) stable all eigenvalues should be negative ξ < 0 for all perturbations. Equation (35) gives In (36) the cosmological constant enters implicitly, since: The boundary conditions of (36) are as given in (31) We developed an algorithm that can determine eigenvalues and eigenstates for the boundary value problem defined by equations (36), (37), (38). The main difficulty is that 8 in V enters the unknown function F ξ . We resolve the problem as follows. For a given range of ξ, the problem is solved for trial values of V , call them V T , and then the integral (37) is calculated, which gives some valueṼ . Some value ξ is indeed an eigenvalue, only ifṼ = V T and, in this case, of course V =Ṽ = V T . The algorithm is applied for the dimensionless version of (36), namely: and z, y, y ′ , B are calculated at the equilibrium. Using the algorithm we can determine a turning point, where ξ = 0, an instability ξ > 0 or verify a stable branch of series of equilibria by checking every equilibrium point for a zero eigenvalue and for an as large as possible range of positive eigenvalues. We performed these tasks for every series of equilibria demonstrated in this paper.
In an equilibrium for which ξ = 0, there is a transition from a stable branch to an unstable branch or from an unstable branch to a more unstable branch (or vice versa, of course). Suppose you approach the turning point from a stable series (all eigenvalues negative), then at the next equilibrium point after the transition, one eigenvalue becomes positive. In appendix Appendix D we show how one can determine the branch with the additional positive eigenvalue (instability) near a turning point, from the previous analysis.
Performing similar calculations for the second variation of free energy J, it is straightforward to find the corresponding eigenvalue problem for the canonical ensemble: Compared to eigenvalue equation (36), we see that the difference is only the absence of the term containing the derivative of the potential. This difference changes the onset of the instability for the two ensembles.
There is a way to study stability without solving an eigenvalue problem, due to a classical result of Poincaré [39]. In thermodynamics of self-gravitating systems, it was for the first time applied by Lynden-Bell & Wood [2]. Briefly, it states that a change of stability can only occur at a point, where two or more series of equilibria have a common point or where they merge into each other. As indicated by Katz [40], practically this means, that in the microcanonical ensemble, the change of stability happens in the point where β(E) has infinite slope, while in the canonical where β(E) has extrema. Equivalently this means that E in the microcanonical or β in the canonical ensemble, has an extremum with respect to some other variable (e.g. the density contrast log(ρ 0 /ρ R ) or z in dimensionless variables) at the turning point. We determine the point of change of stability in the dimensionless variables by finding an extremum of ER/GM 2 in the microcanonical case and GM β/R in the canonical case, with respect to the logarithm of the density contrast log ρ0 ρR .

Temperature and Energy
At the thermodynamic equilibrium, the gas sphere has the same temperature β everywhere, so that it is called an isothermal sphere. We define the dimensionless inverse temperature [2] which can be calculated with use of the dimensionless variables (19), by integrating the Poisson with Λ equation (18): The kinetic energy per particle is To calculate the total energy, we will use the Virial theorem, so that to avoid performing one more numerical integration to calculate the Newtonian potential energy and therefore improving the computer's performance. However, in the Appendix Appendix E we perform a straightforward calculation of the energy in order to numerically cross-check the two expressions for some cases, who are proven to give identical results. This confirms the fact that bounded isothermal spheres are virialized in the presence of the cosmological constant. The Virial theorem for a discrete collection of matter: can be generalized in our case as where V = 4 3 πR 3 is the volume of the shell and P the pressure exerted on the gas by the walls. The term 3P V arises simply since F · R = P · 4πR 2 · R = 3P V . The contribution of the Newtonian potential to the right-hand side is just the Newtonian potential energy U N , while for the cosmological potential φ Λ we get the term Substituting everything in equation (44) we get the following form of the Virial theorem: Using (19) we write the potential U Λ in the dimensionless variables Let us calculate the term 3P V in dimensionless variables We define the dimensionless energy We can calculate Q from the virial equation (46) RE For simplicity let us focus in the microcanonical ensemble for the moment. The selfgravitating gas is characterized by two instabilities [3]; a 'strong' instability that is associated with complete absent of equilibria and a 'weak' instability which refers to equilibria that are unstable, i.e. entropy, although it is an extremum, is not a local maximum. The equilibrium point at which the weak instability sets in, called the turning point, is therefore the point at which the second variation of entropy becomes zero, as calculated in the previous section. However, Poincaré's theorem insures [39,40] that this point is the same with the marginal point of the strong instability. In Figure 1, where the series of equilibria β = β(E) is drawn for the flat Λ = 0 case, this is point B. The strong instability corresponds to the region at the left of the vertical dashed line that crosses B. The weak instability refers to the branch BS, where S is the focal point of the spiral. Thus, point B is simultaneously the turning point of stability (from stable branch AB to the unstable branch BS) and the marginal point of the strong instability. This strong instability leads to a core-halo structure as verified by monte-carlo simulations [6,15,16,17,18] and is associated with a collapse phase transition [6,7,8,29], while the weak instability leads to a fractal structure and is associated with a fragmented collapse, called a clumping phase transition [5,6,7,8,29]. This fractal structure is due to the secondary instabilities that set in at points S 1 , S 2 , etc. All of the above hold in the canonical ensemble, as well, where the axes in Figure 1 should be interchanged and entropy should be replaced with free energy. 11 Figure 1: The series of equilibria β = β(E) for Λ = 0. In the microcanonical ensemble, AB is the stable branch and B is the turning point of stability. As we can see, it is also the marginal point for which equilibrium states do exist.

The series of equilibria and asymptotic behavior
We want to solve the Emden-Λ equation (20) for various ρ Λ keeping M constant and for various isothermal spheres with radius R. The cosmological constant introduces a mass scale We define the dimensionless mass is the mean density of matter. This gives Equation (52) implies that in order to keep m fixed, λ has to be different at each z.
We developed an algorithm that solves the Emden-Λ equation keeping the quantity m constant; for each z, some λ values are iterated until the requested value of m, calculated 12 by equation (52), is found for a relative tolerance predetermined by the user. For the various plots in this paper, we used relative tolerance between 10 −7 − 10 −11 depending on the needs of each case. From equation (51), it is evident that solving for various fixed m can be interpreted as varying ρ Λ and/or R for a fixed M . Therefore we can determine how various quantities change with respect to ρ Λ by solving for various m.
For the AdS case (ρ Λ < 0) we find that for each ρ Λ only one series of equilibria exists, likewise flat (ρ Λ = 0) case. However, for dS (ρ Λ > 0), we find that for a fixed ρ Λ , there exist more than one series of equilibria. This is evident in Figure 2 where m is plotted w.r.t. λ, z. We see that the intersection of a plane m = const with the m-surface defines various different curves in the (λ, z) space.
The Emden equation (20) for ρ λ = 0, i.e. λ = 0, is well known to have an exact but singular solution with infinite density at the origin since e −ys = 2 z 2 . It can be shown [2] that for z → ∞ the series of equilibria approach the singular solution. Therefore, the singular solution in the flat case is the focus point of the central spiral β(E) in Figure 3. We can see in this picture that in dS and AdS cases there seems to exist an equivalent to the flat singular solution, which can be identified with the focal points of the spirals (note that in dS, not all series of equilibria do form a spiral β(E)). Unfortunately, no analytical solution of equation (20) is known. However, following the flat case paradigm we can determine asymptotically the equivalent of the singular solution in dS and AdS. Applying the Emden-Λ equation becomes This is the equation of an oscillator in the potential V (u) = e u − 2u with external force F ext = λe 2ζ − u ′ where prime denotes differentiation w.r.t. ζ. The external force is damping in AdS case (λ < 0) and forced, damping in dS (λ > 0). The potential has a minimum u 0 = log 2. In dS case, if the external force is damping, i.e. the term −u ′ dominates, then for ζ → ∞, we have u → u 0 . Therefore, in this case, we can make an expansion of u about u 0 . The damping dominates if: where y ′ denotes differentiation w.r.t. z and we have used equation (52). This limit means that we are considering very small ρ Λ , just about the flat case, since λz 2 → 0 has to be taken together with z → ∞. The two limits are consistent with each other, since by equation (52) we see that as z → ∞ it must be λ → 0 so as λz 2 to remain finite in order for m to be finite (and B is bounded for the equilibria we are considering).
Making the transformation u = log2 + u 1 14 equation (55) becomes Provided condition (56) holds and for ζ → ∞, u 1 is small and we can expand the exponential keeping the first two terms, to get: The solution of the corresponding homogeneous equation is well known to be and one solution of the full non-homogeneous equation can be found easily to be u p = λ 8 e 2ζ . We get in total the asymptotic behavior of the Emden-Λ equation for z → ∞ and for small λz 2 is The density ρ is given by the exponential e −y . We get the asymptotic behavior We see that the equivalence to the singular solution (53) of the flat case, is given in dS and AdS by the 'asymptotic singular' solution e yAS = 2 z 2 e λ 8 z 2 that is This solution y AS , is easy to check that indeed satisfies the Emden-Λ equation (20) to first order in λz 2 . By equations (42) and (49) we can calculate the temperature and energy of the asymptotic singular solution to be: Recall that m =ρ/2ρ Λ is used to vary the cosmological constant. So that m → ∞ corresponds to the flat case, m < 0 to AdS and m > 0 to dS. We find that the singular point determined analytically by equation (60)

Homogeneous solution in dS
For a positive cosmological constant, the Emden-Λ equation has solution with a uniform density (for r < R), we call homogeneous solution. For ρ = const. the Poisson equation gives φ N = 2πG 3 ρr 2 , so that the full potential with the cosmological constant is The Poisson with Λ equation (8)  Since φ ′ = 0 the homogeneous solution has the same turning points in the two ensembles, because the eigenvalue equations (36) and (40) are the same in this case. Therefore, although the following analysis is done in terms of the microcanonical ensemble, the results hold for the canonical ensemble, as well. So that all turning points for the various solutions presented in this section are identical in the two ensembles. The radius R H of the homogeneous solution is independent of energy and temperature and is given for a fixed mass M and cosmological constant ρ Λ by equation: The homogeneous solution resembles the Einstein's static universe ( [41,42] has found the Einstein's static universe to be a local entropy maximum among other possible universes) in the Newtonian limit. However, to be more precise there are many homogeneous configurations (uniform density) with different temperatures, from which only one is literally static, that is the solution with T = 0 (β → ∞). The question is whether all, some or none of these solutions are stable. The static solution is surely thermodynamically unstable since it corresponds to only one microstate, therefore it has the minimum possible 17 entropy, i.e. zero entropy. On the other hand the solution β = 0, which is allowed from (63), is stable since it behaves as an ordinary gas. More formally, equation (32) can be written as which gives δ (2) S < 0 for all perturbations q, for β = 0 (T → ∞). In dimensionless variables equation (66) reads which gives δ (2) S < 0 for B = 0, as well. Therefore, there exists a point of change of stability somewhere between B = 0 and B → ∞. From equation (64) we see that the energy has not an extremum. However, Poincaré's criterion does not exclude the possibility of having a change of stability at a point other than an energy extremum [40]. This is our case. The differential equation (39) for the homogeneous solution y = y ′ = 0 (and m = λ = 1) and for ξ = 0 becomes We want the smallest z that satisfies the second boundary condition F (z) = 0, that is the change of stability happens at this z that is a solution of the equation The solution, call it z 0 , can be found graphically (see Figure 4) to obtain z 0 ≃ 4.4934 which corresponds to B 0 ≃ 6.73 from equation (63). Therefore, the homogeneous solution is stable for temperature T > T 0 and unstable for T < T 0 , with As we can see in Figure 2 there are infinite series of solutions for the homogeneous radius R H that corresponds to m = 1 (the homogeneous solution corresponds to m = λ = 1). In Figure 5 is drawn the dimensionless energy Q = RH E GM 2 versus the density contrast log(ρ 0 /ρ R ) for two of these series. We see that there are solutions with ρ 0 < ρ R that continuously turn to solutions with ρ 0 > ρ R . At the point ρ 0 = ρ R , which corresponds to a homogeneous solution, there occurs a change of stability for the two solutions. That is because the corresponding energies Q 0 = −0.6771, Q 1 = −0.8246 correspond to the two first zero eigenvalues z 0 = 4.4934, z 1 = 7.7251 of the homogeneous 18 log(ρ0/ρR) solution as can easily be verified by equation (64). We numerically determined the unstable branch to be the one with ρ 0 > ρ R . One is forced to conjecture that this pattern of different series, for R = R H , is infinitely continued (as indicated by Figure 2) for lower and lower energies. The change of stability at ρ 0 = ρ R should always correspond to Q > −0.9 as indicated by equation (64).

Microcanonical ensemble
In this section, we review our results of Ref. [12]. Let investigate the solutions for R < R H and R > R H in the microcanonical ensemble. The energy versus the density contrast is plotted in Figure 6. The upper series labelled 1 (R < R H ) and 4 (R > R H ) correspond to monotonically changing density. For series 1 it is decreasing (ρ 0 > ρ R ), while for 4 increasing (ρ 0 < ρ R ). Series 1 suffers a change of stability at point B 1 (stable branch is A 1 B 1 ), while series 4 is stable everywhere and does not suffer any change of stability. This is proved as follows: for this series the limit E → ∞ does exist, which corresponds to β = 0. By equation (66) for β = 0 we get δ (2) S < 0. In addition, the energy does not have an extremum (where a transition to instability could occur) and (to be sure) the whole series is numerically checked at each point for a zero eigenvalue. No one is found. Every such solution (series 4) corresponds to configurations somewhat hollow at the center with matter concentrated mainly at the edge. The next series at more negative energies have one density extremum and at the next, one more is added and so on. At points B, except B 4 an instability sets in. For series A 2 B 2 and A 6 B 6 we have strong numerical evidence that are stable, while series 3 and 5 are found to be unstable. Series like A 2 B 2 and A 6 B 6 correspond to solutions diluted at the center with periodic condensations away from the center. One would normally expect this pattern in Figure 6 to continue as one finds series with more and more negative energy.

Critical quantities
A general result is that as the cosmological constant increases, isothermal spheres exist at even lower temperatures and energies. This can be seen in Figure 3, where the classic spiral β(E) is drawn for some positive and some negative cosmological constant. One could say that AdS destabilizes, while dS stabilizes the system.
In the flat case, isothermal spheres exist only for This means that for a fixed radius, there exists a minimum critical energy E cr = −0.335GM 2 /R down to which, equilibria do exist. We want to determine how this critical energy changes with respect to ρ Λ . The answer lies in Figure 7. The critical energy decreases for increasing cosmological constant. In AdS the critical energy becomes positive for ρ Λ −4.2ρ, whereρ is the mean density. Assuming the energy E and radius R are fixed at values that respect (69) in flat case, does not guarantee that the equilibrium is stable. This depends on the density contrast ρ 0 /ρ R , that is the ratio of the central density versus the edge density. The situation is similar in the presence of ρ Λ . For ρ Λ = 0 it is (ρ 0 /ρ R ) cr = 709 [1] with the unstable branch being the one with ρ 0 /ρ R > 709. In Figure 8 we see how this number changes w.r.t. the cosmological constant. The critical density contrast decreases for increasing cosmological constant. The instability occurs in AdS at more condensed configurations, while in dS at less condensed configurations.
Assume that the energy is negative and fixed at some value E in the flat case. Then, inequality (69) shows that there is a maximum critical radius R A = (−0.335/E)GM 2 that constrains the existence of an equilibrium. For R > R A there are no equilibria. In  Figure 8: The critical density contrast versus ρ Λ for fixed M , R in the microcanonical ensemble, whereρ is the mean density of matter. An instability sets in for ρ 0 /ρ R > (ρ 0 /ρ R )cr (except when (ρ 0 /ρ R )cr < 1) where a clumping phase transition occurs.
AdS this radius decreases as ρ Λ attains more negative values. In dS this radius increases as the cosmological constant increases and in addition there appears a second critical radius, we call R IA , that constrains the existence of equilibria from below. At R A and R IA a collapse phase transition takes place. That is, no equilibria exist for R A < R < R IA as can be seen in Figure 9 and the system lies in a collapsed phase. This is a typical reentrant behavior, that is common to statistical systems [32,33,34,35], whenever competing interactions are present. Beyond the marginal value ρ dS Λ ≃ 7.14(3|E| 3 /8πG 3 M 5 ) there can always be found equilibrium states in dS case.
In region I of Figure 9 there exist series 1 of equilibria of Figure 6(a) and in region II all equilibria of Figure 6(b). In the small upper gray shaded region there are the rest series of Figure 6(a). These type of equilibria exist only for values of the cosmological constant greater than a minimum value ρ min Λ . This is the smallest value for which the cosmological force can keep marginally the whole matter at the edge. It can easily be calculated by equating the forces at the edge, assuming all matter is concentrated at R: whereρ is the mean density.

Comparison with Schwartzschild-dS space
The two critical radii in Figure 9 resemble the two horizons of Schwartzschild-dS space, where the role of the cosmological horizon is played by R IA . The Schwartzschild- dS metric can be written as: where ρ Λ is the 'mass' density of the cosmological constant (the energy density is κ = ρ Λ c 2 = Λc 4 8πG ). This metric has two horizons for ρ Λ > 0; the black hole horizon R BH and the cosmological horizon R C . Both are defined as the real roots R H of the third order polynomial equation which, for various values of the cosmological constant ρ Λ , are plotted in Figure 10. The resemblance with Figure 9 is too much striking to be considered accidental! It seems as though the reentrant phenomenon of Figure 9 is the closest Newtonian analogue of the horizons ( Figure 10) in Schwartzschild-dS space. However, there is a big deference between the Schwartzschild-dS space and the Newtonian reentrant phenomenon. It is the opposite sense of the inequality for the instability, i.e. the stable region in the Newtonian case corresponds to the unstable region (R < R BH and R > R C ) of the Schwartzschild-dS space.   Figure 12: The critical density contrast versus ρ Λ for fixed M , R in the canonical ensemble, whereρ is the mean density of matter. An instability sets in for ρ 0 /ρ R > (ρ 0 /ρ R )cr (except when (ρ 0 /ρ R )cr < 1), where a clumping phase transition occurs.
In the canonical ensemble [13] the instability sets in at different points than in the microcanonical ensemble. The series of equilibria expressed as function of temperature with respect to the density contrast can be seen in Figure 11 in case of dS. In AdS, just like flat case, there exist only one series of equilibria for a fixed cosmological constant, that is similar in form to series 1 of Figure 11(a). We note that Figure 11 of temperature is very similar to Figure 6 of the energy. At points B, except B 4 , of Figure 11, an instability sets in. Series A 1 B 1 and A 4 B 4 are stable. A safe conclusion on the stability of the green series A 2 B 2 , A 3 B 3 , A 5 B 5 and A 6 B 6 could not be reached, although the positive specific heat indicates stability, since we work in the canonical ensemble. The rest series are unstable.
In Figure 12 we can see that the critical density contrast for series 1 of Figure 11(a), is decreasing for increasing cosmological constant, just like in the microcanonical ensemble. However, the instability sets in earlier in the canonical ensemble. It is well known in the flat case, and this fact is true with a cosmological constant, too, that the instability in the canonical ensemble sets in when the specific heat becomes negative. This negative specific heat region is stable in the microcanonical ensemble, while the instability sets in when the negative specific heat becomes positive again, in this ensemble (see [2]).
The canonical ensemble is completely different then the microcanonical ensemble as far as the stability of the system is concerned. We alway consider the mass to be fixed in both ensembles. As can be seen in Figure 13, the critical radius is decreasing with increasing cosmological constant and the region of the instability is now for radii smaller (and not bigger) then the critical radius R < R cr . The critical radius is decreasing for increasing ρ Λ , because for an increase in ρ Λ the temperature is decreased, so that one should compress the system to balance out this destabilizing temperature decrease. The region of the instability changes, because in a compression although the pressure gradient is increased, it is weakened compared to the microcanonical ensemble, due to the heat transferred to the heat bath. In AdS case there is a marginal value of the cosmological The critical temperature is decreasing with increasing cosmological constant, as can be seen in Figure 14. Increase of the cosmological constant acts as a stabilizer on the system due to the increase of the repelling force (or the decrease of the attracting force in case of AdS), enabling the system to be stable at lower temperatures. As the cosmological constant increases, it reaches a value for which the system can marginally be in static dynamical equilibrium. At this state all matter is still, i.e. T = 0, and is concentrated at the edge. This is point A in Figure 14. We have calculated this point earlier in equation (70) and found ρ min Λ =ρ/4. For greater values the outward pointing cosmological force is increasing, enabling the mass to approach towards the center and to greater temperatures. The system undergoes a reentrant phase transition in the canonical ensemble. For a fixed cosmological constant at this region (after point A), there are metastable states for low temperatures up to some maximum critical value T 1 , where a collapse phase transition occurs. For greater temperature, there exist no metastable states and the system suffers isothermal collapse. This happens up to some second critical temperature T 2 . For even greater temperatures, the equilibria are restored.

Conclusions
In the presence of the cosmological constant we find that, in the microcanonical ensemble, the critical radius, namely the Antonov radius, above which the system becomes unstable is increasing with increasing (negative or positive) cosmological constant, while in the canonical ensemble it is decreasing. The critical energy and temperature (less than which the system becomes unstable), as well as the critical density contrast (ρ 0 /ρ R ) cr , are decreasing in the two ensembles.
In dS case a new phenomenon is discovered, namely a reentrant phase transition; in microcanonical ensemble there emerges a second critical radius above which metastable states are restored, while in the canonical there appears a second critical low temperature, lower than which equilibria are restored. We stress out the similarity of the behavior, we have discovered, of the two critical radii in dS case in the microcanonical ensemble with the two horizons of relativistic Sschwartzscild-dS space. It seems that our 'dS case' is a Newtonian analogue of Schwartzschild-dS system.
Another interesting feature of dS is the turning point of stability for the homogeneous solution, which resembles the Einstein's static universe, and the infinite numbers of nonuniform solutions at the homogeneous radius which suffer a transition from stability to instability when passing from solutions with density contrast ρ 0 /ρ R < 1 to the ones with ρ 0 /ρ R > 1. In addition, there exist multiple series of equilibria for a given positive cosmological constant and fixed radius and mass.
We stress out, that in this work the non-equivalence of ensembles in gravitating systems is confirmed in the most dramatic way. For a fixed mass, in the microcanonical ensemble the instability occurs for radii larger than a critical value, while in the canonical the instability occurs for radii smaller than a critical value. That is because in the canonical ensemble, the pressure gradient that balances gravitation, is not drastically increased during a compression of the system, due to the heat transfer to the heat bath. In contrast, in the microcanonical ensemble where there is no energy loss, the pressure gradient is increased during a compression, drastically enough to balance gravity.
Regarding applications to the physical world (positive cosmological constant), we state two interesting issues raised by our work. In the microcanonical ensemble, our analysis could have implications to the evolution of galaxy clusters, since many of them present a core in their centre and others a supermassive black hole. A quick calculation shows that, for the observed value of ρ Λ , the relevant to our stability analysis, dimensionless quantities 2ρ Λ /ρ and 8πG 3 M 5 ρ Λ /3|E| 3 are of order unity for some typical values of regular galaxy clusters [12]. This implies that the cosmological constant could influence the onset of the instability of galaxy clusters. In the canonical ensemble, the cosmological constant could have an effect on the fractal structure of the Universe, since it is connected with the secondary instabilities [5]. Our analysis on the asymptotic behavior of the equilibria (section 5), can be used to calculate the secondary instabilities.
Let us close, noting that, even though there might be phenomenological connection with the physical world, our work is mainly focused and wishes to contribute on the theoretical understanding of the impact of an arbitrary cosmological constant term to the stability of gravitational systems.
Let 0 < α < 1 and αM be the mass uniformly distributed inside a sphere with radius r 1 and (1 − α)M the mass uniformly distributed inside a second sphere with radius r 2 . The distance of the spheres' centers is r 12 . The distribution functions are constants where υ 1 , υ 2 the corresponding to the two spheres velocity bounds and the velocities υ are equally probable by construction. Let us calculate the entropy of the system: The Newtonian dynamical energy U 1 of the first sphere is and similarly for the second sphere Let our coordinate system be centered at the center of the first sphere and the center of the second sphere be at z = r 12 . The cosmological dynamical energy for the first sphere is For the second sphere we use the cylindrical coordinates (τ, φ, z) (see Figure A.15). Its cosmological dynamical energy is The Newtonian dynamical energy between the two spheres is The kinetic energy of the first sphere is  gives S → ∞. Therefore there exists a configuration with finite energy for which entropy is not bounded from above, or equivalently, matter can always be redistributed in such a way keeping energy fixed and increasing the entropy.

Appendix B. Poisson equation with Λ
Let us calculate the equation of the gravitational potential φ in the presence of a cosmological constant Λ, in the Newtonian limit. The Einstein's equations are We assume ρ Λ to be fixed. Using (C. In the third raw we used the identity 3) The constraint δE = 0 gives by use of (C.2): Appendix D. How to determine the unstable branch near a turning point We search for a solution of the problem (36) for ξ = 0. LetT be the operator Let V T be the trial value we use to solve (36). If it does not correspond to the solution F 0 , it would correspond to some other solution F n for a different eigenvaluê F 0 will correspond to a solution with a zero eigenvalue only if

Near the turning point it is
R 0 dr F 0 F n ≃ R 0 dr F 2 0 so that an instability (ξ n > 0) sets in whenṼ < V T .