Entropy Maximization, Cutoff Distribution, and Finite Stellar Masses

Abstract

Conventional equilibrium statistical mechanics of open gravitational systems is known to be problematical. We first recall that spherical stars/galaxies acquire unbounded radii, become infinitely massive, and evaporate away continuously if one uses the standard Maxwellian distribution 𝑓^𝐵 (which maximizes the usual Boltzmann-Shannon entropy and hence has a tail extending to infinity). Next, we show that these troubles disappear automatically if we employ the exact most probable distribution 𝑓 (which maximizes the combinatorial entropy and hence possesses a sharp cutoff tail). Finally, if astronomical observation is carried out on a large galaxy, then the Poisson equation together with thermal de Broglie wavelength provides useful information about the cutoff radius 𝑟_𝐾, cutoff energy 𝜀_𝐾, and the huge quantum number 𝐾 up to which the cluster exists. Thereby, a refinement over the empirical lowered isothermal King models, is achieved. Numerically, we find that the most probable distribution (MPD) prediction fits well the number density profile near the outer edge of globular clusters.

1. Introduction

It is well known that standard methods [1, 2] of equilibrium statistical mechanics run into conceptual difficulties when applied to gravitational bound systems [3–11]. Basically, these troubles arise due to the peculiar behaviour of the gravitational interaction (either the pair potential or the mean field) at short or long distances. The aim of the present paper is to focus attention on the triple problems, namely, unbounded radius, infinite mass, and continuous evaporation of every stellar/galactic system described by the conventional Maxwell-Boltzmann (hereinafter referred to as the MB) distribution.

Section 2 below points out that since the MB function 𝑓^𝐵 maximizes only the simple-minded Boltzmann-Shannon entropy, its tail becomes illogical in the energy cells of small occupancy. The ensuing problems of the Maxwellian distribution cannot be really overcome by using ad hoc prescriptions such as enclosing the system in a hypothetical box [7] or modifying the Maxwellian form empirically by invoking gravitational tidal cutoff [4, 8]. Next, Section 3 presents a detailed derivation of our most probable distribution (MPD) f by taking hints from a preliminary investigation by Menon and Agrawal [12] in the molecular context and by Menon et al. [13] in the cosmological context. Such f maximizes rigorously the more sophisticated combinatorial entropy and the corresponding variational conditions dictate that f must possess a sharply truncated tail. Next, Section 4 demonstrates how our MPD idea applied to cosmology resolves the aforesaid troubles of the MB formalism, and how the Poisson equation brings additional features into our theory. We feel that the MPD philosophy may have bright applicational prospects in fitting cosmological data such as the classic study of stellar number densities in globular clusters done by King [14, Figure 2] and the important measurements performed by van Loon et al. [15, Figure 6] showing velocity distribution on the post-mail-sequence stars in 𝜔 Centauri. Finally, the paper ends by presenting several concluding remarks in Section 5 where some other approaches to the subject (namely, self-consistent Hartree calculations, incomplete relaxation in low-density tail, canonical ensemble treatment of virialization, occurrence of a stellar mass spectrum in real gravitating systems, etc.) are also mentioned.

Some related aspects of algebraic interest are reported in two useful appendices. Careful study of Sections 2 and 3 will reveal that quantum mechanical discretization of the single-particle levels is very convenient for setting up the combinatorial entropy and in finding the cutoff number; hence for the sake of ready reference we collect in Appendix A several known formulae concerning semiclassical one-body spectrum as well as the energy cell occupation number 𝜈. Also, a detailed treatment of our variational conditions in Section 4 requires that 𝜈! be replaced by Γ(𝜈+1) everywhere (even for 𝜈≪1); hence Appendix B tells why derivatives of factorials or gamma functions can be readily taken even in the cells of small occupation numbers.

2. Difficulties with the MB Distribution

2.1. Preliminaries

This section begins by quickly recalling a standard derivation of the famous Maxwell-Boltzmann distribution 𝑓^𝐵 in equilibrium statistical mechanics. Particles are assumed to be moving in 𝐷 spatial dimensions at temperature 𝑇 under the influence of a mean field potential energy 𝑊(𝑟). The one-body energy spectrum is divided into 𝐽 cells, particles are distributed at random over these, those in the 𝑗th cell are regarded as mutually identical, and the simple-minded Boltzmann-Shannon entropy functional 𝑆^𝐵=−𝑘^𝐽_𝑗=1𝑔_𝑗𝑓^𝐵_𝑗ln𝑓^𝐵_𝑗−𝑓^𝐵_𝑗+𝑘𝛼^𝐽_𝑗=1𝑔_𝑗𝑓^𝐵_𝑗−𝑁−𝑘𝛽^𝐽_𝑗=1𝑔_𝑗𝑓^𝐵_𝑗𝜀_𝑗−𝐸 (1) is set up. Here 𝑘 is the Boltzmann constant, 𝑔_𝑖 the cell degeneracy, 𝑁 the total number of particles, 𝐸 the total energy, and 𝛼 and 𝛽 are Lagrange multipliers (see Appendix A for precise definitions of various symbols). Next, one maximizes 𝑆^𝐵 with respect to 𝑓^𝐵_𝑗,𝛼, and 𝛽 to arrive at the MB solution 𝑓^𝐵=𝑒^{𝛼−𝛽𝜀};𝜈^𝐵=𝑔𝑓^𝐵:𝜀₀≤𝜀≤∞, (2) where the index j has been dropped in the quasicontinuum limit, 𝜀₀ is the ground level, and the upper end of the simple-particle energy spectrum has been extended to +∞ both for confining as well as nonconfining potentials 𝑊(𝑟). Although (2) has been widely applied [1, 2] to gases/liquids kept in the laboratory, yet its application to open astronomical systems leads to the following serious conceptual puzzles.

2.1.1. Entropy

In the case of gravitational systems, one always looks for the local (not global) maxima of the entropy functional. The MB solution (2) does this job exactly for the Boltzmann-Shannon entropy defined by (1), but only approximately for the more sophisticated combinatorial entropy defined by (10) later. It will be shown in Section 3 that the tail of the MB solution becomes illogical in the energy cells of small occupation numbers.

2.1.2. Density

If (2) is inserted back into the general expression (A.7) of Appendix A for the mass density 𝜌(𝑟), one obtains the famous Boltzmann barometric formula 𝜌^𝐵(𝑟)∝exp−𝑊(𝑟)𝑘𝑇∶𝑟₀≤𝑟≤∞. (3) The attractive short-distance behaviour of 𝑊(𝑟) cannot pose a real problem because the size 𝑟₀ of the quantum ground level 𝜀₀ is finite [6]. But the long-distance behaviour of 𝑊(𝑟) is problematical as regards astrophysics in 𝐷=3 dimensions. Indeed, for a dilute gaseous star [2, page 114] without the Poisson equation constraint, one finds asymptotically 𝑊(𝑟)∼_𝑟→∞−0;𝜌^𝐵(𝑟)∼_𝑟→∞constant∶Gaseous. (4) Also, for the isothermal Emden sphere [3, 8] subject to the Poisson equation constraint, one knows that 𝑊(𝑟)∼_𝑟→∞𝑟ln²𝑎²;𝜌^𝐵(𝑟)∼_𝑟→∞𝑟⁻²∶Emden (5) with a being the isothermal length scale. Clearly, as 𝑟→∞,the nil/slow decrease of 𝜌^𝐵 in (4) and (5) and the logarithmic increase of W in (5) are unphysical.

2.1.3. Radius

