© ESO, 2014

1. Introduction

It was shown by Chernin (2001, 2008) that outer parts of galaxy clusters may be under strong influence of dark energy (DE), which was discovered by observations of supernovae (SN) Ia at redshift z ≤ 1 (Riess et al. 1998; Perlmutter et al. 1999) and in the spectrum of fluctuations of cosmic microwave background (CMB) radiation (see e.g. Spergel et al. 2003; Tegmark et al. 2004). To investigate these effects on the gravitational equilibrium of the clusters, solutions for polytropic configurations in presence of DE have been obtained by Balaguera-Antolínez et al. (2006, 2007) and Merafina et al. (2012). We here derive a Boltzmann-Vlasov kinetic equation in presence of DE and gravity in a Newtonian regime. The solutions generalize those obtained by Bisnovatyi-Kogan et al. (1993, 1998) for the kinetic equation without DE. The Newtonian approximation is chosen because the observed chaotic velocities of galaxies inside a cluster are much lower than the velocity of light.

The general relativistic solution in presence of DE can be applied for equilibrium configurations of point masses of some exotic particles that only interact gravitationally. In early stages of the Universe expansion, before and during the inflation stage, these particles may form gravitationally bound configurations that collapse during the inflation, when anti-gravity decreases. As a result of this collapse, such hypothetical objects may be transformed into primordial black holes that appear after the end of inflation. The relativistic kinetic equation and its solutions in presence of DE will be considered elsewhere.

2. Newtonian approximation in description of DE

The substance that is called now DE was first introduced by Einstein (1918) for a stationary universe in the form of the cosmological constant Λ during his unsuccessful attempts to construct a solution for a stationary universe. Shortly before this, de Sitter (1917) had shown that in presence of Λ the solution for an empty space describes an exponential expansion. Friedmann (1922, 1924) was the first to obtain exact solutions for the expanding universe that contained matter in presence of the cosmological constant Λ. Another exact solution for the metric in presence of Λ, around the gravitating point mass, was obtained by Carter (1973). This solution is a direct generalization of the Schwarzschild solution for a black hole (BH) in vacuum with a metric $\begin{array}{ccc}\mathrm{d}{\mathit{s}}^{\mathrm{2}}\mathrm{=}{\mathit{g}}_{\mathrm{00}}{\mathit{c}}^{\mathrm{2}}\mathrm{d}{\mathit{t}}^{\mathrm{2}}\mathrm{-}{\mathit{g}}_{\mathrm{11}}\mathrm{d}{\mathit{r}}^{\mathrm{2}}\mathrm{-}{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathrm{d}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{d}{\mathit{\varphi }}^{\mathrm{2}}\mathrm{)}\mathit{,}& & \end{array}$(1)of the form $\begin{array}{ccc}{\mathit{g}}_{\mathrm{00}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{g}}_{\mathrm{11}}}\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathit{GM}}{{\mathit{c}}^{\mathrm{2}}\mathit{r}}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{r}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{2}\mathit{GM}}{{\mathit{c}}^{\mathrm{2}}\mathit{r}}\mathrm{-}\frac{\mathrm{8}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}{\mathit{r}}^{\mathrm{2}}}{\mathrm{3}{\mathit{c}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(2)where the density of DE ρ_Λ is connected with Λ as $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{=}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi G}}\mathrm{·}& & \end{array}$(3)A transition to the Newtonian limit, where DE is described by the anti-gravity force in vacuum, was made by Chernin (2008). In the limit of a weak gravity (v² ≪ c², GM/r ≪ c²) the metric coefficients are connected with a gravitational potential Φ_g as (Landau & Lifshitz 1962) $\begin{array}{ccc}{\mathit{g}}_{\mathrm{00}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{1}\mathrm{+}\frac{{\mathrm{\Phi }}_{\mathrm{g}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(4)Then, using the Eqs. (4) and (2) at Λ = 0, we obtain the expression for the Newtonian potential Φ_g and the Newtonian gravity force acting on the unit mass F_g$\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{\mathit{r}}\mathit{,} {\mathit{F}}_{\mathrm{g}}\mathrm{=}\mathrm{-}\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{g}}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(5)For the Schwarzschild-de Sitter metric (2) we have in the Newtonian limit $\begin{array}{ccc}\mathrm{\Phi }\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{\mathit{r}}\mathrm{-}\frac{\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}{\mathit{r}}^{\mathrm{2}}}{\mathrm{3}}\mathit{,F}\mathrm{=}{\mathit{F}}_{\mathit{g}}\mathrm{+}{\mathit{F}}_{\mathrm{\Lambda }}\mathrm{=}\mathrm{-}\frac{\mathit{GM}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{8}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathit{r}}{\mathrm{3}}\mathrm{·}& & \end{array}$(6)In this way, the cosmological constant creates a repulsive (anti-gravity) force between a BH and a test particle in vacuum, which force increases linearly with a distance between them. The normalization of the potential here is chosen so that Φ_g = 0 at r = ∞, and Φ_Λ = 0 at r = 0.

We now consider the equilibrium of a self-gravitating object in presence of DE. In general relativity the equations describing the equilibrium in a spherically symmetric configuration in vacuum (without DE) have been derived by Oppenheimer & Volkoff (1939)$\begin{array}{ccc}& & \\ & & \frac{\mathrm{d}{\mathit{M}}_{\mathit{r}}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{4}\mathit{\pi \rho }{\mathit{r}}^{\mathrm{2}}\mathrm{·}\end{array}$(7)

Here ρ, P are the total density and total pressure of the matter, and M_r is the total (gravitating) mass, including a gravitationally binding energy, inside a radius r in the Schwarzschild-like metric $\begin{array}{ccc}& & \\ & & \\ & & \end{array}$

In presence of DE, ρ and P are represented as $\begin{array}{ccc}\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{m}}\mathrm{+}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{=}{\mathit{\rho }}_{\mathrm{m}}\mathrm{+}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi G}}\mathit{,} \mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{m}}\mathrm{+}{\mathit{P}}_{\mathrm{\Lambda }}\mathrm{=}{\mathit{P}}_{\mathrm{m}}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{4}}}{\mathrm{8}\mathit{\pi G}}\mathrm{·}& & \end{array}$(11)We consider a Newtonian limit when $\begin{array}{ccc}{\mathit{P}}_{\mathrm{m}}\mathrm{\ll }{\mathit{\rho }}_{\mathrm{m}}\hspace{0.17em}{\mathit{c}}^{\mathrm{2}}\mathit{,} \mathit{r}\mathrm{\gg }\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(12)Here we used a definition ${\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}\mathrm{=}\mathrm{4}\mathit{\pi }{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{r}}{\mathit{\rho }}_{\mathrm{m}}{\mathit{r}}^{\mathrm{2}}$dr. In the Newtonian limit we have from Eq. (7) $\begin{array}{ccc}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{{\mathit{\rho }}_{\mathrm{m}}\mathrm{(}\mathrm{3}\mathit{G}{\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}\mathrm{-}\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}}\mathrm{)}}{{\mathit{r}}^{\mathrm{2}}\left(\mathrm{3}\mathrm{-}\mathrm{\Lambda }{\mathit{r}}^{\mathrm{2}}\right)}\mathrm{·}& & \end{array}$(13)We estimate the last term in the denominator. For an equilibrium configuration with a finite radius to exist, we need a positive sign of the numerator; from condition (12) we have $\begin{array}{ccc}\mathrm{\Lambda }{\mathit{r}}^{\mathrm{2}}\mathit{<}\frac{\mathrm{3}\mathit{G}{\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{\mathit{r}{\mathit{c}}^{\mathrm{2}}}\mathrm{\ll }\mathrm{1}\mathrm{·}& & \end{array}$Therefore, in the denominator we have Λr² ≪ 3 and we can neglet the term with Λ. In the Newtonian approximation, in presence of DE, we obtain the following equilibrium equation $\begin{array}{ccc}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}{\mathit{\rho }}_{\mathrm{m}}\left(\frac{\mathit{G}{\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}\mathit{r}}{\mathrm{3}}\right)\mathrm{=}\mathrm{-}{\mathit{\rho }}_{\mathrm{m}}\left(\frac{\mathit{G}{\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{8}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathit{r}}{\mathrm{3}}\right)\mathit{,}& & \end{array}$(14)with ρ_Λ given by definition (3), which was used without derivation by Merafina et al. (2012). On the other hand, we can write the Poisson equation for the gravity of the matter together with the hydrostatic equilibrium equation $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{m}}\mathit{,} \frac{\mathrm{\nabla }\mathit{P}}{{\mathit{\rho }}_{\mathrm{m}}}\mathrm{=}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{\Lambda }}\mathit{,}& & \end{array}$(15)and then, the potential created by DE in the vacuum, taking into account that P = P_m + P_Λ, satisfies the Poisson equation $\begin{array}{ccc}{\mathrm{\nabla }}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{\Lambda }}\mathrm{=}\mathrm{-}\mathrm{8}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathit{,} {\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{=}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi G}}\mathrm{·}& & \end{array}$(16)This equation, together with the Poisson equation for the gravity of the matter fully describes a static gaseous equilibrium configuration in presence of DE. Similarly, we can write the hydrodynamic Euler equation in presence of DE as $\begin{array}{ccc}\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{+}\frac{\mathrm{\nabla }\mathit{P}}{{\mathit{\rho }}_{\mathrm{m}}}\mathrm{=}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{-}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{\Lambda }}\mathrm{·}& & \end{array}$(17)

