Journal of Modern Physics
Vol.10 No.07(2019), Article ID:92917,12 pages

Statistical Induced Dynamic of Self-Gravitating System

B. I. Lev

Bogolyubov Institute for Theoretical Physics, NAS of Ukraine, Metrolohichna, Ukraine

Copyright © 2019 by author(s) and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

Received: April 29, 2019; Accepted: June 3, 2019; Published: June 6, 2019


The self-gravitating systems a priory are non-equilibrium. The new approach, based on non-equilibrium statistical operator, which allows taking into account inhomogeneous distribution of particles and temperature, has been proposed. Such method employs a saddle point procedure to find dominant contributions to the partition function and permits to obtain all thermodynamic parameters of the system. Statistical induced dynamic and behavioral peculiarities of self-gravitating systems for different conditions have been predicted.


Nonequilibrium Statistical Operator, Self-Gravitating System, Inhomogeneous

1. Introduction

The statistical description of a system of interacting particles has a fundamental physical interest. In that context, the interacting system provides a model of fundamental interest for which ideas of statistical mechanics and thermodynamics can be tested and developed. For self-gravitating system the thermodynamical ensemble can be nonequivalent [1] [2] . In particular, since energy is not-additive, and cannot use the canonical ensemble to study the system with long-range interaction. Since equilibrium states are only local entropy maximum.

Two type approaches (statistical and thermodynamic) have been developed to determine the equilibrium states of the interacting system and describe possible phase transition. It is generally believed that mean field theory is exact for systems. Since in mean field theory any thermodynamical function depends on the thermodynamical variables only through the dimensionless combinations and the system is thermodynamically stable but the thermodynamical limit does exist [3] .

Formation of the spatially inhomogeneous distribution of interaction particles is a typical problem in condensed matter physics and requires non-conventional methods of the statistical description of the system tailored to gravitational interacting particles with regard for an arbitrary spatially inhomogeneous particles distribution. Statistical description must employ the procedure to find dominant contributions to the partition function and to avoid entropy divergences for infinite system volume. Formation of the inhomogeneous distribution of interaction particles requires a non-conventional method, which was proposal in [4] [5] [6] [7] [8] . This method is based on Hubbard-Stratonovich representation of statistical sum [9] which is now extended and applied to the self-gravitating system to find solution for particles distribution without using spatial box restrictions. It is important that this solution has no divergences in thermodynamic limits. For this goal using the saddle point approximation took into account the conservation number of particles in limiting space and provided to nonlinear equation. The partition function in the case of homogeneous distribution of the particles, and in case of inhomogeneous distribution was obtained in [5] [6] [10] , but in this approach described only condition in the formation of the possible inhomogeneous distribution of the interacting particles.

Systems with long-range interactions, such as self-gravitating system do not relax to the usual Boltzmann-Gibbs thermodynamic equilibrium, but become trapped in quasi-stationary states, the lifetime of which diverges with the number of particles. The instability threshold was predicted for spontaneous symmetry breaking for a class of d-dimensional systems [11] . Non-equilibrium stationary states of systems were described in article [12] where concluded that three-dimensional systems do not evolve to thermodynamic equilibrium but trap in non-equilibrium quasi-stationary states.

Mostly due to the fact that self-gravitating system exist in a states that far from equilibrium and the relaxation time to equilibrium state is very long. The homogeneous particle distribution in a self-gravitating system is not stable. The particles distribution in such a system is inhomogeneous in space, from the beginning. There is standardized behavior of the interacting system, it is described differently for different equilibrium ensembles. That approach does not seem to be very consistent, due to the fact that, the state equation has to come from the definition of partition function, but the definition of it is unknown for the space inhomogeneous systems [13] [14] . Therefore, there is a dilemma, either to take the postulates of equilibrium statistical mechanics to obtain only instability criteria or not to try to account for the space in homogeneity using different approach. Such inhomogeneous distribution of the particles, temperature and chemical potential can be accounted for in non-equilibrium statistical operator approach [15] , where possibility of local change of thermodynamic parameters is included. This system is non-equilibrium and inhomogeneous distribution of particle can motivate the inhomogeneous distribution of temperature and chemical potential and other thermodynamic parameter.