From the MB density (3), one computes the mean size ⟨𝑟⟩^𝐵 of the system via ⟨𝑟⟩^𝐵=∫^∞₀𝜕𝑟𝑟³𝜌^𝐵∫^∞₀𝜕𝑟𝑟²𝜌^𝐵=∞ (6)whichdiverges both for gaseous stars (4) and Emden spheres (5).

2.1.4. Mass

The total mass of the MB system is calculated from 𝑀^𝐵=4𝜋^∞₀𝜕𝑟𝑟²𝜌^𝐵(𝑟)=∞ (7) which also diverges for the two cases mentioned above. Thus, in the Boltzmann-Shannon view, the most likely state of an isotropic stellar system has infinite mass.

2.1.5. Evaporation

Since all regions of the phase space up to 𝑟=∞ are allowed an open MB system, for example, the dilute gaseous star goes on evaporating with time, producingas a thereby a net outgoing flux 𝜙^𝐵 of particles [2, page 114] at every positive thermodynamic temperature T: 𝜑^𝐵>0. (8) Of course, the isothermal sphere can be stable against evaporation [9] but its mean field growing like ln(𝑟²/𝑎²) up to 𝑟=∞ is unphysical.

King-like Lowered Isothermal Models
In the conventional literature, the above difficulties are usually circumvented by enclosing the system within a hypothetical box of some radius 𝑅_box [7], or by modifying the original distribution heuristically into non-Maxwellian form 𝑓^Low such as 𝑓^Low−𝜀∝exp−𝜀𝑘𝑇−1,∝exp+𝜀𝑘𝑇𝑘𝑇−1, (9) and so forth, holding in the range 𝜀<0 and vanishing elsewhere [4, 8, 16–19]. In particular, King [4] and Wilson [19] appealed to the tidal force field of the galaxy for physically setting its outer boundary and assumed the velocity distribution of stars to be cut off at the local escape velocity. Lowered isothermal prescriptions such as (9) are often employed by astronomers to fit data.

Physical Motivation for 𝑓^Low
If one takes a stellar cluster in an original Liouville collisionless state then the cluster will start evolving in space-time through trajectory mixing and stellar encounters which are most frequent in the core region. Mathematically, the complicated dynamics of such a nonequilibrium system is governed by the coupled Fokker-Planck and Poisson equations [17]. Physically, this evolution will involve momentum/mass/heat flow, tide generation, and entropy production. At equilibrium, the macroscopic flows will stop, tides will stabilize, and the entropy would become maximum. Naturally, von Hoerner [20] and King [14] realized that a finite boundary to the star cluster is set up by the tidal force of the galaxy, that is, the cutoff tail in the essentially classical stellar systems can be ascribed to the physical outcome of the boundary conditions and/or constraints (independent of the Plank constant).

3. Our Most Probable Distribution (MPD)

3.1. Preliminaries

We adopt the view that the above-mentioned King models can be refined further by utilizing the following facts. (i) At equilibrium, the entropy of a multiparticle thermal system should become a (local) maximum. Of course, the Boltzmann-Shannon definition of 𝑆^𝐵 in (1) will not serve the purpose due to the difficulties of the Maxwellian; we shall show in (11) and (12) below that a more suitable candidate is the so-called combinatorial entropy S that counts the number of microstates in energy cells corresponding to specified total particle number N and total energy E. (ii) The resulting most probable distribution f should develop a tail which is automatically truncated at a finite energy 𝜀_𝐾. This is because a star moving in the mean potential field 𝑊(𝑟) will have a farthest turning at distance 𝑟_𝐾 satisfying 𝑊(𝑟_𝐾)=𝜀_𝐾, where 𝑟_𝐾 may now be identified with the classical King radius of the galaxy. (iii) By Bohr’s correspondence principle, classical motion is the limiting case of quantum motion in states of very large quantumnumbers. The cutoff quantum number K and cutoff energy 𝜀_𝐾 should be determinable from the variational constraint equations of our MPD theory provided that h is brought into the picture explicitly. (iv) Our MPD solution for f should be able to provide a theoretical justification (or better characterization) of the lowered isothermal Maxwellian models (9). Now we shall demonstrate how such a task is accomplished in practice.

Gibbs Combinatorial Entropy
We follow the basic theme of Huang [1, page 182], and a preliminary investigation by Menon and Agrawal [12] as well as by Menon et al. [13]. The single-particle spectrum is divided into J cells into which the particles are distributed at random such that the 𝑗th cell has central energy 𝜀_𝑗, width Δ_𝑗, degeneracy 𝑔_𝑗, occupation number 𝜈_𝑗, and occupation probability per state 𝑓_𝑗=𝜈_𝑗/𝑔_𝑗 defined by Appendix A. Next, treating the particle in the 𝑗th cell as indistinguishable, a Gibbs combinatorial entropy functional S is constructed via 𝑆=𝑘ln^𝐽_𝑗=1𝑔^𝜈𝑗_𝑗𝜈_𝑗!. (10)

Gamma Function Form
We deliberately rewrite (10) in the equivalent form 𝑆=𝑘^𝐽_𝑗=1𝜈_𝑗ln𝑔_𝑗−lnΓ(𝜈_𝑗+1). (11) The replacement of factorials by gammas has several algebraic advantages. (i) The equality 𝜈_𝑗!=Γ(𝜈_𝑗+1) is exact at the integer values 𝜈_𝑗=0,1,2,3,…,∞. (ii) The asymptotic behaviour, namely, (2𝜋)^1/2𝜈^𝜈+1/2𝑒^−𝜈 of both 𝜈! and Γ(𝜈+1) are the same as 𝜈→∞. (iii) Hence, by a theorem due to Carlson [21, 22], Γ(𝜈+1) provides the most economical, essentially unique continuation of 𝜈! to all continuous values throughout the range 0≤𝜈≤∞. (iv) While setting up the variational conditions, later we shall need to replace the derivative 𝜕lnΓ(𝜈_𝑗+1)/𝜕𝜈_𝑗 evaluated at integer values by the digamma function 𝜓(𝜈+1) [23] computed at general continuous values. This problem of integer programming is handled in Appendix B by using an efficient finite-difference package for all natural numbers up to 4. (v) Appendix B also shows that the numerical differentiation of lnΓ(𝜈+1) can be readily done even at small values𝜈=0.2,0.4, and so forth, giving results in good agreement with 𝜓(𝜈+1).

Exact Variational Conditions
Next, we consider the following objective functional to be maximized: 𝑆^∗=𝑆+𝑘𝛼_𝑗𝜈_𝑗−𝑁−𝑘𝛽_𝑗𝜈_𝑗𝜀_𝑗−𝐸, (12) where 𝛼 and 𝛽 are unknown Lagrange multipliers. Equating to zero, the partial derivatives 𝜕𝑆^∗/𝜕𝜈_𝑗,𝜕𝑆^∗/𝜕𝛼, and 𝜕𝑆^∗/𝜕𝛽 lead to the following set of exact variational conditions still using the discrete index 𝑗 [12]: 𝜓(𝜈_𝑗+1)=ln𝑔_𝑗+𝛼−𝛽𝜀_𝑗∶1≤𝑗≤𝐽,_𝑗𝜈_𝑗=𝑁;_𝑗𝜈_𝑗𝜀_𝑗=𝐸. (13a)(13b)