3. Kinetic equation for a self-gravitating cluster in presence of DE

The kinetic Boltzmann-Vlasov equation for a distribution function f of non-collisional gravitating points of equal mass m in spherical coordinates (r,θ,ϕ) is written as $\begin{array}{ccc}\frac{\mathit{\partial f}}{\mathit{\partial t}}\mathrm{+}{\mathit{v}}_{\mathit{r}}\frac{\mathit{\partial f}}{\mathit{\partial r}}\mathrm{+}\frac{{\mathit{v}}_{\mathit{\theta }}}{\mathit{r}}\frac{\mathit{\partial f}}{\mathit{\partial \theta }}\mathrm{+}\frac{{\mathit{v}}_{\mathit{\varphi }}}{\mathit{r}\sin \mathit{\theta }}\frac{\mathit{\partial f}}{\mathit{\partial \varphi }}& & \end{array}$$\begin{array}{ccc}\mathrm{+}\left(\frac{{\mathit{v}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}}{\mathit{r}}\mathrm{-}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}\right)\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{v}}_{\mathit{r}}}\hspace{0.17em}\mathrm{+}\left(\mathrm{-}\frac{{\mathit{v}}_{\mathit{r}}{\mathit{v}}_{\mathit{\theta }}}{\mathit{r}}\mathrm{+}\frac{\cot \mathit{\theta }\hspace{0.17em}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial \theta }}\right)\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{v}}_{\mathit{\theta }}}& & \end{array}$$\begin{array}{ccc}\mathrm{+}\left(\mathrm{-}\frac{{\mathit{v}}_{\mathit{r}}{\mathit{v}}_{\mathit{\varphi }}}{\mathit{r}}\mathrm{-}\frac{\cot \mathit{\theta }\hspace{0.17em}{\mathit{v}}_{\mathit{\varphi }}{\mathit{v}}_{\mathit{\theta }}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\mathit{r}\sin \mathit{\theta }}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial \varphi }}\right)\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{v}}_{\mathit{\varphi }}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(18)where, in presence of DE, we have Φ = Φ_g + Φ_Λ. In a spherically symmetric stationary cluster, we have ∂Φ /∂t = 0 and Φ = Φ(r). Moreover, the kinetic Eq. (18) has four first integrals, written in Cartesian coordinates (x,y,z) as $\begin{array}{ccc}\frac{\mathit{E}}{\mathit{m}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{\Phi }\mathit{,}\\ \frac{{\mathit{L}}_{\mathit{x}}}{\mathit{m}}& \mathrm{=}& \end{array}$(19)In spherical coordinates, these integrals can be expressed by $\begin{array}{ccc}\frac{\mathit{E}}{\mathit{m}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{\Phi }\mathit{,}\hspace{0.17em}\frac{{\mathit{L}}_{\mathit{x}}}{\mathit{m}}\mathrm{=}\mathrm{-}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{\theta }}\sin \mathit{\varphi }\mathrm{-}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{\varphi }}\cos \mathit{\theta }\cos \mathit{\varphi ,}\\ \frac{{\mathit{L}}_{\mathit{y}}}{\mathit{m}}& \mathrm{=}& \end{array}$(20)where E and L_i(i = x,y,z) are the energy and the projection of the angular momentum on the corresponding axis. From the last three integrals follows the conservation of the absolute value of the angular momentum L, written in the form $\begin{array}{ccc}\frac{{\mathit{L}}^{\mathrm{2}}}{{\mathit{m}}^{\mathrm{2}}}\mathrm{=}{\mathit{r}}^{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{·}& & \end{array}$(21)Then, the solution of the kinetic Eq. (18) is an arbitrary function of the first integrals (20). We restrict ourselves to an isotropic distribution function f(E). For a uniform DE, a normalization of its energy at r = ∞ is not possible, therefore we choose Φ_Λ = 0 at r = 0 as the most convenient one (Merafina et al. 2012). Thus, from Eqs. (15), we have $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{\Lambda }}\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi G}}{\mathrm{3}}{\mathit{\rho }}_{\mathrm{\Lambda }}{\mathit{r}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}}{\mathrm{6}}{\mathit{r}}^{\mathrm{2}}\mathrm{·}& & \end{array}$(22)Following Zel’dovich & Podurets (1965) and Bisnovatyi-Kogan et al. (1993, 1998), we consider a Maxwell-Boltzmann distribution function with a cut-off $\begin{array}{ccc}\{\begin{array}{c}\\ & \\ & \\ & \end{array}& & \end{array}$(23)where the cut-off energy E_cut is given by $\begin{array}{ccc}{\mathit{E}}_{\mathrm{cut}}\mathrm{=}\mathrm{-}\frac{\mathit{\alpha T}}{\mathrm{2}}& & \end{array}$(24)and α is the so-called cut-off parameter, while T is the temperature in energy units. The total energy is $\begin{array}{ccc}\mathit{E}\mathrm{=}\frac{\mathit{m}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{m}\mathrm{\Phi }\mathrm{=}\frac{\mathit{m}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\mathit{m}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{-}\frac{\mathit{m}\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}}{\mathrm{6}}\mathit{,}& & \end{array}$(25)where the total potential Φ and the velocity v are given by $\begin{array}{ccc}\mathrm{\Phi }\mathrm{=}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{\Lambda }} \mathrm{and} \mathit{v}\mathrm{=}\mathrm{(}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(26)The constant B in the first of Eqs. (23) depends on the total potential Φ and therefore is different for each model. To consider a unique distribution function for all the equilibrium configurations, following Merafina & Ruffini (1989), we must choose a different normalization by introducing a new constant A connected with B through the following relation¹$\begin{array}{ccc}\mathit{B}\mathrm{=}\mathit{A}{\mathrm{e}}^{\mathit{m}{\mathrm{\Phi }}_{\mathit{R}}\mathit{/}\mathit{T}}\mathit{,}& & \end{array}$(27)with Φ_R the value of the total potential Φ at r = R. In this way, the expression of the distribution function (23) for E ≤ E_cut becomes $\begin{array}{ccc}\mathit{f}\mathrm{=}\mathit{A}\hspace{0.17em}\exp \left[\frac{\mathit{m}{\mathrm{\Phi }}_{\mathit{R}}}{\mathit{T}}\mathrm{-}\frac{\mathit{m}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}\mathit{T}}\mathrm{-}\frac{\mathit{m}}{\mathit{T}}\left({\mathrm{\Phi }}_{\mathrm{g}}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}}{\mathrm{6}}\right)\right]\mathrm{·}& & \end{array}$(28)The maximum kinetic energy ϵ_c is connected with the potential Φ by the relation $\begin{array}{ccc}{\mathit{ϵ}}_{\mathrm{c}}\mathrm{=}\mathit{m}\mathrm{(}{\mathrm{\Phi }}_{\mathit{R}}\mathrm{-}\mathrm{\Phi }\mathrm{)}\mathit{.}& & \end{array}$(29)Then the distribution function can be rewritten as $\begin{array}{ccc}\mathit{f}\mathrm{=}\mathit{A}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{ϵ}\mathrm{-}{\mathit{ϵ}}_{\mathrm{c}}\mathrm{)}\mathit{/}\mathit{T}}\mathit{,}& & \end{array}$(30)where ϵ = mv²/ 2 is the kinetic energy of the single-point mass.