In this article are proposed the new approach, based on non-equilibrium statistic operator [15] , which are more suitable for the description of a self- gravitating systems. The equation of state and all thermodynamic characteristics needed are defined by the equations which contribute the most to the partition function. Therefore, there is no need to introduce additional hypothesis about temperature on density dependence. This dependence is an outcome of solving corresponding thermodynamic relations, which describe extreme of non-equilibrium partition function. The possible space inhomogeneous distributions of particles and temperature have been obtained for the simple cases. For the equilibrium case the well-known result [16] [17] of partition function has been reproduced. Has been shown, that used this approach can describe inhomogeneous distribution of particle and determine necessary parameter in the self -gravitating system. The main idea of this paper is to provide a possible detailed description of self -gravitating system based on the principles of non-equilibrium statistical mechanics and obtain possible distribution particle by fixed number of particle and energy of the system. After that we can determine the dynamic of self -gravitating system which induced the statistical principles and low increasing the entropy in the case existence the thermodynamic limit.

2. Nonequilibrium Statistical Operator

Phenomenological thermodynamic based on the conservation lows for average value of physical parameter as number of particles, energy and impulse. Statistical thermodynamic non-equilibrium system based too on conservation lows not the average value dynamic variables but in particular for this dynamic variable. It present local conservation lows for dynamic variables. For the determination thermodynamic function of non- equilibrium system are teed use the presentation of corresponding statistical ensembles which take into account the non-equilibrium states of this systems. The conception of Gibbs ensembles can brings to description non-equilibrium stationary states of system. In this case can determine non-equilibrium ensemble as totality of system which be contained in same stationary external action. This system has same character of contact with thermostat and possess all possible value macroscopic parameters which compatibility present conditions. In system, which are in same stationary external condition will are formed local equilibrium stationary distribution. If external condition will be depend from time that local equilibrium distribution are not stationary. For exactly determination local equilibrium ensemble must accordingly determine the distribution function or statistical operator of the system [15] , finally, can recall that the stable states on the series of equilibrium of classical self-gravitating particles are only metastable because they correspond to local maximal of entropy from which can determine behavior of the system.

If assume that non-equilibrium states of system can determine through inhomogeneous distribution energy H ( r ) and number of particles (density) n ( r ) the local equilibrium distribution function for classical system can write in the form [15] :

f l = Q l 1 exp { ( β ( r ) H ( r ) η ( r ) n ( r ) ) d r } (1)


Q l = D Γ exp { ( β ( r ) H ( r ) η ( r ) n ( r ) ) d r } (2)

present determination of statistical operator local equilibrium distribution. The integration in present formula must take over all phase space of system. Must note, that in the case local equilibrium distribution Lagrange multipliers β ( r ) and η ( r ) are function of spatial point. The microscopic density of particles can present in standard form:

n ( r ) = i δ ( r r i ) (3)

The introduction local equilibrium distribution are possible if relaxation time in all system is more as relaxation time on local macroscopically area as part of this system.

After determination of non-equilibrium statistical operator can obtain all thermodynamic parameter non-equilibrium system. For this can determine thermodynamic relation for inhomogeneous systems. The variation of statistical operator by Lagrange multipliers can write necessary thermodynamic relation in the form [15] : δ ln Q l δ β ( r ) = H ( r ) l and δ ln Q l δ η ( r ) = n ( r ) l . This relation is natural general prolongation well-known relation which take place in the case equilibrium systems, on the case inhomogeneous system. The conservation number of particles and energy in system can present in form natural relations n ( r ) d r = N and H ( r ) d r = E

For further statistical description of non-equilibrium system is necessary determine Hamiltonian of system. In general case Hamiltonian of system interacting particle can present as:

H = i p i 2 2 m 1 2 i , j W ( r i r j ) (4)

where W ( r i r j ) determine the attractive interaction energy. Further will be use only density of energy which for self-gravitating system can write in the following form:

H ( r ) = p 2 ( r ) 2 m n ( r ) 1 2 W ( r , r ) n ( r ) n ( r ) d r (5)

This density energy of a system is possible to use if we smash to equal bits all space with equal mass and consider moving in phase space not compressed gravitational fluid. In our case interacting system can write are non-equilibrium statistical operator in the form:

Q l = D Γ exp { ( β ( r ) p 2 ( r ) 2 m η ( r ) ) n ( r ) d r + 1 2 β ( r ) W ( r , r ) n ( r ) n ( r ) } (6)

