Copyright 1980-2015 John Bruce Davies

This research was awarded an Honorable Mention in the 1980 Gravity Research Foundation essay competition.


We shall approach the production of an equation of state from a purely continuum model with no assumption about its constituents and no assumption of a specific model of matter. We instead impose the condition that any real equation of state must govern a body whose internal energy is at a minimum and thus at equilibrium. By considering an isothermal self-gravitating body in equilibrium and thus at a minimum of internal energy, a new field equation is obtained governing pressure, density, and the self-gravitational field. An asymptotic method is used for solution of this new equation when combined with the hydrostatic equilibrium equation and the Poisson equation of self-gravitation. The solution yields the unique determination of the classical equation of state, the only necessary information input being the density, cohesive energy and bulk modulus of the material at one particular pressure.


Ever since Newton’s associate Boyle experimented on gases, the effect of changing pressure and temperature on volume has been of paramount importance. Boyle showed that his equation of state for small changes gave volume inversely proportional to pressure and directly proportional to temperature. As the atomic description of matter progressed, physicists developed models of equations of state based on the atomic and molecular structures. At constant temperature, the pressure is just the negative derivative of the substance’s energy with respect to volume. Thus, by invoking an atomic description with electrons at various radii around nuclei, energy as a function of volume was specified in such models. This gave a quantum mechanical description of the equation of state of the particular model of the matter under consideration.

It will be shown that the continuum self-gravitational field of matter can be used to obtain the classical equation of state which relates the macrophysical parameters of pressure and density at constant temperature. This approach is vitally different as the equation of state of a material is usually imposed by the physicist's atomic or molecular model of the body. For isotropic, isothermal bodies of a single constituent it will be shown that the equation of state, in the Newtonian non-relativistic domain, is uniquely determined from gravitational considerations and the minimum energy properties of the equilibrium state. The classical equation of state is applicable to equilibrium; a necessary though not sufficient condition for such being satisfaction of the equation of hydrostatic equilibrium. This hydrostatic equation relates pressure gradients to the gravitational field in a self-gravitating body. We will mathematically construct a self-gravitating sphere which obeys the Poisson equation that relates the gravitational field to the density of matter. The three independent variables, pressure, density and gravitational field, obey the above two equations. They are usually obtained as functions of radius of the sphere by imposition of an equation of state hypothesized for that particular material.

However, the fundamental property of an equilibrium state is that the total internal energy is a minimum at constant entropy. Solution of this variational problem allows a third equation to be derived relating pressure, density, and the gravitational field. Solution of these three equations yields the unique determination of the classical equation of state, pressure versus density, given necessary and sufficient boundary conditions. We thus show that it is not necessary to hypothesize the classical equation of state--it is uniquely determined from self-gravitational and equilibrium considerations alone.

Governing Equations

In the non-relativistic domain, we shall investigate an isothermal, isotropic, spherically symmetric material body of a single constituent under hydrostatic pressure and self-gravitating forces. This obeys the necessary, though not sufficient, condition of hydrostatic equilibrium:

(1) dP/dr = ρ dψ/dr

where P is hydrostatic pressure, ρ is density, ψ is gravitational potential and r is the independent variable, the radius. The equation of self-gravitation is:

(2) d2ψ/dr2 + 2/r dψ/dr = -4πGρ

where G is gravitational constant. The body will have a particular mass, M, which is taken as constant, so that:

M = R 4πρ r2.dr

where R is the unknown radius of the body. As this body is in equilibrium, the total internal energy is a minimum. The total energy consists of the gravitational potential energy, the strain energy and other constant energies.

The potential energy is given by:

(3) EG = R 4πρ r3r.dr

and the strain energy is:

(4) ES = 0R 4πρr2.ε(r).dr

where ε(r) is strain energy per unit mass given by the work done in compression:

(5) ε(r) = ρ(0)ρ(r) P(ρ).dρ/ρ2