The Poisson Eq. (15) in a spherical symmetry applied to a gravitational field is given by $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left({\mathit{r}}^{\mathrm{2}}\frac{\mathrm{d}{\mathrm{\Phi }}_{\mathrm{g}}}{\mathrm{d}\mathit{r}}\right)\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{m}}\mathit{,}& & \end{array}$(31)with the boundary conditions Φ_g(0) = Φ_g0 and ${\mathrm{\Phi }}_{\mathrm{g}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$. The matter density can be expressed as $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\mathrm{4}\mathit{\pi m}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{p}}_{\max }}\mathit{f}{\mathit{p}}^{\mathrm{2}}\mathrm{d}\mathit{p,} \mathrm{with} \mathit{p}\mathrm{=}\mathit{mv,}& & \end{array}$(32)where the expression of the maximum momentum p_max is given by $\begin{array}{ccc}{\mathit{p}}_{\max }\mathrm{=}\sqrt{\mathrm{2}\mathit{m}\left(\mathrm{-}\mathit{m}\mathrm{\Phi }\mathrm{-}\frac{\mathit{\alpha T}}{\mathrm{2}}\right)}\mathrm{=}\sqrt{\mathrm{2}\mathit{m}\left(\mathrm{-}\mathit{m}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{+}\frac{\mathit{m}\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}}{\mathrm{6}}\mathrm{-}\frac{\mathit{\alpha T}}{\mathrm{2}}\right)}\mathrm{·}& & \end{array}$The cluster with a finite radius is possible only when the following condition is satisfied: $\begin{array}{ccc}\frac{\mathit{\alpha T}}{\mathrm{2}\mathit{m}}\mathit{<}\mathrm{-}{\mathrm{\Phi }}_{\max }\mathit{,} \mathrm{with} {\mathrm{\Phi }}_{\max }\mathit{<}\mathrm{0.}& & \end{array}$Then, by using the form of the distribution given in Eq. (30), we can finally rewrite the matter density as $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\mathrm{4}\sqrt{\mathrm{2}}\mathit{\pi A}{\mathit{m}}^{\mathrm{5}\mathit{/}\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{ϵ}}_{\mathrm{c}}}{\mathrm{e}}^{\mathrm{-}\mathrm{(}\mathit{ϵ}\mathrm{-}{\mathit{ϵ}}_{\mathrm{c}}\mathrm{)}\mathit{/}\mathit{T}}\sqrt{\mathit{ϵ}}\hspace{0.17em}\mathrm{d}\mathit{ϵ,} \mathrm{where} {\mathit{ϵ}}_{\mathrm{c}}\mathrm{=}\frac{{\mathit{p}}_{\max }^{\mathrm{2}}}{\mathrm{2}\mathit{m}}\mathrm{·}& & \end{array}$(33)Introducing dimensionless variables $\begin{array}{ccc}\mathit{W}\mathrm{=}\frac{{\mathit{p}}_{\max }^{\mathrm{2}}}{\mathrm{2}\mathit{mT}}\mathrm{=}\frac{{\mathit{ϵ}}_{\mathrm{c}}}{\mathit{T}} \mathrm{and} \mathit{x}\mathrm{=}\frac{{\mathit{p}}^{\mathrm{2}}}{\mathrm{2}\mathit{mT}}\mathrm{=}\frac{\mathit{ϵ}}{\mathit{T}}\mathit{,}& & \end{array}$(34)we obtain W = m(Φ_R − Φ) /T, and the expression of matter density ρ_m becomes $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\mathrm{4}\sqrt{\mathrm{2}}\mathit{\pi A}{\mathit{m}}^{\mathrm{5}\mathit{/}\mathrm{2}}{\mathit{T}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{W}}{\mathrm{e}}^{\mathit{W}\mathrm{-}\mathit{x}}\sqrt{\mathit{x}}\mathrm{d}\mathit{x,}& & \end{array}$(35)where, as usual, at W = 0 we have ρ_m = 0, being Φ = Φ_R.