Comments
Without making any assumption concerningthe largeness or smallness of 𝜈_𝑗, we can rewrite (13a) in the compact form𝜓(𝜈_𝑗+1)=_exactln𝜈^𝐵_𝑗;𝜈^𝐵_𝑗≡_def𝑔_𝑗𝑒^{𝛼−𝛽𝜀𝑗}, (14a) where the symbol 𝜈^𝐵_𝑗 was already encountered earlier in (2). In principle, (14a) can be solved for the desired cell occupation numbers 𝜈_𝑗 in terms of 𝛼,𝛽,𝜀_𝑗. Thereafter, the Lagrange multipliers 𝛼 and 𝛽 can be determined from the constraints (13b). Equivalently, the chemical potential 𝜇 and thermodynamic temperature T may be introduced via 𝜇𝑘𝛼=𝑇=−𝜕𝑆^∗1𝜕𝑁;𝑘𝛽=𝑇≡𝜕𝑆^∗.𝜕𝐸 (14b) Finally, if the total number J of levels is very large, we are permitted to take the quasilimit (A.4) leading to a continuous distribution for 𝜈 versus 𝜀 by dropping the index j and converting sums into integrals. Let us derive several interesting properties of our most probable distribution (MPD) defined by (13a), (13b) and (14a), (14b) with the suffix j omitted.

Location of The Peak
Differentiating (14a) with respect to 𝜀, we get 𝜕𝜓(𝜈+1)⋅𝜕𝜈𝜕𝜈=1𝜕𝜀𝜈^𝐵⋅𝜕𝜈^𝐵.𝜕𝜀 (15) Clearly, the MPD occupancy 𝜈 and MB occupancy 𝜈^𝐵 are both peaked at a common energy 𝜀_𝑝 which satisfies𝜈̇(𝜀_𝑝)=0;𝜈̇^𝐵(𝜀_𝑝̇)=0,𝑔(𝜀_𝑝)𝑔(𝜀_𝑝)−𝛽=0, (16) where the dot stands for derivative with respect to 𝜀. Typical algebraic estimates of 𝜀_𝑝 for the soluble potential models will be reported later in (28).

Large 𝜈 Region
In the so-called head region of the continuous distribution, the cells have large occupancy 𝜈≫1 so that the Stirling’s approximation 𝜓(𝜈+1)≈ln𝜈+𝑂(1/𝜈) holds in the fundamental equation (14a). Hence the MB solution is roughly retrieved, namely,𝜈≈𝜈^𝐵≫1;𝑓≈𝑒^{𝛼−𝛽𝜀}:Large𝜈, (17)but it must be violated in the cells where the occupancy becomes comparable to, or less than, unity.

The Tail Region
On the other extreme lies the tail region of the continuous distribution where the cell occupancy becomes small, that is, 𝜈≪1. Then the digamma function possesses a Taylor expansion 𝜓(𝜈+1)≈_𝜈≪1𝜓(1)+𝜁(2)𝜈+𝑂(𝜈²):Tail, (18) where 𝜓(1)=−0.577 is the negative of Euler’s constant and 𝜁(2)=𝜋²/6 is a Riemann zeta value. Substitution of the expansion (18) into the fundamental equation (13a) leads to the following three surprising yet important observations.
(i)The tail of the distribution intersects the energy axis at a cutoff point 𝜀_𝐾 such that𝜈_𝐾=0;𝜓(1)=ln𝜈^𝐵_𝐾≡ln𝑔_𝐾+𝛼−𝛽𝜀_𝐾, (19) where the suffix K refers to energy 𝜀_𝐾.(ii) The said intersection happens linearly because, in its neighbourhood, the occupancy 𝜈≈(𝜀_𝐾−𝜀)𝑔̇𝜁(2)𝛽−^𝐾𝑔_𝐾:Tail, (20) where 𝑔̇_𝐾 stands for 𝜕𝑔/𝜕𝜀 evaluated at 𝜀_𝐾.(iii)Extension of the graph of 𝜈 versus 𝜀 beyond the cutoff point is not allowed because that would tend to make 𝜈 negative in (13a), that is,𝜈_𝑗≡0;𝑓(𝜀)≡0:𝑗≥𝐾,𝜀≥𝜀_𝐾, (21)implying that the original occupied spectrum (A.2) has shrunk below J or 𝜀_𝐽 due to strict entropy maximization under stable equilibrium. (The possibility 𝐾>𝐽,𝜀_𝐾>𝜀_𝐽 would correspond to unstable equilibrium, that is, continuous evaporation of the system.) Schematic plots of 𝑓^𝐵 and 𝑓 versus 𝜀 are shown in Figure 1. Typical algebraic estimates of the cutoff energy 𝜀_𝐾 and quantum number K for the solvable models will be reported later in (28).

Figure 1 

Schematic plots (not to scale) of the distribution functions in the MB approach (cf. (2)) and MPD theory (cf. (23)). (Plot of the occupancy function would have shown a peak in both cases.)

Compact Solution for 𝜈(𝜀)
Our equation (14a) is a transcendental equation in 𝜈 and its precise analytical solution in closed form is not known. Fortunately, there exists an ansatz 𝜓(𝜈+1)=_ansatzln{𝜈+𝜃};𝜃≡𝑒^𝜓(1)=0.562, (22) which works excellently throughout 0≤𝜈≤∞ as shown graphically in Figure 2. Combining (14a) and (22), we obtain a very compact, quite accurate, MPD solution valid in all energy cells of relevance as𝜈=𝜈^𝐵−𝜃;𝑓=𝑓^𝐵−𝜃𝑔:𝜀₁≤𝜀≤𝜀_𝐾. (23) It is interesting to note that if 𝑔were replaced by a constant in (23), our MPD solution 𝑓 would agree with the first line of (9) implying a sort of justification for the lowered isothermal Maxwellian models. Actually, our numerically accurate solution (23) should be regarded as a better characterization since the degeneracy function 𝑔(𝜀) is strongly energy-dependent.

Figure 2 

The exact digamma function 𝜓(𝜈+1) and three numerical approximations to thesame used in the context of (22), (17), and (18). The values of three important constants are 𝜓(1)=−0.577,𝜃=𝑒𝜓(1)=0.562, and 𝜍(2)=𝜋2/6=1.646.

Compact Number Condition
Combining the number constraint (13b) with the general solution (23), we can define an effective number 𝑁^eff through 𝑁^eff≡𝑁+𝐾𝜃=^𝐾_𝑗=1(𝜈_𝑗+𝜃)=^𝐾_𝑗=1𝜈^𝐵_𝑗=^𝜀𝐾_𝜀0𝜕𝜀Δ⋅𝑔𝑒^{𝛼−𝛽𝜀}=𝑒^𝛼⋅𝑁^∗. (24a) Here we have employed the quasicontinuum limit (A.4) and introduced 𝑁^∗≡^𝜀𝐾_𝜀0𝜕𝜀𝑤(𝜀)𝑒^−𝛽𝜀;𝑒^𝛼=𝑁^eff𝑁^∗. (24b)This gives a formal expression for the Lagrange multiplier (or reduced chemical potential) 𝛼 provided that the underlying mean field W or its reduced degeneracy function 𝑤=𝑔/Δ is known.

Compact Cutoff Condition
Lastly, we convert the cutoff criterion (19) into 𝜃≡𝑒^𝜓(1)=𝜈^𝐵_𝐾=𝑔_𝐾𝑒^{𝛼−𝛽𝜀𝐾}. (25a) Eliminating 𝑒^𝛼 with the help of (24b), we find 𝑒^𝛽𝜀𝐾=𝑁^eff𝑔_𝐾𝑁^∗𝜃 (25b)which yields a formal expression for the cutoff energy 𝜀_𝐾=𝜀(𝐾) whose functional dependence can be inverted to specify also the number K of levels. The sharply cutoff tail of (21) will play a crucial role in the cosmological application to be discussed later in Section 4.