where ρ(0) is the surface density at radius R. The remaining energies are the cohesive, thermal and other energies which are constant and thus independent of radius. Differentiating (5) gives:

(6) P.dρ/dr = ρ2.dε/dr

New Field Equation

The total internal energy, comprising the gravitational potential energy (3) plus the strain energy (4) is now minimized, subject to the subsidiary conditions (1), (2) and (6) and with constant mass. Using Lagrange Undetermined Multipliers λ1(r), λ2(r), λ3(r) and μ, the variational statement is:

(7) δ.0R4πρr3ψr+4πρr2ε+λ1(Pr-ρψr)+λ2(Pρr2εr)+λ3rr+2/r.ψr+4πGρ)+μ4πρr2.dr = 0

This states that equilibrium is achieved due to the balance of gravitational potential energy and strain energy subject to satisfaction of the governing differential equations. The necessary condition for a minimum is the Euler-Lagrange equations i.e. for

δ.0RF.dr = 0

where F contains everything under the integral in (7), where the independent variables are ρ, ψr, ε and P, so that the equations are:

δF/δP - d/dr.(δF/δPr ) = 0

δF/δρ - d/dr.(δF/δρr ) = 0

δF/δψr - d/dr.(δF/δψrr ) = 0

δF/δε - d/dr.(δF/δεr ) = 0

Performing these operations and removing the Lagrangian Undetermined Multipliers by algebraic methods results in the fundamentally new equation:

(8) εr = -K1G/2r2 - ψrK2G/2πr3

where K1 and K2 are integration constants. Using equations (6) and (8), we obtain the new original equation relating P, ρ and ψr in isothermal equilibrium systems:

(9) ψr = Ar/(B+r3) + 3r3/4(B+r3).Pρr2

where A and B are constants with B=3/8.GK1.

We thus have three equations (1),(2) and (9) for three dependent variables, which, although non-linear, will yield unique values of P(r), ρ(r) and ψ(r), given necessary and sufficient boundary conditions.

As the gravitational potential enters only in its gradients, we can put φ = ψr and rearrange the three equations in the standard form for a set of ordinary differential equations. From (1):

(10) dP/dr = ρφ

From (2):

(11) dφ/dr = -2/r.φ - 4πGρ

From (9), our original equation now becomes:

(12) dρ/dr = ρ2/P.4/3r3.{(B+r3)φ – Ar}

This is the fundamentally new equation that governs a self-gravitating body in equilibrium and specifies its density gradient in terms of density, pressure and self-gravitational force.


These three equations comprise a non-linear system; no exact solutions are presently known and thus approximate solutions will be obtained. First normalize the dependent variables with respect to their surface values, P(0), φ(0), ρ(0), and the independent variable r is also normalised so that s = r/R where s ranges from 0 to 1. Thus the new normalised variables are taken as:

Q(s) = P/P(0), g(s) = φ/φ(0), w(s) = ρ/ρ(0).

The hydrostatic equation (10) then becomes:

(13) dQ/ds = w.g.L1


(14) L1 = ρ(0).φ(0).R/P(0)

The self-gravitation equation (11) becomes:

(15) dg/ds = -2/s.g – w.L2


(16) L2 = 4πGρ(0)/φ(0)

and our new equation (12) becomes:

(17) dw/ds = w2/Q .4/3s3.L1{(σ + s3)g – σL3s}


(18a) σ = B/R3

(18b) σ L3 = Aρ(0)/{R.P(0)}

Asymptotic Solutions

In the case , B >>R3 , then equations (15) and (17) show that this implies σ -> infinity. Thus an asymptotic method immediately suggests itself by expansion of the dependent variables in inverse orders of σ.

(19) w = w0 + 1/σ.w1 + 1/σ2.w2 + ...

Q = Q0 + 1/σ.Q1 + 1/σ2.Q2 + ...

g = g0 + 1/σ.g1 + 1/σ2.g2 + ...