From the Poisson Eq. (31) we can deduce the equation describing the structure of the Newtonian configurations in presence of DE by also considering the potential Φ_Λ. In fact, inserting the expression of the gravitational potential Φ_g = Φ − Φ_Λ into Eq. (31) and using Eq. (22), we obtain $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left({\mathit{r}}^{\mathrm{2}}\frac{\mathrm{d\Phi }}{\mathrm{d}\mathit{r}}\right)\mathrm{=}\mathrm{4}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{m}}\mathrm{-}\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}\mathit{,}& & \end{array}$(36)where the potential Φ now includes all the contributions. Then, by considering the first relation in Eq. (11), the equilibrium equation becomes $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left({\mathit{r}}^{\mathrm{2}}\frac{\mathrm{d\Phi }}{\mathrm{d}\mathit{r}}\right)\mathrm{=}\mathrm{4}\mathit{\pi G}\mathrm{(}{\mathit{\rho }}_{\mathrm{m}}\mathrm{-}\mathrm{2}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{)}\mathit{.}& & \end{array}$(37)Now, we have to consider the boundary conditions for the potential Φ with respect to the conditions given for the potential Φ_g in Eq. (31). Starting from Eq. (22), we can write $\begin{array}{ccc}\mathrm{\Phi }\mathrm{=}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}}{\mathrm{6}}{\mathit{r}}^{\mathrm{2}} \mathrm{and} {\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{g}}^{\mathrm{\prime }}\mathrm{-}\frac{\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}}{\mathrm{3}}\mathit{r,}& & \end{array}$(38)and therefore, for r = 0, we have Φ(0) = Φ_g0 and ${\mathrm{\Phi }}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{g}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$.