The integration over phase space can present as D Γ = 1 ( 2 π ) 3 i d r i d p i .

In order to perform a formal integration in second part of this presentation, additional field variables can be introduced making use of the theory of Gaussian integrals [8] [9] :

exp { ν 2 2 β ( r ) ω ( r , r ) n ( r ) n ( r ) d r d r } = D σ exp { ν 2 2 ( β ( r ) ω ( r ) ) 1 σ ( r ) σ ( r ) d r d r ν β ( r ) σ ( r ) n ( r ) d r } (7)

where D σ = s d σ s det 2 π β ω ( r , r ) and ω 1 ( r , r ) is the inverse operator which satisfies the condition ω 1 ( r , r ) ω ( r , r ) = δ ( r r ) , that means: the interaction energy presented the Green function for this operator and ν 2 = ± 1 depending on the sign of the interaction or the potential energy. After present manipulation the field variable σ ( r ) contains the same information as original distribution function, i.e. all information about possible spatial states of the systems.

After this manipulation the statistical operator can rewrite in the form:

Q l = D Γ D ψ exp { ( β ( r ) p 2 ( r ) 2 m η ( r ) β ( r ) ψ ( r ) ) n ( r ) d r } Q i n t (8)

where part which cam from interaction

Q i n t = exp { 1 2 ( β ( r ) W ( r , r ) ) 1 ψ ( r ) ψ ( r ) d r d r } (9)

In this general functional integral can be provide the integration on phase space. If use the definition of density we can rewrite the non-equilibrium statistical operator as:

Q l = D ψ 1 ( 2 π ) 3 N ! i d r i d p i ξ ( r i ) exp { ( β ( r i ) p i 2 2 m β ( r i ) ψ ( r i ) ) } Q i n t (10)

where introduce the new variable ξ ( r ) exp η ( r ) which can interpreter as chemical activity. Now can make integration over impulse. Real part of non-equilibrium statistical operator take the form:

Q l = D φ D ψ Q i n t 1 N ! i d r i ξ ( r i ) ( 2 π m 3 β ( r i ) ) 3 2 exp ( β ( r i ) ψ ( r i ) ) (11)

Q l = D φ D ψ Q i n t N 1 N ! ( d r ξ ( r ) ( 2 π m 3 β ( r ) ) 3 2 exp ( β ( r ) ψ ( r ) ) ) N (12)

After that the non-equilibrium statistical operator take the simple form:

Q l = D φ D ψ Q i n t exp { [ ξ ( r ) ( 2 π m 3 β ( r ) ) 3 2 exp β ( r ) ψ ( r ) ] d r } (13)

In definition Q i n t are the reverse operator of interaction energy. For further description we must determine this inverse operator. For general case of long range interaction as Coulomb like or Newtonian gravitational interaction in continuum limit should be treated in the operator sense, i.e.

W 1 ( r , r ) = 1 4 π G m 2 Δ r δ ( r r ) (14)

where Δ r is Laplace operator in real space. Gravitation interaction energy in three dimensional case can write in well-known form

W ( r , r ) = G m 2 | r r | (15)

G is gravitational constant and m is mass of separate area. For long range interaction between particle we can rewrite the non-equilibrium statistical operator in the form:

Q l = D ψ exp { [ 1 2 r ψ [ ( ψ ) 2 ] ξ Λ 3 exp β ψ ] d r } (16)

where all function β , φ , ψ are dependence from spatial point. In this presentation was using the definition thermal de-Broglie wavelength and definition of interaction length as

Λ 1 ( r ) = ( 2 m 2 β ( r ) ) , r ψ ( r ) = 4 π G m 2 β ( r ) , (17)

In general case the non-equilibrium statistical operator can rewrite in the form:

Q l = D ψ exp { S ( ψ ( r ) , ξ ( r ) , β ( r ) ) } (18)

where effective non-equilibrium “local entropy” take the form:

S = [ 1 2 r ψ ( ψ ( r ) ) 2 ξ ( r ) Λ 3 exp β ( r ) ψ ( r ) ] d r (19)