Illustration for the (Truncated) Oscillator Well
The above methodology may be illustrated in the case of the truncated harmonic oscillator potential listed in Table 1: 𝑊(𝑟)=−|𝑊₀|𝑟1−²𝑅²⋅step(𝑅−𝑟). (26) Before going ahead with the algebra, the following important remarks should be kept in mind. (a) If the well was untruncated, that is, the step function in (26) was absent then all particles would remain truly confined, the Boltzmann mass density (3) would vanish asymptotically, and the MB distribution would not be problematical. (b) However, if the well is truncated by the use of the step function in (26), then the potential vanishes for 𝑟>𝑅, particles can be ejected into the continuum, the MB distribution becomes problematical, and the MPD philosophy becomes very useful. (c) Near the origin the oscillator potential is rather flat, that is, smoothly varying so that it can approximately mimic the realistic mean field in the core region of astronomical galaxies. In sharp contrast, the truncated-linear and Coulomb-like potentials of Table 1 cannot do so since these vary rather quickly as 𝑟→0. (d) As they stand, the depth 𝑊₀ and range 𝑅 are only illustrative parameters introduced in (26). However, when we come to cosmological applications in Section 4 (especially the Poisson equation), it will be found that these parameters are directly related to the physical mass M and observed radius 𝑟_𝐾 of the galaxy. We are now ready to apply the MPD program to (26).

Table 1 

Properties of one-body semiclassical spectrum labeled by the integer j for four solvable potentials 𝑊(𝑟) in D dimensions. The well-depth 𝑊₀≡−|𝑊(0)| and the range R of the rectangular, linear, and oscillator wells are made finite using the unit step function. The Coulomb well is left unrestricted over all r. For other notations, see Appendix A. (We do not tabulate the rigid box potential model which has the upper end of the energy spectrum at 𝜀_𝐽=𝐸−(𝑁−1)𝑊₀≈+∞. Such a model is of little interest in cosmology In usual physical applications, one puts 𝐷=3.)

The Tilde Nomenclature
First, we read off the symbols 𝐽,𝜀,Δ,𝐵, and 𝑔 from the fourth column of Table 1. Next, for algebraic convenience, the following dimensionless quantities are defined along with the thermal de Broglie wavelength 𝜆: 𝐽≡𝐽|𝛽𝑊₀|^1/2=(𝑚|𝑊₀|/2)^1/2𝑅/ℏ(𝛽|𝑊₀|)^1/2=𝜋𝑅𝜆≫1,𝜆=𝜋ℏ2𝛽𝑚^1/2;𝜀≡𝛽(𝜀−𝑊₀𝑗)=𝐽,𝛼≡𝛼−𝛽𝑊₀,𝜀_𝐾≡𝛽(𝜀_𝐾−𝑊₀𝐾)=𝐽,𝑗=𝐽𝜀,𝑔=(𝐷𝐵𝐽^𝐷−1)𝜀^𝐷−1,𝜈^𝐵=𝑔𝑒^{𝛼−𝛽𝜀}=𝑔𝑒^𝛼−𝜀,𝑁≡𝑒^𝛽𝑊0𝑁^∗≡^𝐾₀𝑑𝑗(𝐷𝐵𝑗^𝐷−1)𝑒^𝐽−𝑗/=𝐷𝐵𝐽^𝐷𝛾_𝐾,𝛾_𝐾≡𝛾(𝐷,𝜀_𝐾). (27) A few remarks are in order concerning these definitions. The Planck constant ℎ or thermal de Broglie wavelength 𝜆≡ℎ/(2𝑚𝑘𝑇)^1/2 has appeared in the value of the symbol 𝐽 and the inequality 𝜆/𝑅≪1 is essential for the validity of classical motion (cf. (A.8)). The function 𝜀 measures the single-particle energy from the ground level in terms of 𝑘𝑇. The symbols 𝛼 and 𝜀_𝐾 may be called the dimensionless chemical potential and dimensionless cutoff energy, respectively, whose fixation using MPD constraints is yet to be done. The integral 𝑁 will play a crucial role below with 𝛾 being the incomplete gamma function.

Use of MPD Conditions
Remembering the tilde quantities, we can readily evaluate the conditions (16), (24b), and (25b). This yields the peak location 𝜀_𝑝, peak height 𝜈^𝐵_𝑝, dimensionless chemical potential 𝛼, and dimensionless cutoff energy 𝜀_𝐾 through 𝜀_𝑝=𝐷−1∼1;𝜈^𝐵_𝑝=𝑔_𝑝𝑒^{𝛼−𝜀𝑝}∼𝑁^eff𝐽𝛾_𝐾,𝑒^𝛼∼𝑁^eff𝐽^𝐷𝛾_𝐾;𝑒^𝜀𝐾𝜀_𝐾^𝐷−1∼𝑁^eff𝐽𝛾_𝐾𝑁^eff≡𝑁+𝐾𝜃;𝐷>1, (28) where the wavy symbol ~ implies the order of magnitude, and the multiplicative factors of order unity have been suppressed. We still have to show that the formal equations (28) do admit valid, that is, self-consistent MPD solutions under suitable restrictions. For this purpose, we consider below two cases in which the parameter 𝜀_𝐾=𝐾/𝐽 has markedly different behaviours.

Case 1 (well depth large compared to 𝑘 times temperature). For the truncated oscillator potential (26), we recall the tilde notations (27) and impose the following inequalities: 𝜀_𝐾=𝐾𝐽≫1;|𝛽𝑊₀|=|𝑊₀|𝐽∼𝑅𝑘𝑇≫1,𝜆≫1;𝐽≪𝑁^eff≪𝐽^𝐷,𝐽=(𝛽|𝑊₀|)^1/2𝐽>𝐾≫𝐽≫1. (29) The physical meaning of these restrictions is as follows. The inequalities 𝐽≫1,𝐾≫1,𝐽≫1 guarantee the validity of classical dynamics in states of large quantum numbers, the condition 𝐽>𝐾 ensures that the 𝐾th level lies below the ionization threshold for stable MPD, the assumption 𝜀_𝐾≫1 in Case 1 implies that the actual cutoff energy 𝜀_𝐾 is several 𝑘𝑇 above the ground level, |𝛽𝑊₀|≫1 means that the well depth is large compared to 𝑘𝑇, the condition 𝐽≪𝑁^eff implies that the effective number 𝑁^eff/𝐽 of particles grossly counted per cell is much more than unity, and 𝑁^eff≪𝐽^𝐷 means that the system is dilute or nondegenerate (because the packing fraction 𝑁^eff/𝐽^𝐷∼𝑁^eff𝜆^𝐷/𝑅^𝐷≪1, that is, the average number of particles contained inside a 𝐷-dimensional sphere of radius 𝜆 is small compared to unity). Then the incomplete gamma function 𝛾_𝐾→Γ(𝐷) and consistent handling of (28) leads to the estimates 𝜀_𝑝∼1;𝜈^𝐵_𝑝∼𝑁^eff𝐽𝑁≫1,𝛼∼ln^eff𝐽^𝐷;𝜀_𝐾𝑁∼ln^eff𝐽,𝐾∼𝐽𝜀_𝐾∼𝑅𝜆𝑁ln^eff𝜆𝑅. (30) The present case should apply to usual gases/liquids contained in the laboratory and we have independently verified that the functional forms of (29) and (30) are very rugged, that is, they hold for all the soluble models reported in Table 1.