To write the dimensionless form of the equilibrium equation, we can express the radial coordinate as and, using the definition W = m(Φ_R − Φ) /T, the equilibrium equation can be rewritten as $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r̂}}\left(\mathit{r̂}\mathrm{2}\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\mathit{r̂}}\right)\mathrm{=}\mathrm{-}\frac{\mathrm{4}\mathit{\pi Gm}{\mathit{\eta }}^{\mathrm{2}}}{\mathit{T}}\mathrm{(}{\mathit{\rho }}_{\mathrm{m}}\mathrm{-}\mathrm{2}{\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{)}\mathrm{·}& & \end{array}$(39)In the same way, following Merafina & Ruffini (1989), we can introduce the expression of dimensionless densities by defining the following quantities $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathit{G}{\mathit{\eta }}^{\mathrm{2}}}\mathit{\rho ̂}\mathrm{m} \mathrm{and} {\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{\mathit{G}{\mathit{\eta }}^{\mathrm{2}}}\mathit{\rho ̂}\mathrm{\Lambda }\mathit{,}& & \end{array}$(40)where σ² = 2T/m. Thus, the dimensionless form of the equilibrium equation will be given by $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r̂}}\left(\mathit{r̂}\mathrm{2}\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\mathit{r̂}}\right)\mathrm{=}\mathrm{-}\mathrm{8}\mathit{\pi }\mathrm{(}\mathit{\rho ̂}\mathrm{m}\mathrm{-}\mathrm{2}\mathit{\rho ̂}\mathrm{\Lambda }\mathrm{)}\mathit{,}& & \end{array}$(41)with the boundary conditions W(0) = W₀ and W′(0) = 0. Moreover, it is important to note that the relation must be satisfied at the centre of the equilibrium configuration to obtain the condition of initial decreasing density W′′(0) < 0. However, this is a necessary but not sufficient condition for the existence of the equilibrium solution, because the presence of the DE can enable conditions of increasing density (W′> 0) to be reached at other values of the radial coordinate.

