Data Analysis of Globular Cluster Harris Catalogue in view of the King models and their dynamical evolution. I. Theoretical model

We discuss the possibility to analyze the problem of gravothermal catastrophe in a new way, by obtaining thermodynamical equations to apply to a selfgravitating system. By using the King distribution function in the framework of statistical mechanics we treat the globular clusters evolution as a sequence of quasi-equilibrium thermodynamical states.


Introduction
Globular clusters (GCs) are stellar systems with masses within the interval 10 4 − 10 6 M ⊙ , containing a number of stars of the order of 10 5 . They are considered as nearly spherical systems due to their low values of eccentricity e; at least 50% of GCs have e < 0.1 and there are no clusters with e > 0.2. The core radius r c , namely the radial coordinate at which the brightness becomes one half of the corresponding value at the center of the system, is almost 10 pc, whereas the tidal radius r t , which is the biggest spatial extension of the cluster allowed by the external tidal field, is typically around 50 pc.
For their symmetry and age, there is the possibility to test the evolution of a GCs studying a classical single mass King model (King, 1966) in relation to thermodynamical instability phenomena. In fact, in the analysis of the evolution of GCs, stellar encounters strongly contribute in phase space mixing of stellar orbits. In this scenario, thermodynamics plays a centrale role in the gravitational equilibrium and stability of these clusters, being the average binary relaxation time shorter than their old absolute age which ranges between 10 to 13 Gyr.
This means that Fokker-Planck approximation, which takes into account the nature of collisions in globular clusters, can determine the distribution function relevant for obtaining the equilibrium configurations of these systems, whereas the tidal effects due to the presence of Galactic gravitational potential are responsible of the confination of the cluster.
On the other hand, the observations of the luminosity profiles of different GCs (King, 1962) show similar curves depending only on different values of the star concentration, giving the possibility to fit them by an empirical law and suggesting a unique distribution function for the whole sample of clusters (King, 1966). This effect can be described as a change of the main parameters of the cluster during the dynamical evolution (Horwitz & Katz, 1977), which maintains the form of the distribution like in a sort of reversible trasformation of a gas in thermodynamic equilibrium, also in accordance to numerical simulations existing in literature which result in keeping unchanged the distribution of velocities of stars for a wide range of values of the concentration during the time evolution driven by the Fokker-Planck equation.
Therefore, the evolution of globular clusters can be studied by considering small thermodynamic transformations which keep constant the functional form of the velocity distribution of stars like in the framework of Boltzmann statistical mechanics. It is important to note that while the equilibrium is given by the form of the distribution which depends on the Fokker-Planck equation and consider the real nature of collisions, thermodynamics plays a role in the tidal effects acting on the confination of the system, due to a two competitive phenomena: one given by stellar encounters which tend to refresh the tail of high velocities in the distribution and one due to evaporation of stars which prevents the formation of it, maintaining the system in a sort of thermodynamical equilibrium with the same distribution function even if in presence of a cutoff in the velocity of the stars.

The effective potential
The King DF characterizing the energy distribution of stars with the same mass m in a model that describes a spherically symmetric system with isotropic velocity distribution may be written as Here ψ = m (ϕ R − ϕ) is the energy cutoff, corresponding to the maximum kinetic energy that a star can have at a given radial coordinate r, while ϕ is the gravitational potential. This energy is sufficient to reach the border of the equilibrium configuration r = R, being also the difference between the value of the gravitational potential at the edge of configuration and the same quantity evaluated at a generic distance r from the center. The quantity θ is the thermodynamic temperature of the system while B is a constant of normalization. The behavior of equilibrium solutions for King models has been also analyzed by Merafina & Ruffini (1989) by solving the Poisson equilibrium equation in Newtonian regime. We can see the presence of a maximum value of the total mass M at increasing values of the central gravitational potential W 0 , which denotes the arising of a sort of thermodynamic instability at W 0 = 1.35 (Fig.1). Over this value, we can think the system can evolve towards the loss of thermodynamical equilibrium (gravothermal catastrophe), in accordance to the expected evolution of Lynden- Bell & Wood (1968).
In order to consider thermodynamic transformations in the framework of statistical mechanics, it is possible to describe the King DF like a Maxwell-Boltzmann one by introducing an effective potential.
In this way the evolution of the King models can be treated as a succession of quasi-equilibrium stages by a thermodynamic theory formally equivalent to the classical one.
Then, the expression of the effective potential is given by (2) and the distribution function can be expressed as where H = ε + mϕ + φ is the single particle Hamiltonian of the system which includes also the gravitational energy of the single star. The effective potential is a screen potential which restricts the phase space of the available velocities for the stars and takes into account the effect of the tidal forces on the system. In this way, the kinetic temperature T connected with the average velocity of the stars, depending on the radial coordinate r, becomes distinguished from the thermodynamic temperature θ, constant all over the equilibrium configuration. From the modified Boltzmann DF of Eq.3, we can deduce the generalized thermodynamical quantities, as the energy U , the thermodynamical pressure Π and the entropy S, related to a shell with radial coordinate r. We get where we have replaced the costant B with A, being B = Ae α/kθ and now f = Ae (α−H)/kθ , while α = µ + mϕ is the chemical potential in presence of the gravitational potential ϕ.
In this way we can rewrite the first law of thermodynamics and obtain a new form for the Eulero expression, that include the extensive and intensive quantities. We can get also an equation of state formally equivalent to classical one which involves the thermodynamical quantities, valid for a shell with radial coordinate r. We have and, for kinetic quantities like temperature T and pressure P , Finally, by integrating the expression of U containing the single particle Hamiltonian H (Eq.5) all over the configuration, we can find an additional term E ef f in the expression of the total energy E tot of the system, called effective energy where E kin and E gr are the total kinetic energy and the total gravitational energy, respectively. The partecipation of the effective potential in the total energy corresponds to the account of the tidal potential which determines a finite radius of the cluster. Moreover, Eq.5 defines the energy of the test shell but, for calculating E gr , we need to use the expression