Substituting these expansions into the three equations (13), (15) and (17), equating terms with the same order of σ , and solving the resulting equations, gives to order 1/σ2 :

(20) w0 = -3L3/L1L2

w1 = 6L3/L1L2.s3

w2 = [-21/2.N0/L3 s4 +45/4.L3/L1L2.s6]

(21) g0 = L3/L1.s

g1 = -L3/L1.s4

g2 = [-5/4.L3/L1 s5 +3/2.N0L2/L3.s7]

(22) Q0 = N0 -3/2.(L3)2/L2L1.s2

Q1 = N1 -9/5.(L3)2/L2L1.s5

Q2 = N2 -5/2.N0s6 +9/8.(L3)2/L2L1.s8

where N0, N1 and N2 are integration constants.

Boundary Conditions

We will now impose surface boundary conditions in order to obtain the relationship between density ρ(c) and pressure P(c) at the center of the sphere, in terms of the parameter set. The constants in equations (20), (21), and (22) will thus be determined by the surface values.

At the center, s = 0, we have from (20) and (22):

(23) ρ(c)/ρ(0) = -3L3/L1L2

(24) P(c)/P(0) = N0 + 1/σ.N1 + 1/σ2.N2

Similarly at the surface, s = 1, another set of equations is determined, which, after some manipulation, relates L1, L2, L3, to σ, P(c) and P(0). The constants, L1 and L2, are dependent on φ(0), the force of gravity at the surface. This can be obtained in terms of the known bulk modulus K(0) of the material at the surface by defining:

(25) D = K(0)/P(0) = ρ(0)/P(0).[dP/dρ]0

and using (17) and (13), we obtain:

(26) L3/L1 = 1 + [1-3/4D]/σ

From equation (23), the ratio of central to surface density can be obtained using the surface boundary conditions where (20) and (21) are calculated at s = 1:

(27) ρ(c)/ρ(0) = L3/3L1.[7 - L3/L1.{4-1/σ+2.5/σ2}]

Similarly, the ratio of surface to central pressure is obtained using (24) and the surface boundary conditions:

(28) P(0)/P(c) = 1 - 1.5/σ2.[L3/4.{3-9/5σ-9/8σ2}

+5/3.[1-L3/L1.{1-1/5σ+5/4σ2 }] / [1-L3/L1.{1-1/σ+5/4σ2 }]

The density ρ(0) is the density of the matter at the surface conditions.

The surface pressure P(0) is not however the applied pressure at the surface alone; it also incorporates the effective pressure due to the cohesive energy of the material. Thus when the particles of matter are brought from infinity to form the solid (or fluid), a cohesive energy of the material is implicit in the total strain energy.

The effective internal pressure P(i) is thus the negative of the tension necessary to break the bonds of attraction, so that if the externally applied pressure is P(ext), then:

P(0) = P(i) + P(ext)

In practice, this P(i) is taken as the cohesive energy per unit volume:

P(i) = CE.ρ(0)/Mol

where CE is the cohesive energy/mole and Mol is the gram molecular weight of the material.


We compare the predicted isothermal equation of state for the elements, Copper and Aluminum, with that derived from shock-wave measurements in Figure 1, Clarke (1966). The two approaches are in good agreement, though it must be remembered that shock wave derived information loses accuracy on reduction to the isotherm.

Also, different material samples of the same element do not have identical values of density and bulk modulus at the zero external pressure thereby supplying additional error.


We have shown that mathematical construction of an isothermal self-gravitating body in equilibrium yields a new field equation relating pressure, density and the

self-gravitational field. This, together with the hydrostatic equation and Poisson's equation yield three equations which are solved asymptotically to give the unique determinaton of the classical equation of state. The only information necessary is the

density and bulk modulus at the surface pressure, and the cohesive energy of the particular material.


Clark, S.P. ed., 1966. Geological Soceity of America Memoir 97.