We still need to define the expression of the dimensional quantity η. To derive the result, we can use the relations (35) and (11) for the densities ρ_m and ρ_Λ, respectively, and compare them with the definitions (40). We obtain $\begin{array}{ccc}\mathit{\eta }\mathrm{=}{\mathrm{(}\mathit{A}{\mathit{m}}^{\mathrm{4}}\mathit{G}{\mathit{\sigma }}^{\mathrm{)}}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(42)with $\begin{array}{ccc}\mathit{\rho ̂}\mathrm{m}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{W}}{\mathrm{e}}^{\mathit{W}\mathrm{-}\mathit{x}}\sqrt{\mathit{x}}\mathrm{d}\mathit{x} \mathrm{and} \mathit{\rho ̂}\mathrm{\Lambda }\mathrm{=}\frac{\mathrm{\Lambda }{\mathit{\eta }}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}\mathit{,}& & \end{array}$(43)where is given by the value of Λ. The total mass M^(m) at radius R is given by $\begin{array}{ccc}{\mathit{M}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}{\mathit{\rho }}_{\mathrm{m}}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathit{\eta }}{\mathit{G}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R̂}}\mathrm{4}\mathit{\pi }\mathit{\rho ̂}\mathrm{m}\mathit{r̂}\mathrm{2}\mathrm{d}\mathit{r̂}\mathit{,}& & \end{array}$(44)where $\begin{array}{ccc}\mathit{M̂}\mathrm{\left(}\mathrm{m}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R̂}}\mathrm{4}\mathit{\pi }\mathit{\rho ̂}\mathrm{m}\mathit{r̂}\mathrm{2}\mathrm{d}\mathit{r̂} \mathrm{and} {\mathit{M}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathit{\eta }}{\mathit{G}}\mathit{M̂}\mathrm{\left(}\mathrm{m}\mathrm{\right)}\mathrm{·}& & \end{array}$(45)Finally, to make the dependence of the dimensional quantities on the velocity σ explicit, we can introduce the quantity $\begin{array}{ccc}\mathit{\zeta }\mathrm{=}\mathit{\eta }\hspace{0.17em}{\mathit{\sigma }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{(}\mathit{A}{\mathit{m}}^{\mathrm{4}}\mathit{G}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(46)and the dimensional quantities can be rewritten as $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{3}}}{\mathit{G}{\mathit{\zeta }}^{\mathrm{2}}}\mathit{\rho ̂}\mathrm{m} \mathrm{and} {\mathit{\rho }}_{\mathrm{\Lambda }}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{3}}}{\mathit{G}{\mathit{\zeta }}^{\mathrm{2}}}\mathit{\rho ̂}\mathrm{\Lambda }& & \end{array}$(47)and $\begin{array}{ccc}{\mathit{M}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{\zeta }}{\mathit{G}}\mathit{M̂}\mathrm{\left(}\mathrm{m}\mathrm{\right)} \mathrm{and} \mathit{R}\mathrm{=}\frac{\mathit{\zeta }}{{\mathit{\sigma }}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{R̂}\mathrm{·}& & \end{array}$(48)Turning to the condition (24) on the energy E_cut, we can express the cut-off parameter α by using the condition at the edge of the configuration $\begin{array}{ccc}\frac{\mathit{\alpha }}{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{\mathit{m}{\mathrm{\Phi }}_{\mathit{R}}}{\mathit{T}}\mathrm{=}\mathrm{-}\frac{\mathit{m}}{\mathit{T}}\mathrm{(}{\mathrm{\Phi }}_{\mathrm{g}}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{\Lambda }}{\mathrm{)}}_{\mathit{r}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{R}}\mathrm{·}& & \end{array}$(49)Thus, because Φ_g(R) = −GM^(m)/R and Φ_Λ(R) = −Λc²R²/ 6, we obtain $\begin{array}{ccc}\mathit{\alpha }\mathrm{=}\frac{\mathrm{2}\mathit{Gm}{\mathit{M}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{\mathit{RT}}\mathrm{+}\frac{\mathit{m}\mathrm{\Lambda }{\mathit{c}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{2}}}{\mathrm{3}\mathit{T}}& & \end{array}$(50)and, finally, by using dimensionless quantities (47), (48) and relation (3), we have $\begin{array}{ccc}\mathit{\alpha }\mathrm{=}\frac{\mathrm{4}\mathit{M̂}\mathrm{\left(}\mathrm{m}\mathrm{\right)}}{\mathit{R̂}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{4}\mathit{\pi }\mathit{\rho ̂}\mathrm{\Lambda }\mathit{R̂}\mathrm{3}}{\mathrm{3}\mathit{M̂}\mathrm{\left(}\mathrm{m}\mathrm{\right)}}\right)\mathrm{·}& & \end{array}$(51)For low values of the cut-off parameter α, maintaining a finite value of αT that corresponds to high values of the temperature T, the distribution function (23) may be taken as a constant (Bisnovatyi-Kogan et al. 1998). Then, the solutions only exist for low values of W₀ and , assuming a more simplified form that converges to a limiting sequence. In the limit of W → 0, the dimensionless density ρ_m can be expressed as $\begin{array}{ccc}\mathit{\rho ̂}\mathrm{m}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{W}}\sqrt{\mathit{x}}\mathrm{d}\mathit{x}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{W}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(52)whereas the equilibrium Eq. (41) becomes $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{r̂}\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r̂}}\left(\mathit{r̂}\mathrm{2}\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\mathit{r̂}}\right)\mathrm{=}\mathrm{-}\frac{\mathrm{32}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{3}}{\mathit{W}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{16}\mathit{\pi }\mathit{\rho ̂}\mathrm{\Lambda }\mathrm{·}& & \end{array}$(53)Expressed in dimensional terms, the density can be written as $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}\frac{{\mathit{\sigma }}^{\mathrm{3}}}{\mathit{G}{\mathit{\zeta }}^{\mathrm{2}}}{\mathit{W}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{=}{\mathit{\rho }}_{\mathit{p}}{\mathit{W}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(54)where $\begin{array}{ccc}{\mathit{\rho }}_{\mathit{p}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}\frac{{\mathit{\sigma }}^{\mathrm{3}}}{\mathit{G}{\mathit{\zeta }}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(55)Moreover, by imposing a change of radial coordinate from to y for which $\begin{array}{ccc}\mathit{r̂}\mathrm{=}\mathit{y}\hspace{0.17em}{\left(\frac{\mathrm{3}}{\mathrm{32}{\mathit{\pi }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(56)the dimensionless equilibrium equation can be rewritten as $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{y}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{y}}\left({\mathit{y}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\mathit{y}}\right)\mathrm{=}\mathrm{-}{\mathit{W}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}\mathit{\pi }}\mathit{\rho ̂}\mathrm{\Lambda }\mathrm{·}& & \end{array}$(57)We can also substitute the density by using Eqs. (47) and (55) and finally obtain $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{y}}^{\mathrm{2}}}\frac{\mathrm{d}}{\mathrm{d}\mathit{y}}\left({\mathit{y}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\mathit{y}}\right)\mathrm{=}\mathrm{-}{\mathit{W}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\frac{\mathrm{2}{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\mathit{\rho }}_{\mathit{p}}}\mathit{,}& & \end{array}$(58)which corresponds, if we take ρ_p = ρ_m0 and W ≡ θ, exactly to the equilibrium equation for a polytropic configuration with index n = 3/2 in presence of DE introduced by Merafina et al. (2012) in accordance with the dimensionless Emden variables and the initial conditions W(0) = θ(0) = 1 and W′(0) = θ′(0) = 0. Therefore, the polytropic configurations calculated by Merafina et al. (2012) in hydrostatic approach can be used to describe clusters of gravitating point masses with distribution function (28) with the energy cut-off (29), at low values of Λ and high values of T.

4. Numerical results

The dimensionless equilibrium Eq. (41) depends on two parameters: the gravitational potential at the centre of configurations W₀ and , which determines the intensity of DE through the value of the cosmological constant Λ. Different values of these parameters give a two-dimensional family of equilibrium solutions. The set of solutions for at different values of W₀ was obtained by Bisnovatyi-Kogan et al. (1998).

We solved numerically the Poisson equation for gravitational equilibrium at different values of the two parameters (W₀, ) mentioned above. First of all, we focused our attention on the matter density profiles ρ_m(r) of the equilibrium configurations; in detail, we investigated how they change for increasing values of the dimensionless DE density at fixed values of the dimensionless gravitational potential W₀. We chose three values of W₀ and four values of , which are the basis of pairs of parameters (W₀, ) that do not allow equilibrium solutions. This peculiarity is clearly represented in the matter density profiles shown in Figs. 1–3. For conciseness, we define the following quantitites in the figures: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{\ast }}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{3}}}{\mathit{G}{\mathit{\zeta }}^{\mathrm{2}}}\mathit{,} {\mathit{M}}_{\mathrm{\ast }}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{\zeta }}{\mathit{G}}\mathit{,} {\mathit{R}}_{\mathrm{\ast }}\mathrm{=}\frac{\mathit{\zeta }}{{\mathit{\sigma }}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,}& & \end{array}$(59)and, therefore, the dimensionless quantities introduced in Eqs. (47) and (48) can be rewritten as $\begin{array}{ccc}\mathit{\rho ̂}\mathrm{m}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{\rho }}_{\mathrm{\ast }}}\mathit{,}\mathit{\rho ̂}\mathrm{\Lambda }\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{\Lambda }}}{{\mathit{\rho }}_{\mathrm{\ast }}}\mathit{,}\mathit{M̂}\mathrm{\left(}\mathrm{m}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{{\mathit{M}}_{\mathrm{\ast }}}\mathit{,}\mathit{R̂}\mathrm{=}\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{\ast }}}\mathit{,}\mathit{r̂}\mathrm{=}\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{·}& & \end{array}$(60)For each value of the central potential W₀ there is one value of the parameter after which the matter density profile does not converge to zero, but oscillates indefinitely. If we assert that the radius R of an equilibrium configuration is defined as the value of the radial coordinate r at which the matter density ρ_m(r) becomes zero, it is clear that every time this does not occur, we are unable to estimate the radial extension of the system. All configurations with a given value of W₀ and that correspond to oscillating density profiles cannot be considered in gravitational equilibrium.

Fig. 1 Dimensionless matter density profiles for equilibrium configurations with W0 = 4 and (dashed line), 0.5 (dotted line), 0.9 (dash-dotted line), for which the matter density converges to zero, and (solid line), chosen inside the region of non-equilibrium solutions.