The statistical operator allows obtain use the of efficient methods developed in the quantum field theory without imposing additional restrictions of integration over field variables or the perturbation theory. The functional S ( φ ( r ) , ξ ( r ) , β ( r ) ) depends on the distribution of the field variables φ ( r ) , the chemical activity ξ ( r ) and inverse temperature β ( r ) . The saddle point method can now be further employed to find the asymptotic value of the statistical operator Q l for increasing number of particles N to . The dominant contribution is given by the states which satisfy the extreme condition for the functional. It’s easy to see that saddle point equation present thermodynamic relation and it can write in the other form as: equation for field variable δ S δ ψ ( r ) = 0 for the normalization condition δ S δ ( η ( r ) ) = δ S δ ( ξ ( r ) ) ξ ( r ) d r = N and for the conservation the energy of the system δ S δ ( β ( r ) ) ξ ( r ) d r = E Solution of this equation fully determine all thermodynamic parameter and describe general behavior of interacting system, whether this distribution of particles is spatially inhomogeneous or not. The above set of equations in principle solves the many-particle problem in thermodynamic limit. The spatially inhomogeneous solution of this equations correspondent the distribution of the interacting particles. Such inhomogeneous behavior is associated with the nature and intensity of the interaction. In other words, accumulation of particles in a finite spatial region (formation of a cluster) reflects the spatial distribution of the fields, the activity and temperature. Very important note, that only in this approach can take into account the inhomogeneous distribution temperature, which can depend from spatial distribution of particle in system.

3. Saddle-Point Equation for Self-Gravitating System

The solution of saddle-point equation completely determines all the thermodynamic parameters and describes the general behavior of a self-gravitating system both for spatially homogeneous and inhomogeneous particle distributions. The above set of equations in principle solves the many-particle problem in the thermodynamic limit. The spatially inhomogeneous solution of these equations corresponds to the distribution of the interacting particles. Such inhomogeneous behavior is associated with the nature and intensity of the interaction. In other words, accumulation of the particles in a finite spatial region (cluster formation) reflects the spatial distribution of the field, activity, and temperature. It is very important to note that only this approach makes it possible to take into account the inhomogeneity of temperature distribution that may depend on the spatial distribution of the particles in the system. In other approaches [2] the dependence of the temperature on a spatial point is introduced through the polytrophic dependence of temperature on particle density in the equation of state [18] . In the present approach, such dependence follows from the necessary thermodynamic condition, and can be found for various particle distributions. Now we derive the saddle-point equation for the extreme of the local entropy functional S ( ψ , ξ , β ) . The equation for the field variable δ S δ ψ = 0 in the case λ = 0 and absent repulsive interaction φ = 0 yields.

1 r m Δ ψ ( r ) + ξ ( r ) ( 2 π m 2 β ( r ) ) 3 2 β ( r ) exp ( β ( r ) ψ ( r ) ) = 0, (20)

where the notation r m 4 π G m 2 is introduced. The normalization condition may be written as

ξ ( r ) ( 2 m 2 β ( r ) ) 3 2 exp ( β ( r ) ψ ( r ) ) d r = N (21)

and the equation for the energy conservation in the system is given by

1 2 ( 2 π m 2 β ( r ) ) 3 2 ξ ( r ) β ( r ) ( 3 β ( r ) ψ ( r ) ) exp ( β ( r ) ψ ( r ) ) d r = E . (22)

To draw more information on the behavior of a self-gravitating system, we introduce new variables. The normalization condition ρ ( r ) d r = N yields the definition for the density function, i.e.,

ρ ( r ) ( 2 π m 2 β ( r ) ) 3 2 ξ ( r ) exp ( β ( r ) ψ ( r ) ) , (23)

which reduces the equations to a simpler form. The equation for the field variable is given by

Δ ψ ( r ) + r m β ( r ) ρ ( r ) = 0. (24)

In the case constant temperature and chemical activity, this equation transform into an equation for gravitational potential ψ = β ( r ) ψ in the well-known form:

Δ ψ ( r ) = 4 π G m 2 β ρ ( r ) . (25)

The equation for energy conservation takes the form

1 2 ρ ( r ) β ( r ) ( 3 β ( r ) ψ ( r ) ) d r = E . (26)

Hereinafter we obtain the chemical activity in terms of chemical potential ξ ( r ) = exp ( μ ( r ) β ( r ) ) . Having differentiated the equation for energy conservation over the volume, we obtain interesting relation for the chemical potential