Case 2 (well depth comparable to 𝑘 times temperature). Again we recall the tilde notations (27) and impose the orders of magnitude 𝜀_𝐾=𝐾𝐽𝑊∼1;⁰𝑘𝑇∼1. (31) Then the incomplete gamma function 𝛾_𝐾∼1 and (28) are found to admit the self-consistent estimates 𝜀_𝑝∼1;𝜈^𝐵_𝑝∼𝑁^eff𝐽𝑁∼1,_eff𝐽∼𝑅∼𝑁∼𝐾∼𝜆1≫1,𝛼∼ln𝐽^𝐷−1;𝜀_𝐾𝑁∼ln^eff𝐽∼1. (32) The physics of (31) and (32) is as follows. The statement |𝑊₀/𝑘𝑇|∼1 applies to gravitational systems obeying virialization, 𝑁∼𝐾 tells that the total number of particles is of the same order as the number of MPD cells, and 𝜈^𝐵_𝑝∼1 signifies that the cell occupancies have become comparable to unity with ℏ again playing a role through the symbol 𝐽. The present case should correspond to open astronomical systems and the ruggedness of the results (32) can be verified also for the other solvable models in Table 1.

4. Conceptual Application of MPD to Cosmology

We are now ready to resolve the conceptual difficulties of the MB distribution mentioned already in Section 2 by employing the MPD solution obtained in Section 3.

4.1. Entropy

The Boltzmann-Shannon entropy 𝑆^𝐵 of (1) is simple-minded, its maximization leads to the MB solution 𝜈^𝐵 in (2) with untruncated tail, and its generalization to quantum statistics is difficult. In sharp contrast, the combinatorial entropy S of (11) is sophisticated, its maximization leads to our MPD solution 𝜈 in (23) with a truncated tail, and its generalization to quantum statistics as straightforward.

4.2. Density

If the MPD information (21) is inserted back into the general expression (A.7) for the local mass density 𝜌, based on the transformation 𝜎=𝜀−𝑊(𝑟), we obtain 𝜌(𝑟)∝^𝜎𝐾₀𝜕𝜎𝜎^𝐷/2−1𝑓(𝜀);𝜎_𝐾≡𝜀_𝐾−𝑊(𝑟), (33) where 𝜌 surprisingly vanishes if W(r) equals 𝜀_𝐾. This is explained by remembering that since no particle in MPD is allowed to have an energy more than 𝜀_𝐾,there exists a largest classical turning point at 𝑟_𝐾 beyond which the density must become zero identically, that is, 𝑊(𝑟_𝐾)=𝜀_𝐾;𝜌(𝑟)≡0for𝑟≥𝑟_𝐾, (34) in sharp contrast to the MB density profiles (3)–(5). We can also find the rate at which 𝜌(𝑟) approaches zero as 𝑟 tends to 𝑟_𝐾. However, (20) has already told us that 𝑓(𝜀)∝𝜀_𝐾−𝜀∝𝜎_𝐾−𝜎 in the tail region. Hence, (33) yields the leading behaviour 𝜌(𝑟)∝_𝜎𝐾→0^𝜎𝐾₀𝜕𝜎𝜎^𝐷/2−1⋅𝜎∝_{𝑊(𝑟)→𝜀𝐾}{𝑊(𝑟_𝐾)−𝑊(𝑟)}^𝐷/2+1. (35) Since in a “good” MPD solution 𝜀_𝐾 and 𝑟_𝐾 are finite, our result (35) tells that the mass density obeys a (𝑟_𝐾−𝑟)^𝐷/2+1 law near the edge of the system.

4.3. Radius

Clearly, the distance 𝑟_𝐾 in (34) is the upper bound on the size of our galactic system and, for binding, we must have 𝑟_𝐾<𝑟_𝐽 with 𝑟_𝐽 being the turning point just before the continuum starts. (In the soluble models of Table 1, this 𝑟_𝐽 was called 𝑅). Since the density 𝜌 vanishes beyond 𝑟_𝐾 the MPD integral defining the average size⟨𝑟⟩ will also converge, that is, ⟨𝑟⟩<∞, (36) in sharp contrast to the MB mean radius (6).

4.4. Mass

By the same token, the MPD integral defining the total mass of the stellar system also exists, that is, 𝑀<∞, (37) in sharp contrast to the MB mass (7).

4.5. Nonevaporation

As is well known if an attractive mean field 𝑊(𝑟) vanishes asymptotically, then the energy 𝜀=0 is called the ionization threshold. Hence, our galaxy will be stable against evaporation if the MPD cutoff energy 𝜀_𝐾 happens to be negative at the given thermodynamic temperature 𝑇. Consequently, for 𝜀_𝐾<0, there is no net outgoing particle flux, that is, 𝜑=0, (38) in sharp contrast to the MB result (8).

Comment
Of course, a galaxy which is observed experimentally to evaporate is not in true equilibrium. Then simplifying restrictions like (31) may not hold, that is, the cutoff conditions (19) and (32) will admit a positive root for 𝜀_𝐾.

Poisson Equation Implications
So far in our treatment, the detailed algebraic form of the mean field W(r) was not required explicitly for self-gravitating systems. Actually this is a tough problem theoretically/numerically because one must solve the coupled equations for the distribution function f and mass density 𝜌 in accordance with the Poisson equation in 3 dimensions ∇²1𝑊≡𝑟²𝜕𝑟𝜕𝑟²𝜕𝑊𝜕𝑟=4𝜋𝐺𝑚𝜌;𝐷=3, (39) where 𝐺 is the gravitational constant. Our limited aim in the present paper will, however, be served by noting the following features.

4.6. Features

(a) Since the density 𝜌(𝑟) is sharply cutoff at 𝑟_𝐾, by Gauss theorem, the exact potential energy and force at exterior points become 𝑊(𝑟)=−𝐺𝑀𝑚𝑟;𝐹(𝑟)=−𝐺𝑀𝑚𝑟²:𝑟≥𝑟_𝐾. (40)

(b) At the edge 𝑟_𝐾 itself, the potential energy becomes equal to the cutoff energy, namely, 𝜀_𝐾=𝑊(𝑟_𝐾)=−𝐺𝑀𝑚𝑟_𝐾. (41)

(c) In the interior region, the mean field may get smoothened so as to yield a finite depth 𝑊₀≡𝑊(0)∼−𝑘𝑇 (42) by virtue of the gravitational virial theorem.

(d) At interior points, the exact profile of the mean field is not known a priori since it has to be, in general, computed numerically by solving (39). However, for the purposes of illustration, we can represent it by an oscillator form 𝑊(𝑟)=𝑊₀(1−𝑟²/𝑅²) if 𝑟≤𝑟_𝐾 with unknown phenomenological constants 𝑊₀ and R. The corresponding interior potential energy and force at the system edge 𝑟_𝐾 then become 𝑊(𝑟_𝐾)=𝑊₀𝑟1−^2𝐾𝑅²;𝐹(𝑟_𝐾)=2𝑊₀𝑟_𝐾𝑅². (43) Matching these to the exterior values given by (40) at 𝑟_𝐾, we identify 𝑊₀=3𝑊(𝑟_𝐾)2;𝑅=√3𝑟_𝐾. (44) Thus, 𝑊₀ is deeper than 𝑊(𝑟_𝐾) and R is larger than 𝑟_𝐾 (although orders of magnitude are the same).

Suggested Procedure for Cosmologists
Suppose a practical astronomical observation has been made on a cluster of 𝑁∼10⁵ stars. For utilizing our MPD theory with respect to his collected data, the cosmologist should proceed through the following steps.