Fig. 2 Dimensionless matter density profiles for equilibrium configurations with W0 = 8 and 0 (dashed line), 0.3 (dotted line), 0.5 (dash-dotted line), for which the matter density converges to zero, and 0.8 (solid line), chosen inside the region of non-equilibrium solutions.

Fig. 3 Dimensionless matter density profiles for equilibrium configurations with W0 = 12 and (dashed line), 0.3 (dotted line), 0.5 (dash-dotted line), for which the matter density converges to zero, and (solid line), chosen inside the region of non-equilibrium solutions.

Moreover, the calculation of the total radius R of an equilibrium configuration is strictly connected to the calculation related to the total mass M^(m). Following Eq. (45) and expressed in terms of dimensionless quantities, the mass $\mathit{M̂}\begin{array}{c}\mathrm{\left(}\mathrm{m}\mathrm{\right)}\\ \mathit{r}\end{array}$ within the radius is given by $\begin{array}{ccc}\mathit{M̂}\begin{array}{c}\mathrm{\left(}\mathrm{m}\mathrm{\right)}\\ \mathit{r}\end{array}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{r̂}}\mathrm{4}\mathit{\pi }\mathit{\rho ̂}\mathrm{m}{\mathit{\xi }}^{\mathrm{2}}\mathrm{d}\mathit{\xi }\mathit{.}& & \end{array}$(61)As a consequence, the non-equilibrium solutions for which the total radius R cannot be defined do not even allow evaluating the total mass M^(m) of the system.

Bisnovatyi-Kogan et al. (1998) found the set of solutions at Λ = 0 for and curves in the Newtonian case. These curves are shown in Figs. 4 and 5 (continuous line) together with the curves given for different values of (Λ ≠ 0). Analysing Fig. 4, when the parameter is different from zero, and for increasing values of this parameter, the curves are no longer continuous and the absolute maximum of the mass disappears. Within the interval , the curves present several branches (in the figure, only the branches for are shown to be able to clearly understand the different behaviours). Out of the interval the branches reduce to a unique curve and, in particular, for the curve becomes gradually shorter at increasing values of until it reaches the critical value , when the curve reduces to a unique point (we discuss this critical value below). It is clear that the unusual behaviour of the is related to the density profiles of the non-equilibrium solutions. To analyse Fig. 5, by considering Eq. (51), we can conclude that the parameter α is also connected to the values of the total radius and the mass . Consequently, it is easy to show that for non-equilibrium solutions it is not possible to calculate the cut-off parameter α. Therefore we can expect the existence of different branches of solutions here as well. Moreover, the behaviour of the curves at different values of extends the range of solutions at values of α higher than the critical value (α = 2.87) valid for Λ = 0 (Bisnovatyi-Kogan et al. 1998). As previously underlined, it is possible to distinguish several branches of solutions with a limiting value of α that changes in dependence of the value of . This limiting value, systematically higher than 2.87, increases at increasing values of until the absolute limiting value α_lim ≃ 3.42.

Fig. 4 Dimensionless mass as a function of the dimensionless central matter density for 0 (solid line), 0.5 (dashed line), 0.8 (dotted lines), 1.3 (dash-dotted line). The discontinuity of the curves for some values of the central matter density arises because the condition R ≤ RΛ that is necessary for the existence of the equilibrium configurations is not fulfilled (see Figs. 1–3 and Eqs. (62), (63)).

Fig. 5 Dimensionless mass as a function of the cut-off parameter α for (solid line), 0.5 (dashed line), 0.8 (dotted lines), 1.3 (dash-dotted line). The discontinuity of the curves for some values of the cut-off parameter arises because the condition R ≤ RΛ that is necessary for the existence of the equilibrium configurations is not fulfilled (see Figs. 1–3 and Eqs. (62), (63)).

The DE background in which all the bodies of the Universe are embedded produces the anti-gravity that changes their gravitational equilibrium, acting in contrast to the matter gravity. To establish when we have found configurations for which the presence of the DE can change the gravitational equilibrium, following Bisnovatyi-Kogan & Chernin (2012), we introduce the so-called zero gravity radius R_Λ. This is a physical parameter that is defined as the distance from the centre of the system where the matter gravity and DE anti-gravity balance each other exactly. We consider the total force acting on the unit mass $\begin{array}{ccc}\mathit{F}\mathrm{=}{\mathit{F}}_{\mathrm{g}}\mathrm{+}{\mathit{F}}_{\mathrm{\Lambda }}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{8}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{\Lambda }}}{\mathrm{3}}\mathit{r,}& & \end{array}$(62)where, differently from Eq. (6), this relation is also valid within the matter and not only in the vacuum. Then, the total force F defined in Eq. (62) and, consequently, the acceleration, are both zero at a distance $\begin{array}{ccc}{\mathit{R}}_{\mathrm{\Lambda }}\mathrm{=}{\left[\frac{\mathrm{3}{\mathit{M}}_{{\mathit{R}}_{\mathrm{\Lambda }}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}}{\mathrm{8}\mathit{\pi }{\mathit{\rho }}_{\mathrm{\Lambda }}}\right]}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}& & \end{array}$(63)where the zero-gravity radius depends on the total mass of the equilibrium configurations if R ≤ R_Λ, while if the condition F = 0 is satisfied inside the configuration, we have no equilibrium, and the mass to consider is ${\mathit{M}}_{\mathit{r}}^{\mathrm{\left(}\mathrm{m}\mathrm{\right)}}$ with r = R_Λ.

This means that every cluster has its zero-gravity radius. This definition allows us to identify a gravitationally bound system only if it is enclosed within the sphere of radius R_Λ, namely only if its total radius is smaller than its zero-gravity radius (R<R_Λ). Galaxies in the external regions where r ≥ R_Λ can flow out from the centre of the cluster under the action of the DE anti-gravity force.