The gravothermal catastrophe
Thermodynamical instability of a selfgravitating spherical system was first studied by Lynden-Bell & Wood (1968), by considering an isothermal sphere (core) confined in a spherical box. Using the classical form of the virial theorem, including a boundary term due to spatial truncation of the density prophile, it is possible, for that system, to calculate the critical value of the central gravitational potential W 0 = 6.55 after that thermodynamical instability, known as gravothermal catastrophe, onsets. It is important to note that such instability takes place only in presence of an external thermal bath exchanging heat with the core and driving the system towards the dynamical collapse.
With the introduction of the effective potential, we can repeat this analysis for King models and get another critical value for the central gravitational potential W 0 = 6.9, which differs from the one obtained by Katz (1980) W 0 = 7.4, due to the additional term in the total energy E tot (see Eq.12). The most interesting results concern the profile of specific heat for different values of W 0 (see Merafina et al., in preparation). By analyzing the behavior of the specific heat all over the configuration, we found different results. The expression of the specific heat C V = (dQ/dθ) V arises from Eq.8, being constant N and V , by using the expression • For W 0 < 1.35, we have equilibrium configurations with positive heat capacity all over the system. There are not existing conditions for an evolution of the system towards the critical value corresponding to the onset of the gravothermal catastrophe (W 0 = 6.9). Further, this particular value (W 0 = 1.35) corresponds to one concerning the first maximum mass we found among the equilibrium solutions (see Fig.1).
• For W 0 > 1.35, the system shows an external halo with negative heat capacity and an internal core with a positive value. The system can evolve by increasing the value of W 0 until reaching the critical value in which the gravothermal instability onsets. These evolution can take place without the necessity of the presence of an external thermal bath, differently from the previously requested condition in the Lynden Bell & Wood model.
Results showing the specific heat profiles in function of the radial coordinate for different values of W 0 are summarized in Fig.2.

Preliminar observational evidences
The stability of the King models was analyzed in detail by Katz in 1980, with the same investigation carried out by Lynden-Bell & Wood for the isothermal sphere. Katz introduced a new parameter K, which corresponds essentially to the ratio between the escape velocity and the dispersion velocity, both calculated at the center of the cluster. This parameter is directly connected with W 0 . Calculations performed by Katz showed that models become thermodynamically unstable over the value K = 8.1 (W 0 = 7.4). But, analyzing the sample of data coming from Peterson & King (1975) and Peterson (1976), Katz highlighted an unexplainable gap between the expected value of the sample, K = 7.8 equivalent to W 0 = 6.9, and the one corresponding to the onset of gravothermal instability, resulting at W 0 = 7.4 (see Fig.3). We indeed expect that GCs had enough time to undergo the gravothermal catastrophe and, therefore, the distribution of GCs in terms of K or W 0 should peak exactly in correspondence to the critical value. In fact, the primeval Gaussian distribution, approaching the critical value during the evolution, deforms in a non-symmetric Gaussian curve due to the effect of gravothermal catastrophe which progressively subtracts the collapsed GCs with values of W 0 larger than the critical value. For these reasons, the resulting distribution must present a maximum which corresponds to the critical value.
It is remarkable to note that, with the introduction of the effective potential in the study of the thermodynamical instability, we obtain a critical value, W 0 = 6.9, that bridges this gap and corresponds exactly to the expected value of the sample. This correspondence becomes much more evident by considering the sample of Harris (1996) with 127 clusters (if we exclude the PCC ones), as well becomes more evident the non-symmetric form of the distribution. On the other hand, by making a z-test in order to verify the statistical significance of the gap between the stability limit W 0 = 7.4 expected by Katz and the peak value of the distribution at W 0 = 6.9, it can be shown that these two values are not compatible within a confidence level of 95%.

Conclusions
• The additional (positive) contribution of the effective potential on the total energy, considering also the virial condition 2E kin + E gr = 0, implies that E tot = −E kin + E ef f . This enables us to construct models in which the core has a positive heat capacity, allowing to assume the possibility of a survival of the system from the gravothermal cathastrophe which could explain the existence of post core-collapsed objects (PCC).
• The model is selfconsistent and admits regions with positive and negative heat capacity which can exchange energy and produce gravothermal instability, without the necessity to assume an external bath as in the Lynden-Bell & Wood model.
• We obtain a new critical value for the onset of gravothermal instability by the presence of the effective potential. This value coincides with the value of K (or, equivalently, to W 0 ) corresponding to the peak of the GCs distribution, removing the unexplainable difference outlined by Katz. This is an observational evidence of the effects due to the presence of the effective potential, confirmed in the analysis of data of more than 150 GCs contained in the last version of catalogue recently published by Harris in 2010 (see also Harris, 1996).
Finally, it may be useful to consider some unsolved problems and perspectives in order to develop the analysis of thermodynamical instabilities of GCs.
• The model is not a multimass one and does not take into account the effects in the formation of binary stars. At moment this is a preliminary model which has to be improved.
• The new possibility of measuring transverse velocities of the stars in GCs opens important perspectives on the knowledge of the distribution of the star orbits and their eccentricity, in order to better develop N-body simulations in supporting the validity of the model.