Step1 (characterization parameters). From the observed size 𝑟_𝐾 and the known mass M of the cluster, the MPD cutoff energy 𝜀_𝐾 is immediately given by (41) as 𝜀_𝐾=𝑊(𝑟_𝐾)=−𝐺𝑀𝑚/𝑟_𝐾. Next, the oscillator well-depth 𝑊₀ and the range parameter R for motion inside the cluster are set up from (44) as 𝑊₀=3𝑊(𝑟_𝐾)/2 and 𝑅=√3𝑟_𝐾. Next, according to (32) applicable to cosmology, the MPD parameters have the rough orders of magnitude 𝐽∼𝑅𝐾∼𝑁∼𝜆∼(𝑚|𝑊(𝑟_𝐾)|)^1/2𝑟_𝐾ℏ∼𝐺𝑀𝑚²𝑟_𝐾ℏ^21/2 (45) upon using the value of 𝐽 given by the first line of (27) under virialization. Next, the cosmologist may treat (45) as providing a new mass versus radius relationship 𝑀/𝑟_𝐾∼𝐺𝑚⁴/ℏ² for nondegenerate clusters whose experimental status is, however, not yet studied. Finally, for an accurate interlink among all MPD parameters, the astronomer may like to solve the transcendental equations (27) numerically.

Step2 (MPD density near the edge). Next, the astronomer may look at (35) which gives the leading behaviour of the stellar number density near the cluster’s boundary: 𝑛(𝑟)∝_{𝑟→𝑟𝐾}{𝑊(𝑟_𝐾)−𝑊(𝑟)}^3/2+1∝_MPD1𝑟−1𝑟_𝐾^5/2 (46) since the mean field W(r) becomes Newtonian near the periphery. This can be cast into more convenient form by defining the variable 𝑥=1/𝑟, choosing a normalization point 𝑥₁=1/𝑟₁, and working with the modified function 𝑛^2/5≡_def𝑛(𝑟)𝑛(𝑟₁)^2/5=𝑥−𝑥_𝐾𝑥₁−𝑥_𝐾 (47) which becomes unity at 𝑥₁ but vanishes at 𝑥_𝐾=1/𝑟_𝐾. To test the validity of (47), the astronomer may, for example, concentrate on the star counts made on photographs of the cluster M 15 (see Figure 2 of King [14]) taken with the 48-inch Schimdt camera in the Palomer observatory. The results of 𝑛^2/5 are plotted in Figure 3. Clearly, there is quite good agreement between experimental observation and MPD prediction, although a slight curvature in the data trend may imply the presence of additional weak nonlinear terms on the RHS of (47).

Figure 3 

Normalized stellar number density raised to 2/5 power, that is, 𝑛2/5≡{𝑛(𝑟)/𝑛(𝑟1)}2/5 near the edge of the cluster M 15. The experimental data points are adapted from King [14, Figure 2]. The theoretical MPD prediction is computed from (47) with 𝑥1=0.202,𝑥𝐾=0.008.

Step3 (comparison with King density). Next, it is worthwhile to consider the function 𝑛^1/2 and expand its MPD expression (47) around the matching point 𝑥₁ binomially in the form 𝑛^1/2=_MPD{1+𝛿}^5/4=_𝑥→𝑥151+4𝛿+𝑂(𝛿²);𝛿≡𝑥−𝑥₁𝑥₁−𝑥_𝐾. (48) Dropping the 𝛿² term, the cosmologist retrieves the famous formula proposed empirically by King, namely, 𝑛^1/2=_King𝑥−𝑥_𝑡𝑥₁−𝑥_𝑡;𝑥_𝑡=4𝑥_𝐾+𝑥₁5, (49) whose square gives the King’s profile [14, equation (2)] near the cluster’s periphery as 𝑛(𝑟)∝_King1𝑟−1𝑟_𝑡². (50) It is well known that the phenomenological proposal (49) has been extensively used in the past by astronomers. For example, in context of M15 cluster Figure 4 shows the plot of 𝑛^1/2 near the cluster’s boundary. Clearly, the agreement between experimental observation and King’s parametrization is good, ignoring the slight curvature in the data trend. Incidentally, the qualities of fit seen in Figures 3 and 4 are quite comparable implying that, with the present accuracy of measuring 𝑛, it is not possible to say whether MPD formula (47) or King recipe (49) is superior.

Figure 4 

Normalized stellar number density raised to 1/2 power, that is, 𝑛1/2 near the boundary of the M 15 cluster. The experimental data points are read from King [14, Figure 2]. His model fit is given by (49) with 𝑥1=0.202,𝑥𝑡=0.047.

Step 4 (complete density profile). Finally, the astronomer may like to have an expression for the number density 𝑛(𝑟) valid throughout the range 0≤𝑟≤𝑟_𝐾. In principle, our MPD distribution function f given by (23) yields the formal expression 𝑛(𝑟)=_MPD𝜕^𝐷→𝑝ℎ^𝐷𝑓^𝐵−𝜃𝑔 (51) with the mean potential 𝑊(𝑟) being approximately harmonic oscillator in the interior and Newtonian near the edge. Unfortunately, analytical evaluation of the phase space integral (51) is somewhat tedious and will be dealt with in a future communication. However, the cosmologist should note that the integral (51) is the algebraic difference of two terms which is very satisfying because the empirical full density profile written by King [14, equation (14)] also contains a difference of two terms.

Step 5 (velocity distribution of stars). It is a standard astronomical practice to measure the local radial velocity distributions (along with other properties) of stars in a globular cluster, for example, see the extensive photometric study made by van Loon et al. [15] on the post-main-sequence stars in 𝜔 Centauri (NGC 5139). The cosmologist may ask how well our MPD distribution function f given by (23) fits the observed data. Unfortunately, a straightforward answer to this question is difficult because exact values of the unknown parameters 𝛼 and K must be obtained numerically from the transcendental conditions (24b) and (25b). We plan to accomplish this task in a future communication.

5. Concluding Remarks

The main results of the present work appear in the abstract along with Sections 2–4 and are often emphasized by italics. It is hoped that astronomers will benefit from the algebraic properties of MPD derived in Section 3, its cosmological implications mentioned in Section 4, numerical plots of number density profiles in Figures 3 and 4, and a pointwise comparison between the King model and MPD philosophy made in Table 2. Clearly, both types of theories can be applied to cosmology although our f may be regarded as providing a better characterization from the conceptual viewpoint.

Table 2 

Salient features of King-like lowered isothermal models 𝑓^Low (cf. (9)) and our MPD distribution f (cf. (23)).

The essence of our cosmological discussion in Section 4 is the following. Suppose that an astronomer makes observations on a (quantum mechanically nondegenerate) cluster having N stars, total mass M, and radius 𝑟_𝐾. Then, its MPD solution will be characterized by the cutoff energy 𝜀_𝐾=−𝐺𝑀𝑚/𝑟_𝐾 and cutoff quantum number 𝐾∼𝑁. Before ending the paper, we mention below briefly several important points not discussed explicitly in the earlier sections.

(i) In the mean field description of a multiparticle system, fluctuations arising from short-range pair correlations are usually ignored. The effect of fluctuations is likely to be stronger on the MB solution 𝑓^𝐵 whose tail extends to ∞ in (2). Such effect is likely to be weak on the MPD solution f whose tail gets truncated at 𝜀_𝐾 in (23).

(ii) One may argue that a sharp radius is also known to arise in the method of self-consistent Hartree fields applied to gravitational systems [17]. We stress, however, that the Hartree method is done through a numerical algorithm because the coupled equation for the mean field and distribution function must be solved iteratively on a computer. Therefore, our analytical maximum-entropy treatment of Section 3 still retains its novelty.

(iii) One may also argue that it is not meaningful to demand thermodynamic equilibration in the peripheral region of the galaxy because, due to low densities, relaxation may remain incomplete there. However, it must be kept in mind that since gravitational forces are of long range, the mechanism of collisionless relaxation [9] still operates. Therefore, our assumption of equilibration even in the tail region may remain justified.