In Fig. 6 we have represented the curve of the equilibrium configurations with R = R_Λ, through the behaviour of W₀ as a function of . In addition, we have also shown the curves that represent the families of equilibrium solutions at fixed values of α. When the value of W₀ is kept constant and the value of is increased, we obtain one limiting value located on the curve after which it is no longer possible to obtain equilibrium solutions. Along this limiting curve, which separates two regions (solid line), the equilibrium configurations have the total radius exactly equal to the zero-gravity radius and the matter density profiles vanishing with a minimum in correspondence of the total radius R = R_Λ. This enables defining the region on the right side of the figure in which no gravitational equilibrium can establish and no curves at constant α can lie, corresponding to configurations with matter density profiles that do not converge to zero. In contrast, in the region corresponding to the left side of the figure, we can assert that the force due to the presence of the DE, F_Λ, is weaker than the force due to the gravity, F_g, and gravitational equilibrium can be achieved. In other words, speaking in terms of radial extension, the condition R ≤ R_Λ is satisfied for each configuration belonging to this region, and the matter density profiles are regular and converging to zero in correspondence to the total radius R.

Fig. 6 Limiting curve of the equilibrium configurations with R = RΛ (solid line), expressed in terms of W0 as a funtion of . Labelled curves at constant α (dashed lines) are also considered.

Finally, from Eq. (63) we can see that the zero-gravity radius R_Λ is inversely proportional to the DE density ρ_Λ. Therefore, by decreasing the value of the DE density, the zero-gravity radius increases until the condition R_Λ → ∞ for ρ_Λ = 0. By considering the plane (W₀-) of Fig. 6, the gravitational equilibrium is even achieved more easily and for more values of W₀, when parameter is small. In particular, for , the equilibrium solutions are possible for each value of W₀, recovering the well-known results of Bisnovatyi-Kogan et al. (1998).

5. Conclusions

We have calculated the equilibrium configurations of Newtonian clusters with a truncated Maxwellian distribution function in presence of DE. All clusters that satisfy the condition R ≤ R_Λ have a structural equilibrium and can be considered dynamically stable. On the other hand, there are conditions for which the effects of DE prevail on the gravity, and equilibrium cannot be reached. This occurs whenever the zero-gravity radius lies inside the configuration and divides the inner part, which is dominated by gravity, from the external part in which the expanding forces due to DE are prevalent.

Then we described the density distribution inside galaxy clusters by several phenomenological functions, some of which follow from numerical simulations (see Chernin et al. 2013). Qualitatively, the truncated Maxwellian distribution considered here is similar to the non-singular density distribution suggested by Chernin et al. (2013). It may be used for a more detailed study of the density and velocity distribution on the periphery of rich clusters, where the influence of DE is significant, and their comparison with observations.

The number density of the galaxies located in the outer part of a cluster is less relevant than those of the central region and, in general, these galaxies have smaller masses and lower luminosities in presence of an even weaker relaxation because of the low probability of encounters. Therefore, only the largest

telescopes should be used to search for galaxies in the external cluster regions. Furthermore, the most sensitive X-ray telescopes are needed to detect the hot gas in the low-density regions to reveal the possible outflow in presence of DE that was considered by Bisnovatyi-Kogan & Merafina (2013).

Finally, the evaluation of the parameters that characterize the clusters of galaxies suggests that these systems are collisionless. In fact, if we consider the relaxation time, we obtain higher values than the age of the Universe and, therefore, we can conclude that thermodynamical instabilities are irrelevant in the current evolution of the galaxy clusters. On the other hand, using well-known criteria for identifying the onset of thermodynamical instability (Bisnovatyi-Kogan & Merafina 2006), we can see that the critical point lies far from the first maximum mass of the curve with in Fig. 4, namely, at higher values of the central matter density, as well as in curves with , which allows us to conclude that the larger part of the equilibrium configurations is thermodynamically stable.

Acknowledgments

The work of G.S. Bisnovatyi-Kogan was partly supported by RFBR grant 11-02-00602, the RAN Program Formation and evolution of stars and galaxies, and Russian Federation President Grant for Support of Leading Scientific Schools NSh-5440.2012.2.

References

All Figures

	Fig. 1 Dimensionless matter density profiles for equilibrium configurations with W0 = 4 and (dashed line), 0.5 (dotted line), 0.9 (dash-dotted line), for which the matter density converges to zero, and (solid line), chosen inside the region of non-equilibrium solutions.
In the text

	Fig. 2 Dimensionless matter density profiles for equilibrium configurations with W0 = 8 and 0 (dashed line), 0.3 (dotted line), 0.5 (dash-dotted line), for which the matter density converges to zero, and 0.8 (solid line), chosen inside the region of non-equilibrium solutions.
In the text

	Fig. 3 Dimensionless matter density profiles for equilibrium configurations with W0 = 12 and (dashed line), 0.3 (dotted line), 0.5 (dash-dotted line), for which the matter density converges to zero, and (solid line), chosen inside the region of non-equilibrium solutions.
In the text

Fig. 4 Dimensionless mass as a function of the dimensionless central matter density for 0 (solid line), 0.5 (dashed line), 0.8 (dotted lines), 1.3 (dash-dotted line). The discontinuity of the curves for some values of the central matter density arises because the condition R ≤ RΛ that is necessary for the existence of the equilibrium configurations is not fulfilled (see Figs. 1–3 and Eqs. (62), (63)).

In the text

	Fig. 5 Dimensionless mass as a function of the cut-off parameter α for (solid line), 0.5 (dashed line), 0.8 (dotted lines), 1.3 (dash-dotted line). The discontinuity of the curves for some values of the cut-off parameter arises because the condition R ≤ RΛ that is necessary for the existence of the equilibrium configurations is not fulfilled (see Figs. 1–3 and Eqs. (62), (63)).
In the text

	Fig. 6 Limiting curve of the equilibrium configurations with R = RΛ (solid line), expressed in terms of W0 as a funtion of . Labelled curves at constant α (dashed lines) are also considered.
In the text