1 2 ρ ( r ) β ( r ) ( 3 β ( r ) ψ ( r ) ) = δ E δ V δ V δ N = μ ( r ) ρ ( r ) , (27)

which yields the chemical potential to be given by

μ ( r ) β ( r ) = 3 2 1 2 β ( r ) ψ ( r ) . (28)

and the equation for the energy conservation reduces to very simple form: ρ ( r ) μ ( r ) = E . The equation thus obtained cannot be solved in the general case, but it is possible to analyze many cases of the behavior of a self-gravitating system under various external conditions. Within the context of the expression for density, and the definition of reduced thermal de-Broglie wavelength, and the gravitation length, i.e.,

Λ ( r ) = ( 2 β ( r ) 2 m e ) 1 2 , R g ( r ) = 2 π G m 2 β ( r ) , (29)

we can rewrite all the equations and the normalization condition in terms of density and temperature. Thus we have general equation in the term of concentration and temperature in the form:

Δ ( ln ( Λ 3 ( r ) ρ ( r ) ) β ( r ) ) + R g ( r ) β ( r ) ρ ( r ) = 0. (30)

The chemical potential reduces to

μ ( r ) β ( r ) = 3 2 ln ( Λ 3 ( r ) ρ ( r ) ) . (31)

In this way we can obtain the equation of state for self-gravitating system if we use the thermodynamic relation for conservation energy E,

P = 1 β δ S δ V . (32)

Using the definition of the density of particles, one can obtain the local entropy in the form

S = [ ρ ( r ) ln ( Λ 3 ( r ) ρ ( r ) ) ρ ( r ) ] d r , (33)

from which the local equation of state can present as

P ( r ) β ( r ) = ρ ( r ) ( 1 ln ( Λ 3 ( r ) ρ ( r ) ) ) = ρ ( r ) ( μ ( r ) β ( r ) 1 2 ) (34)

In general case this equation of state has a multiplier which logarithmically depends on the density of the particles. Only in the case Λ 3 ( r ) ρ ( r ) = 1 we obtain the equation of state for an ideal gas P β ρ . We have a sense to talk about pressure in classical cases only. If the concentration is large, and takes place the reverse relation 1 Λ 3 ( r ) ρ ( r ) , the determination of pressure are not correct. In the classical case Λ 3 ( r ) ρ ( r ) 1 . As shown below there is a natural limit of our approach. We cannot describe the processes which can be realize within short distances since because in a system with large concentration quantum effects would take place. In the case of the ideal gas we obtain usual equation of state, because in this case φ ( r ) = 0 and P β = ρ as result absent of interaction. In the case on ideal gas μ ( r ) β ( r ) = 3 2 and equation of state reproduce to the equilibrium relation for the ideal gas. In this case the energy of the system E = 3 2 N k T which corresponds to the previously obtained results [16] [17] . In the next section, we find the classical distributions of particles for various inner and external conditions.

4. Statistical Induced Dynamic for Self-Gravitating System

In system not far from equilibrium, patterns may form as results of the competition long-range interactions operating of different length scales [19] . The latter is countered by a long range gravitation attractive interaction. This is often accompanied by the macroscopical separation transition which leads to spontaneous formation of patterns in ideally homogeneous system upon variation of control parameter. An important class of system with competing interaction is the system in which the long-range interaction is gravitational. The fundamental nature of the gravitational interaction makes this class of system extremely diverse. We stop on the point statistical induced dynamics of self-gravitating system taking into account the possible spatial distribution of particle. For this goal we must determine the dynamic equation for field variable or density of particle. In this sense we can use the Ginsburg-Landau equation for density

ρ ( r , t ) t = 2 γ δ S δ ρ ( r ) , (35)

where γ is dynamic gravitational viscous coefficient [20] . In self-gravitating system we have not any other viscous process without dynamic influence of gravitational action, all systems on local spatial point. This evolution equation is, in fact applicable to a number of systems with non-conserved order parameter. Taking into account the presentation of local entropy in the terms of density, we can write the dynamic equation for density in the form:

ρ ( r , t ) t = 2 γ ln ( Λ 3 ( r ) ρ ( r , t ) ) (36)