(iv) Next, mention must be made of some recent investigations [10, 11] carried out on the question of gravitational galactic clustering, their virialization, and peculiar velocity distribution superposed over the local Hubble flow. These authors start from the N body cosmological canonical partition function 𝑍_𝑁 in a box of large volume 𝑉, perform the individual momentum integrals at the outset over the infinite domain −∞≤𝑝≤∞, write the entropy S as the logarithm of a Gibbs integral over the density of states, and minimize the Helmholtz free energy 𝐸−𝑇𝑆 with respect to the internal energy E. Of course, these investigations are very different from our work because we do not need an enclosing box, momentum integrations over infinite domain are never performed, the entropy functional is combinatorial, and maximization is done with respect to the cell occupation numbers.

(v) Next, suppose that one considers a time span long compared to the two-body relaxation time in a globular stellar cluster. One may argue that a star having energy 𝜀_𝑒=−0 (i.e., arbitrarily close to zero but still negative) will go far away and yet come back. Since the corresponding turning point 𝑟_𝑒 may be arbitrarily big, one expects a very small (but not zero) possibility of the star’s existence even at a very large radius. This logic apparently contradicts the MPD result (34) which had claimed that there is no density outside a finite distance 𝑟_𝐾.

Actually, the above logic has the following very subtle fault. While doing pure dynamics, it is enough to find trajectories and their turning points 𝑟_𝑒; but while doing statistical mechanics, it is essential also to calculate the density profile 𝜌_𝑒(𝑟) and the related total mass 𝑀_𝑒. Now, in direct analogy with (35) but with 𝜀_𝑒=−0, the density profile at large distance and its associated Poisson equation become (in 𝐷=3 dimensions) 𝜌_𝑒(𝑟)∝_𝑟→∞{−𝑊_𝑒(𝑟)}^5/2∝∇²𝑊_𝑒(𝑟);∴𝑊_𝑒(𝑟)∝𝑟^−4/3. (52) This result is physically unacceptable because the gravitational potential due to a finite mass object must fall asymptotically like 1/𝑟. Hence a logic based on 𝜀_𝑒=−0 will not work. In sharp contrast, if the globular cluster has finite experimental mass 𝑀, then it can be easily described by our MPD solution (34) characterized by bounded 𝑟_𝐾 and finite 𝜀_𝐾<−0.

(vi) Finally, a cosmologist may argue that since real gravitating systems have a mass spectrum of stars, the assumption of particles with the same mass m in MPD may not be justified. We wish to point that some workers have attempted to apply hydrodynamical equations to globular clusters employing a phase space density involving the continuous mass [24] as an extra variable. Some other workers have analyzed phenomenologically the mergence of clusters such as Praesepe [25] employing four mass bins. Although, in principle, a multicomponent combinatorial entropy will now replace (10), yet the corresponding variational conditions (13a) and (13b) will be hard to handle analytically because different chemical potentials and different cutoff energies may have to be assigned to various components present in the system. An easy approximation will be to still use the MPD formalism of Section 3 based on the single particle average mass ∑𝑚=^𝑖𝑁_𝑖𝑚_𝑖∑_𝑖𝑁_𝑖, (53) where the suffix 𝑖 runs over different species and there are 𝑁_𝑖 particles of the 𝑖th type. This prescription should be reasonable for those clusters where the mass dispersion is small (in units of the solar mass).

Appendices

A. One-Body Description Recapitulated

A.1. Preliminaries

This section will summarize our notations along with several known formulae dealing with the semiclassical single-particle spectrum/distribution without invoking entropy constraints. Some of these formulae will be used explicitly in Sections 3–5 of the text.

A.2. Assumptions and Notations

Consider the nonrelativistic localized motion of a particle in 𝐷 spatial dimensions under the influence of a smooth attractive central field. Classically, the symbols 𝑚,𝑟,𝑝,𝑊(𝑟),𝐹(𝑟)=−𝑊^′𝑝(𝑟),𝜀=²2𝑚+𝑊(𝑟), (A.1) respectively, denote the mass, distance, absolute momentum, potential energy, applied force, and mechanical energy of this particle. Quantum mechanically invokes the Planck constant ℎ≡2𝜋ℏ and solves the Schrödinger equation for determining the energy spectrum 𝜀₀<𝜀₁<𝜀₂<⋯<𝜀_𝐽, (A.2) where 𝜀₀ is the ground level and 𝜀_𝐽 the highest bound level supported. Of course, solution of the Schrödinger equation for the exact eigenvalues, eigenfunctions, and their degeneracy is generally tedious.

A.3. Sommerfeld Quantization

Perhaps the easiest semiclassical link between the descriptions (A.1) and (A.2) is provided by Sommerfeld’s criterion [26] which says that the phase integral or action variable over a complete oscillation should be an integer multiple of ℎ. Then a discrete level 𝜀_𝑗 in (A.2) corresponds to the classical turning point 𝑟_𝑗, local momentum variable 𝑝, principal quantum number 𝑗, and level spacing Δ_𝑗 given by 𝑊(𝑟_𝑗)=𝜀_𝑗;𝑝={2𝑚(𝜀_𝑗−𝑊(𝑟)}^1/2,2𝑗=∫𝜋ℏ^𝑟𝑗01𝜕𝑟,𝑝=0,1,2,…,Δ_𝑗=𝜕𝑗𝜕𝜀_𝑗=2𝑚∫𝜋ℏ^𝑟𝑗0𝜕𝑟𝑝. (A.3) Since the presence of zero point energy is of little consequence here, hence the more sophisticated WKB quantization [27] will not be needed for our purpose. Also, if the number J of supported levels is very large compared to unity, then the quasicontinuum limit can be taken by writing 𝐽≫1;0≤𝑗≤𝐽;𝜀=𝜀_𝑗,Δ=Δ_𝑗;∑^𝐽_𝑗=0=∫^𝜀𝐽_𝜀0𝜕𝜀Δ. (A.4)

Gibbs’ Prescription
Further information is obtained by imagining a spherical region of range R and remembering that several quantum states of different orbital angular momenta and magnetic projections may be nearly degenerate at a given energy level. Then, the D-dimensional solid angle Ω, total volume V of the region, useful geometrical factor A, Gibbs’ phase space element 𝜕Ξ, accumulated number 𝐼_𝑗 of quantum states below 𝜀_𝑗, local number 𝑤_𝑗 of states per unit energy interval, and the degeneracy 𝑔_𝑗 of the 𝑗th level itself are read off from Ω=2𝜋^𝐷/2Γ𝐷2;𝑉=Ω𝑅^𝐷𝐷,Ω𝐴=²ℎ^𝐷=2^𝐷−2ℏ^𝐷Γ²𝐷2⁻¹,𝜕𝜕Ξ=^𝐷→𝑥𝜕^𝐷→𝑝ℎ^𝐷,𝐼_𝑗=𝐼(𝜀)=𝜕Ξstep(𝜀−𝐻)=𝐴𝐷^−1𝑟𝑗₀𝜕𝑟𝑟^𝐷−1𝑝^𝐷,𝑤_𝑗=𝑤(𝜀)=𝜕𝐼=𝜕𝜀𝜕Ξ𝛿(𝜀−𝐻)=𝑚𝐴^𝑟𝑗₀𝜕𝑟𝑟^𝐷−1𝑝^𝐷−2,𝑔_𝑗=𝑔(𝜀)=𝜕𝐼𝜕𝑗=𝑤Δ. (A.5) Here 𝐻 is the single-particle Hamiltonian, the quasicontinuum limit (A.4) is understood, Γ is the gamma function, step the unit step function, and 𝛿 the delta function.

Solvable Potential Models
The methodology described by (A.3)–(A.5) is best illustrated in the case of 4 soluble models, namely, the rectangular, truncated linear, truncated harmonic oscillator, and Coulomb wells. The results are summarized in Table 1 and the following features are worth noticing.
(i) In the case of the rectangular, linear, and oscillator wells, the range R represents the distance beyond which the particle goes into the continuum. By the same token, the highest level J is fixed through the requirement that 𝜀_𝐽=0.(ii) For the Coulomb well, however, since the bound orbits can have any size, one sets 𝑅=∞. By the same token, the ionization threshold appears at 𝐽=∞,𝜀_𝐽=0.(iii) In every model of Table 1, the semiclassical energy ε increases monotonically with the principal quantum number j, but the trend of the level spacing Δ is not uniform.(iv) In every soluble model, the level degeneracy 𝑔_𝑗=𝐷𝐵𝑗^𝐷−1, where the geometrical factor 𝐵 is of order unity. Hence it is reasonable to expect that, for a more general attractive central field in 𝐷=3 dimensions, 𝑔∝𝑗² at least for large 𝑗.(v) Table 1 does not explicitly treat the infinite rectangular well, that is, rigid box in which particles of any momentum would remain confined. Then, the highest kinematically allowed level would have 𝜀_𝐽=𝐸−(𝑁−1)𝜀₀ using the many-body notation of (A.6) below. Of course, the rigid box model is irrelevant in cosmology.

Multiparticle, Statistical System
In the present paper, we shall not consider pure Bose/Fermi many-body systems where the strict quantum mechanical identity of all particles is crucial. Ours is the so-called Boltzmann system where the one-body spectrum is obtained from the Schrödinger equation/semiclassical quantization but strict identity among all particles is not imposed except within the same energy cell. The mean field 𝑊(𝑟) of (A.1) may be either externally applied or internally generated. Assuming spherical symmetry, independent-particle motion, and ignoring short-range pair correlations, we let the symbols 𝑁𝑁,𝐸,𝑛=𝑉1,𝑇,𝑘,𝛽=ℎ𝑘𝑇,𝜆=(2𝑚𝑘𝑇)^1/2=𝜋ℏ2𝛽𝑚^1/2, (A.6)respectively, denote the specified number of particles, total energy, global average number density, global mean thermodynamic temperature, Boltzmann constant, inverse temperature parameter, and thermal de Broglie wavelength. The one-body phase space may be imagined to be composed of the differential elements 𝜕Ξ (cf. (A.5)) or of J energy cells of successive widths Δ_𝑗 which are arranged in the sequence (A.2). Then a useful transformation 𝜎, single-particle energy 𝜀, one-body distribution function 𝑓_𝑗, cell occupancy 𝜈_𝑗, local number density 𝑛(𝑟), local mass density 𝜌(𝑟), total number N, and total energy E are read off from 𝑝𝜎=²𝑓2𝑚≥0;𝜀=𝜎+𝑊(𝑟),_𝑗=𝑓(𝜀);𝜈_𝑗=𝑔_𝑗𝑓_𝑗=𝑤𝑓Δ,𝑛(𝑟)=𝜌(𝑟)𝑚=𝜕^𝐷→𝑝𝑓ℎ^𝐷=1𝑚Γ(𝐷/2)2𝜋ℏ^{2𝐷/2𝜀𝐽0−𝑊(𝑟)}𝜕𝜎𝜎^𝐷/2−1𝑓,𝑁=^𝐽_𝜈=0𝜈_𝑗=𝜕Ξ𝑓=^𝜀𝐽_𝜀0𝜕𝜀𝑤𝑓,𝐸=^𝐽_𝑗=0𝜈_𝑗𝜀_𝑗=𝜕Ξ𝑓⋅𝜀=^𝜀𝑇_𝜀0𝜕𝜀𝑤𝑓⋅𝜀. (A.7) Two crucial comments are in order at this stage. (i) The functional form of 𝑓(𝜀) is left unspecified at the moment. (ii) Convincing justifications are still needed for retaining ℎ≡2𝜋ℏ in our mechanical as well as statistical expressions (A.3)–(A.7) especially when application to classical galaxies of enormous sizes is being envisaged.

Importance of Planck Constant
(a) By Bohr’s correspondence principle, the motion of a quantum Schrödinger/Sommerfeld particle tends to become classical in the states of large principal quantum numbers. In the notation of (A.4), this requires 𝐽≫1 where J of Table 1 contains ℏ explicitly. (b) Strict Bose/Fermi statistical systems tend to obey classical statistics at low density and high temperature if 𝑛𝜆^𝐷≪1 in the notation of (A.6). This requires that the linear size of the system be large compared to the thermal de Broglie wavelength, that is, 𝑅≫𝜆.Hence the ℏ-dependent dual inequalities𝐽≫1;𝑅≫𝜆 (A.8)tell very precisely when a multiparticle system can be called “classical.” Such a precision would be lacking if ℏ were dropped at the outset in cosmological applications. (c) While the Sommerfeld quantum number j in (A.3) is very suitable for labelling the distinct energy levels, the Gibbs degeneracy g (derived from the phase-space element 𝜕^𝐷→𝑥𝜕^𝐷→𝑝/ℎ^𝐷) in Table 1 is equally convenient to count the precise number of states in any cell. (d) The precise knowledge of a cutoff quantum number 𝐾 and energy 𝜀_𝐾 will be shown to be crucial to find the rigorous most probable distribution 𝑓 in Section 3 which job cannot be done in the cosmological context of Section 4 if ℏ is dropped at the outset (in the classical phase space element 𝜕^𝐷→𝑥𝜕^𝐷→𝑝).

B. Extension from Integer to Continuous 𝜈

In this appendix, we carefully examine the numerical justification of some algebraic manipulations done on the combinatorial entropy S of (10), (11).

Factorials versus Gammas
As is well known, 𝜈_𝑗! identically equals Γ(𝜈_𝑗+1) at all nonnegative integers𝜈_𝑗=0,1,2,3,4,… as seen from the second line of the following brief table. Its third line records the corresponding values of the natural logarithm lnΓ(𝜈_𝑗+1) to be used as the input in Table 3.

Table 3 

B.1. Numerical Differentiation

Next, we address the subtle question of computing 𝜓_num(𝜈)≡𝜕lnΓ(𝜈+1)𝜕𝜈_num:0≤𝜈≤4, (B.1) where the suffix “num” stands for “numerically” and the inequality 0≤𝜈≤4 implies that 𝜈 has become a continuous variable over a test range [0, 4]. This is a problem of integer programming and we tackle it by adopting the following procedure.

(i) First, a finite-difference table was prepared using the above-mentioned data on lnΓ(𝜈_𝑗+1). (ii) Next, at several chosen integral/fractional values of 𝜈, (B.1) was computed employing an efficient package based on Markoff’s version of Newton’s interpolation differentiated [23, page 883]. (iii) Finally, comparison was made with the standard values of the digamma function [23, pages 258, 267, 272] obtained from the “exact” definition 𝜓_exact(𝜈+1)≡𝜕lnΓ(𝜈+1)=1𝜕𝜈Γ(𝜈+1)^∞₀𝑑𝑡𝑡^𝜈𝑒^−𝑡ln𝑡. (B.2)

B.2. Results

The accompanying table shows that 𝜓_num of (B.1) and 𝜓_exact of (B.2) agree within 1% to 5% at the input integer values 𝜈=0,1,2,3,4. We also see that their mutual agreement is good at the small fractional values of 𝜈=0.2,0.4,0.6,0.8. Therefore, taking derivatives of lnΓ(1+𝜈) at all continuous values of 𝜈 (including 𝜈≪1) is mathematically justified in (13a) and (14a) (see Table 4).

Table 4 

Acknowledgment

The authors thank the Council of Scientific and Industrial Research (CSIR), New Delhi, India, for the financial support.