The evolution of patterns in non-equilibrium will generate to be guided by the local entropy landscape and the morphological instabilities of the parameter. The dynamic of the system is dissipative, it will result in decrease of the local entropy of the patterns with time. To mimic this behavior we can use the simple gradient, descent dynamic defined with chemical potential may be considered as generalization of Cahn equation for non-uniform system with an arbitrary concentration gradient which becomes Cahn nonlinear diffusion equation [21] :

ρ ( r , t ) t = γ β μ ( r ) = 2 γ ln ( Λ 3 ( r ) ρ ( r , t ) ) (37)

which is fully equivalent previous dynamic equation in Ginsburg-Landau approach. If we take the solution of this dynamic equation in the form ρ ( r , t ) = exp ( γ r m t ) ρ ( r ) the solution of problem reduces to solve the previous equation

Δ ( ln Λ 3 ρ ( r ) ) + R g ρ ( r ) = 0 (38)

and can be transformed to

Δ ( ln ρ ( r ) ) + R g ρ ( r ) = 0. (39)

for constant temperature, that means, the obtained spatial solution will be existed in all time increasing the density of particles. If we determine the possible structure in distribution of density, this pattern will be conserved in all-time decrease of density. In this approach [7] , the probable behavior of a self-gravitating system can be predicted for any external conditions. On much longer timescales, an evolution towards the true thermal equilibrium is postulated. In this way we can solve the complicated problem of statistical description of self-gravitating systems. Moreover, this method can also be applied for the further development of the physics of self-gravitational and similar systems that are not far from the equilibrium.


The paper is prepared with the financial support of the National Academy of Sciences of Ukraine (the project “Dynamics of structure formation in many-particle systems”) PK 0118U003535.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

Cite this paper

Lev, B.I. (2019) Statistical Induced Dynamic of Self-Gra- vitating System. Journal of Modern Physics, 10, 687-698.


  1. 1. Thirring, W. (1970) Zeitschrift für Physik, 235, 339.

  2. 2. Chavanis, P.-H., Rosier, C. and Sire, C. (2002) Physical Review E, 66, Article ID: 036105.

  3. 3. Laliena, V. (2003) Nuclear Physics B, 668, 403.

  4. 4. Bilotsky, Y.D. and Lev, B.I. (1984) Teoreticheskaya i Matematicheskaya Fizika, 60, 120.

  5. 5. Lev, B.I. and Ya, A. (1998) Physical Review E, 57, 6460-6469.

  6. 6. Lev, B.I. (2011) International Journal of Modern Physics B, 25, 2237-2249.

  7. 7. Lev, S.B. and Lev, B.I. (2017) The European Physical Journal B, 90, 3.

  8. 8. Kleinert, H. (1989) Gauge Field in Condensed Matter. Word Scientific, Singapore.

  9. 9. Stratonovich, R.L. (1958) Sov. Phys. Dokl., 2, 416-419.

  10. 10. Hubbard, J. (1958) Phys. Rev. Lett., 3, 77.

  11. 11. Grigorishin, K.V. and Lev, B.I. (2005) Physical Review, 71, Article ID: 066105.

  12. 12. Pakter, R., Marcos, B. and Levin, Y. (2013) Physical Review Letters, 111, Article ID: 230603.

  13. 13. Benetti, F.P.C., Ribeiro-Teixeira, A.C., Pakter, R. and Levin, Y. (2014) Physical Review Letters, 113, Article ID: 100602.

  14. 14. Baxter, R. (1980) Exactly Solved Models in Statistical Mechanics. Academic Press, New York.

  15. 15. Ruelle, D. (1969) Statistical Mechanics. Rigorous Results, New York, Amsterdam.

  16. 16. Zubarev, D.N. (1974) Non-Equilibrium Statistical Thermodynamics. Consultants Bureau, New York.

  17. 17. Huang, K. (1963) Statistical Mechanics. J. Wiley and Sons, New York.

  18. 18. Isihara, A. (1971) Statistical Mechanics. State University of New York, New York.

  19. 19. Sire, C. and Chavanis, P.-H. (2002) Physical Review E, 66, Article ID: 046133.

  20. 20. Muratov, C.B. (2002) Physical Review E, 66, Article ID: 066108.

  21. 21. Chandrasekhar, S. (1942) An Introduction to the Study of Stellar Structure. Dover Publications, New York.

  22. 22. Onuki, A. (1989) Journal of the Physical Society of Japan, 58, 3065-3068.