# The fundamental equilibrium equation for gaseous stars and the Tolman-
Oppenheimer-Volkoff equation–derivations and applications with emphasis on
optimizational-variational methods.

Steven D. Miller  Rytalix Analytics, Strathclyde, Scotland  stevendM@ed-
alumnus.net

###### Abstract.

Stars are essentially gravitationally stabilised thermonuclear reactors in
hydrostatic equilibrium. The fundamental differential equation for all
Newtonian gaseous stars in equilibrium is

$\displaystyle\frac{dp(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}$

where $p(r),\rho(r)$ are the pressure, density at radius $r$ and
$\mathcal{M}(r)$ is the mass contained within a shell of radius $r$ given by
$\mathcal{M}(r)=\int_{0}^{r}4\pi\overline{r}^{2}\rho(\overline{r})d\overline{r}$,
and $\mathscr{G}$ is Newton’s constant. This simple but crucial differential
equation for the pressure gradient within any star, underpins much of
astrophysical theory and it can derived by various methods:via a simple
heuristic argument; via the Euler-Poisson equations for a self-gravitating
fluid/gas; via a variational method by taking the 1st variation of the sum of
the thermal and gravitational energies of the star; via the 2nd variation of
the Massiue thermodynamic functional for a self-gravitating isothermal
perfect-gas sphere; from conservation of the virial tensor; as the non-
relativistic limit of the Tolman-Oppenheimer-Volkoff equation (TOVE). The TOVE
for equilibrium of relativistic stars in general relativity can in turn be
derived by various methods: from the energy-momentum conservation constraint
on the Einstein equations applied to a spherically symmetric perfect
fluid/gas; via a constrained optimization method on the mass and nucleon
number; via a maximum entropy variational method for a sphere of self-
gravitating perfect fluid/gas or radiation. An overview is given of all
derivations with emphasis on variational methods. Many important applications
and astrophysical consequences of the equilibrium equation are also reviewed.

###### Contents

  1. 1 Introduction
  2. 2 Heuristic Derivation Of The FEHE–Its Applications And Consequences
    1. 2.1 Heuristic Derivation of the FEHE
      1. 2.1.1 Fundamental stellar structure equations and nuclear reactions
      2. 2.1.2 Central pressure and temperature estimates
    2. 2.2 Hydrostatic equilibrium in higher dimensions
      1. 2.2.1 Some basic theorems concerning hydrostatic equilibrium
    3. 2.3 Gravitational and Thermal Energies of a Star in Hydrostatic Equilibrium
    4. 2.4 The stellar virial theorem
    5. 2.5 The isothermal sphere for a perfect gas
      1. 2.5.1 Relevance of the isothermal gas sphere to red giant evolution
  3. 3 Deviations From Hydrostatic Equilibrium and Stellar Timescales
    1. 3.1 Zero-pressure gravitational collapse and the hydrodynamical time
    2. 3.2 The Kelvin-Helmholtz time scale
  4. 4 Derivation Of The Hydrostatic Equilibrium Equation Via Variational Methods
    1. 4.1 2nd variation of the thermodynamic Massieu function for an isothermal self-gravitating sphere of perfect gas
  5. 5 Polytropic Gaseous Stars And The Lane-Emden Equation:Derivation, Existence And Uniqueness Of Solutions, Analytical Solutions And Stability
    1. 5.1 The Lane-Emden equations
      1. 5.1.1 Non-dimensional reduction
    2. 5.2 Existence and uniqueness
    3. 5.3 Transformations of the LE equations and a dynamical systems perspective
    4. 5.4 explicit solutions
    5. 5.5 The isothermal gas sphere again
  6. 6 The Fundamental Hydrostatic Equilibrium Equation (FHEE) And The Euler-Poisson Equations For A Self-Gravitating Gaseous Star
    1. 6.1 Hydrostatic equilibrium via the virial tensor
    2. 6.2 Energy integrals for $EP_{\gamma}$ systems in static equilibrium in all dimensions
    3. 6.3 Bounds on the pressure integral
  7. 7 Perturbations And (Linear) Stability Criteria For Static Gaseous Stars In Hydrostatic Equilibrium
    1. 7.1 Equilibrium as a critical point of the total energy
    2. 7.2
  8. 8 Relativistic Stars And The Tolman-Oppenheimer-Snyder Equation
    1. 8.1 Einstein equations and their solution for a spherically symmetric perfect fluid/gas
    2. 8.2 Some properties of the TOVE
  9. 9 Derivation Of The TOVE By Constrained Optimisational-Variational Methods
  10. 10 Maximum Entropy Derivations Of The Tolman-Oppenheimer-Volkoff Equations For Self-Gravitating Radiation And Relativistic Fluids/Gases
    1. 10.1 Maximum entropy derivation the TOVE for self-gravitating radiation
    2. 10.2 Thermodynamic derivation of the TOVE for a perfect fluid/gas as a constrained optimisational-variational problem
    3. 10.3 Integrability of the TOVE for incompressible stars and the sharp Buchdal bound

##  1\. Introduction

Stars are essentially gravitationally stabilized thermonuclear reactors in
hydrostatic equilibrium. If stars where not in near-perfect hydrostatic
equilibrium for billions of years, stable conditions permitting life to evolve
on Earth (and probably elsewhere in the Universe) would simply not have
existed. Hydrostatic equilibrium is an absolutely necessary condition for
stable luminosity, temperature and thermonuclear fusion over vast periods of
cosmic time. It is a fundamental and crucial property of stars, providing an
exact balance between the gravitational force which attracts the matter toward
the center, and the force due to the gas, the thermal and radiative pressure
pushing outwards, which perfectly counterbalances gravity. Any slight
deviation from this equilibrium will immediately lead to a rapid reaction to
induce restoration of the equilibrium.

Let us suppose, for example, that a gaseous star is arbitrarily compressed to
a smaller radius. The gas becomes hotter rising the internal pressure. The
higher pressure provokes an expansion re-settling the star back to its
equilibrium state. Conversely, an arbitrary extension of the radius would
decrease the internal temperature and pressure, the star would then contract
again due to its self gravitation. Such re-adjustments are very fast,
occurring at the dynamical timescale of the order of half an hour for a
typical main sequence star like the Sun, which is effectively instantaneous
with respect to typical stellar lifetimes. The hydrostatic equilibrium of a
any star therefore governs all its properties. It is always satisfied except
in very short phases, such as the initial collapse of interstellar clouds or
when it has exhausted its nuclear fuel and either collapses or supernovas.

In this article, the fundamental equation of hydrostatic equilibrium (FEHE) is
derived and its consequences are studied in detail. This simple but crucial
equation underpins most of astrophysical theory and in the literature is
usually only derived heuristically by very simple arguments. However, here the
various methods by which the equation can be derived are given in much greater
and more rigorous mathematical detail, with an emphasis on variational methods
whereby HE corresponds to an extremum or critical point of the total energy of
the star. The relativistic extension of the FEHE–the Tolman-Oppenheimer-
Volkoff equation–is also derived by the standard method via the Einstein
equations coupled to a perfect fluid/gas, but also from various constrained
optimisational-variational methods, including maximum entropy methods. The
motivation is due to a dissatisfaction from a lack of mathematical rigour and
detail in most textbooks, papers and online lecture notes on astrophysics. (In
the author’s opinion.) Often derivations and calculations are only roughly
outlined and one is obliged to do the calculation or proof oneself filling in
all the details.

Although this is predominantly an autodidactic-type review article, it does
also contain a few original perspectives and derivations and fills in many
mathematical details omitted within the literature. As such, a great deal of
fundamental astrophysical theory is covered including the theory of polytropic
gaseous stars and the fundamentals of Lagrangian perturbation theory. In
particular, the Euler-Poisson equations for a self gravitating fluid/gas are
considered in detail purely for the equilibrium case. Key applications and
astrophyical consequences of the hydrostatic equilibrium equation are
reviewed: the virial theorem, incompressible stars, stability and
perturbations, lower-bound temperature and pressure estimates in stellar
cores, polytropic gas spheres and the Lane-Emden equation, isothermal gas
spheres and red giants, and the sharp Buchdal bound for incompressible
relativistic stars.

##  2\. Heuristic Derivation Of The FEHE–Its Applications And Consequences

In this section a (standard) heuristic derivation of the equations of
hydrostatic equilibrium for self-gravitating gaseous stars is presented and
its many important and crucial applications and consequences in astrophysical
theory are briefly reviewed. First the relevant variables are defined.

######  Definition 2.1.

The scenario describing the self-gravitating gas/fluid is defined as follows.

  1. (1)

Generally, we use the functions
$\rho:[0,T)\times\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq
0},p:[0,T)\times\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq 0}$ and
$\Phi:[0,T)\times\mathbf{R}^{3}\rightarrow\mathbf{R}$ for some temporal
interval $[0,T)$, to describe the spatio-temporal distributions of mass
density, pressure and Newtonian potential of the gas/fluid. Then
$\rho=\rho(\mathbf{x},t),p=p(\mathbf{x},t),\Phi=\Phi(\mathbf{x},t)$.

  2. (2)

The self-gravitating gas is supported within a spherical region
$\mathbf{D}(\lambda(t))\equiv\mathbf{D}(t)$ defined by

$\displaystyle\mathbf{D}(\lambda(t))=\mathbf{D}(t)=supp~{}\rho(t,\bullet)=\\{(\mathbf{x}\in\mathbf{R}^{3}|\rho(\mathbf{x},t)>0,\lambda(t)\gtrless\lambda(0)=R\\}\subset\mathbf{R}^{3}$

with varying radius $\lambda(t)\gtrless\lambda(0)=R$ depending on whether the
self-gravitating gaseous sphere is collapsing or expanding, or
$\lambda(t)=\lambda(0)=R$ in hydrostatic equilibrium. Then
$\mathrm{V}(\mathbf{D}(t))=\tfrac{4}{3}\pi|\lambda(t)|^{3}$ is the volume at
time t.

  3. (3)

The (moving) boundary is $\partial\mathbf{D}(t)\subset\mathbf{D}$ with surface
area $\mathrm{A}(\mathbf{D}(t))=4\pi|\lambda(t)|^{2}$. The boundary conditions
are that the pressure and density vanish here so that
$\rho(\mathbf{x},t)=p(\mathbf{x},t)=0$ for all $x\in\partial\mathbf{D}(t)$. If
$\mathbf{D}(t)\bigcup\mathbf{D}^{c}(t)=\mathbf{R}^{3}$ then the vacuum region
is defined by $\rho(\mathbf{x},t)=0$ in $\mathbf{D}^{c}(t)$ and the matter or
hydrodynamic region by $\rho(\mathbf{x},t)>0,x\in\mathbf{D}(t)$. The fluid
velocity $\mathbf{u}$ is defined only on $\mathbf{D}(t)$ by $\mathbf{u}\equiv
u^{i}:[0,T)\times\mathbf{D}(t)\rightarrow\mathbf{R}^{3}$ and there is no
definition in $\mathbf{D}(t)$. The space $\mathbf{D}(0)$ is taken to be
precompact. A diffuse boundary condition can also be defined as (ref)

$\displaystyle\lim_{(\mathbf{x},t)\rightarrow(\overline{\mathbf{x}},\overline{t})}\rho(\mathbf{x},t)\nabla_{i}p(\mathbf{x},t))=0$
(2.1)

where $(x,t)\in[0,T)\times\mathbf{D}(t)$ and
$(\overline{x},\overline{t})\in[0,T)\times\partial\mathbf{D}(t)$, where again
$\partial\mathbf{D}(t)$ is the boundary of the volume $\mathbf{D}(t)$, that
is, the moving interface between the fluid/gas matter and the exterior vacuum.
This boundary condition characterizes the gradually vanishing fluid pressure
per mass near the boundary.

  4. (4)

The initial data are

$\displaystyle\rho(0,\mathbf{x})=\rho_{o}(\mathbf{x}),~{}~{}p(0,\mathbf{x})=p_{o}(\mathbf{x})~{}~{}\mathbf{x}\in\mathbf{R}^{3}$
(2.2)
$\displaystyle\mathbf{u}(0,\mathbf{x})=\mathbf{u}_{o}(\mathbf{x}),~{}~{}\mathbf{x}\in\mathbf{D}(0)$
(2.3)

  5. (5)

The Newtonian potential $\Phi$ exists on
$\mathbf{R}^{3}=\mathbf{D}(t)\bigcup\mathbf{D}^{c}(t)$, in both the matter and
vacuum regions and obeys the Poisson equation

$\displaystyle\Delta\Phi=4\pi\mathscr{G}\rho,~{}~{}~{}in~{}~{}\mathbf{R}^{3}$
(2.4)

with the general solution

$\displaystyle\Phi(\mathbf{x},t)=-\mathscr{G}{\int}_{\mathbf{D}(t)}\frac{\rho(\mathbf{y},t)d^{3}\mathbf{y}}{|\mathbf{x}-\mathbf{y}|},~{}~{}~{}~{}for~{}(t,\mathbf{x},\mathbf{y})\in[0,T)\times\mathbf{R}^{3}$
(2.5)

  6. (6)

It is generally assumed that the gas is _isentropic_ having constant entropy
or entropy per nucleon, and obeys a barytropic pressure law of the form
$p=p(\rho)$, typically a polytropic gas law

$\displaystyle
p=\mathsf{K}|\rho|^{\gamma},~{}~{}\mathsf{x}\in\mathbf{D}(t)\in\mathbf{R}^{3}$
(2.6)

(Polytropic stars will be considered in detail in Section (5) The set
$\\{\rho(x,t),p(x,t),\Phi(x,t)\\}$ is also a solution to the Euler-Poisson
system of equations but this will be considered formally in Section (6). For a
gas in hydrostatic equilibrium $\mathbf{u}$=0 and
$\rho:\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq
0},p:\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq 0}$ and
$\Phi:\mathbf{R}^{3}\rightarrow\mathbf{R}$ so that
$\rho=\rho(\mathbf{x}),p=p(\mathbf{x}),\Phi=\Phi(\mathbf{x})$. The volume of
the support is now fixed so that
$\mathrm{V}(\mathbf{D}(t))\mathrm{V}(\mathbf{D}(0))=\tfrac{4}{3}\pi
R^{3}\equiv\tfrac{4}{3}\pi|\lambda(0)|^{3}$.

  7. (7)

For self-gravitating fluid-gas systems like stars in equilibrium at radius
$R$, it is also natural to assume spherical symmetry and coordinates so that
$\rho=\rho(r),p=p(r),\Phi=\Phi(r)$ with central values $\rho(0),p(0)$ and
boundary values $\rho(R)=p(R)=0$. The domain
$\mathbf{D}(\lambda(t))=\mathbf{D}(\lambda(0))\equiv\mathbf{D}(0)=\mathbf{B}(\mathbf{x}_{o},R)$
is now a static ball with centre $\mathbf{x}_{o}$ an radius $R$. Then

$\displaystyle\mathbf{\mathbf{D}}(\lambda(0)=R)=\mathbf{B}(0,R)=\\{\mathbf{x},\mathbf{y}\in\mathbf{R}^{3}|d(x,y)\leq
R\\}$ (2.7)

The mass contained within a radius $r$ is or ball
$\mathbf{B}(0,r)\subset\mathbf{B}(0,R)$ is then $\mathcal{M}(r)$ and
$\mathcal{M}(R)=M$, the total mass contained within $\mathbf{B}(0,R)$ It is
defined as

$\displaystyle\mathcal{M}(r)=\int_{\mathbf{B}(0,r)}4\pi\widehat{r}^{2}\rho(\widehat{r})d{r}\equiv\int_{0}^{r}4\pi\widehat{r}^{2}\rho(\widehat{r})d\overline{r}$
(2.8)

Within this scenario, the fundamental equation(s) of hydrostatic equilibrium
(FEHE) can be derived heuristically.$\bm{[1,2,3]}$ The nonlinear Euler-Poisson
PDEs describing this self-gravitating system can also be formally defined but
this will be done in Section (6).

###  2.1. Heuristic Derivation of the FEHE

The total mass of the gas is $M$ and the mass of fluid/gas contained within a
spherical surface of radius $r$ is $\mathcal{M}(r)$. The mass of gas contained
within a thin shell of thickness $dr=(r+dr)-r$ is $d\mathcal{M}(r)=4\pi
r^{2}\rho(r)$. The mass within the shell experiences a gravitational force
$\mathscr{F}_{grav}$ inward due to the mass of gas $\mathcal{M}(r)$. Applying
Newton’s inverse square law to the shell

$\displaystyle\mathscr{F}_{grav}=-\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r^{2}}=-\frac{\mathscr{G}(4\pi
r^{2}\rho(r)dr)\mathcal{M}(r)}{r^{2}}$ (2.9)

Due to the intrinsic pressure of the gas, the shell also experiences a
buoyancy force pushing radially outwards. The outward pressure on the shell is
the difference in pressures at the interior and exterior boundaries of the
shell at r and r+dr so that

$\displaystyle\mathscr{F}_{buoy}=4\pi r^{2}(p(r)=p(r+dr))-4\pi
r^{2}\frac{dp(r)}{dr}dr$ (2.10)

and $dp(r)/dr=(p(r+dr)-p(r))dr$ for small $dr$. The equilibrium of the shell
requires that these two forces balance so that
$\mathscr{F}_{grav}+\mathscr{F}_{buoy}=0$.

$\displaystyle=-\frac{\mathscr{G}(4\pi
r^{2}\rho(r)dr)\mathcal{M}(r)}{r^{2}}+4\pi r^{2}\frac{dp(r)}{dr}dr=0$ (2.11)

This then gives the fundamental equation of hydrostatic equilibrium for all
Newtonian gaseous stars

$\displaystyle\underbrace{\frac{dp(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}}$
(2.12)

######  Remark 2.2.

Throughout the paper, whenever the FEHE (2.12) is derived for the first time
from a new method or technique, it will be emphasised with an underbrace.

Equations (2.8) and (2.12) are two of the fundamental equations of stellar
structure theory.(REFS) The total pressure $p(r)$ in the star is the sum of
the gas pressure $p_{gas}(r)$ and the photon or radiation pressure
$p_{rad}(r)$ so that $\bm{[1]}$

$\displaystyle\frac{d}{dr}(p_{gas}(r)+p_{rad}(r))=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}$
(2.13)

where for a mixture of radiation and a perfect gas of ionised hydrogen

$\displaystyle p(r)=p_{gas}(r)+p_{rad}(r)=\frac{k_{B}}{\mu
m_{p}}\rho(r)\Theta(r)+\frac{1}{3}a\Theta^{4}$ (2.14)

where $\Theta(r)$ is the temperature of the gas, $k_{B}$ is the Boltzmann
constant,’a’ is the Stefan constant, $\mu$ is the mean molecular weight, and
$m_{p}$ is the mass of the proton. One can also write
$p_{gas}=\mathscr{R}\rho\Theta$ where $\mathscr{R}=\frac{k_{B}}{\mu m_{p}}$ is
the fundamental gas constant. The gradient of the radiation pressure is
determined by the radiant photon flux $\mathcal{L}(r)$ through a spherical
surface by

$\displaystyle\frac{dp_{rad}(r)}{dr}=-\frac{\kappa\mathcal{L}(r)}{4\pi
cr^{2}}\rho(r)$ (2.15)

where

$\displaystyle d\mathcal{L}(r)=4\pi r^{2}\rho\mathfrak{T}dr$ (2.16)

and $\kappa$ is the Krammer’s opacity $\bm{[2,3]}$. Since $\mathfrak{T}$ is
essentially the thermonuclear energy production rate per mass, we can see that
the effect of thermonuclear reactions is to maintain the photon/radiation
pressure gradient and thus hydrostatic equilibrium.

The hydrostatic equilibrium equations can be easily solved for an
incompressible gas star when $\rho(r)$ is constant and for density profiles of
the form$\rho(r)=f(r/R)$.

######  Lemma 2.3.

Let a gaseous star fill a domain $\mathbf{B}\subset\mathbf{R}^{n}$ with
$V(\mathbf{B})=\tfrac{4}{3}pR^{3}$. If $\rho(r)=\rho(0)=\rho_{c}=const.$ then
(2.8) and (2.12) can be solved to give
$\mathcal{M}r)=\tfrac{4\pi}{3}\rho_{c}r^{3}$ and the pressure profile has the
form

$\displaystyle
p(r)=\frac{1}{2}\rho_{c}\beta\left(1-\left|\frac{r}{R}\right|^{2}\right)$
(2.17)

with $\beta=3M/8\pi R^{3}$. The central pressure is
$p_{c}=p(0)=\tfrac{1}{2}\rho_{c}\beta$.

###### Proof.

Equation (2.8) can be easily integrated to give
$\mathcal{M}(r)=\tfrac{4}{3}\pi r^{3}\rho_{c}$ and
$M=\mathcal{M}(R)=\tfrac{4}{3}\pi R^{3}\rho_{c}$. Then the FEHE (2.12) becomes

$\displaystyle\frac{dp(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)|\rho_{c}|^{2}}{r^{2}}=-\frac{4}{3}{\mathscr{G}}\pi
r|\rho_{c}|^{2}$ (2.18)

Integrating

$\displaystyle\int_{p(r)}^{p(R)}dp(r)$
$\displaystyle=-p(r)=\left|-\frac{2}{3}\pi{\mathscr{G}}\rho_{c}r^{2}\right|_{r}^{R}=-2\pi{\mathscr{G}}|\rho_{c}|^{2}R^{2}\left(1-\left|\frac{r}{R}\right|^{2}\right)$
$\displaystyle=\frac{1}{2}|\rho_{c}|^{2}\left(1-\left|\frac{r}{R}\right|^{2}\right)$
(2.19)

∎

######  Theorem 2.4.

The gravitational potential $\Phi(r)$ for a ball of fluid/gas supported within
$\mathbf{D}=\mathbf{B}_{3}(R)\subset\mathbf{R}^{3}$ of radius R of uniform
mass density $\rho$ inside the matter, and outside the ball in vacuum is

$\displaystyle\Phi(x)=\frac{C}{4\pi}\int_{\mathbf{D}}\frac{\rho
d^{3}y}{|x-y|}=\rho\left(\frac{1}{2}R^{2}-\frac{1}{6}\|x|\|^{2}\right),~{}~{}\|x\|\leq
R$ (2.20) $\displaystyle\Phi(x)=\frac{C}{4\pi}\int_{\mathbf{D}}\frac{\rho
d^{3}y}{|x-y|}=\frac{C\rho R^{3}}{3\|x\|},~{}~{}\|x\|>R$ (2.21)

Outside the ball, the potential behaves as thought it is due to a point mass
with all the mass of the ball concentrated at the point of origin. This
property also enables gravitational interactions of stars and planets within
celestial mechanics to be reduced to the problem of interacting point masses.
The result holds for all dimensions $n\geq 4$ and also for logarithmic
potentials of discs within $\mathbf{R}^{2}$. When the density is constant, the
potential across the surface $\partial\mathbf{B}(0)$ is constant.

####  2.1.1. Fundamental stellar structure equations and nuclear reactions

The hydrostatic equilibrium equations (2.8) and (2.12) are two of the four
fundamental first-order differential equations of stellar structure theory,
the other two describing radiative transport. $\mathcal{L}(r)$ is the total
luminosity or radiant energy flux through a spherical surface, $\kappa(r)$ the
Kramers opacity and $\mathfrak{T}(r)$ the thermonuclear energy production rate
per mass. With the chemical composition fixed, then these quantities are taken
to be fixed functions of the density $\rho(r)$ and temperature $\Theta(r)$ of
the gas so that $p(r)=p(\rho(r),\Theta(r)),\mathcal{L}(r)(\rho(r),\Theta(r))$
and $\mathfrak{T}=\mathfrak{T}(r,\Theta(r)$. The four fundamental differential
equations governing the pressure, mass, flux and temperature gradients are
then

$\displaystyle\frac{dp(r)}{dr}=\frac{dp_{gas}(r)}{dr}+\frac{dp_{rad}(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}$
(2.22) $\displaystyle\frac{d\mathcal{M}(r)}{dr}=4\pi\rho(r)r^{2}dr$ (2.23)
$\displaystyle\frac{d\mathcal{L}(r)}{dr}=4\pi r^{2}\mathfrak{T}\rho(r)$ (2.24)
$\displaystyle\frac{d\Theta(r)}{dr}=\frac{3\kappa(r)\rho(r)\mathcal{L}(r)}{4ca|\Theta(r)|^{3}4\pi
r^{2}}$ (2.25)

with BCs $\mathcal{M}(0)=\mathcal{L}(0)=0$ and $\rho(r)=\Theta(r)=0$. The
star’s structure is then determined by four unknowns in term of the variable
$r$, namely $\mathcal{M}(r),\mathcal{L}(r),\rho(r)$ and $\Theta(r)$. However,
in this work we are only concerned with the FEHE equations (2.8)and (2.12)and
their consequences, application and derivation.

The nuclear energy production per mass is expected to obey a power law roughly
of the form

$\displaystyle\mathfrak{T}(\rho,\Theta)\sim\epsilon_{1}|\rho|^{\lambda}(k_{B}\Theta)^{\mu}$
(2.26)

A steady $\mathfrak{T}$ then ensures steady pressure, flux and temperature
gradients within the star. Very light nuclei like
${}^{1}\\!\\!H,^{2}\\!\\!H,^{3}\\!\\!He,^{7}\\!\\!Li$ have less binding energy
per nucleon than medium atomic weight nuclei like nickel or iron. Therefore,
energy can be extracted by the fusion of hydrogen and helium into heavier
elements. While there is hydrogen in the core of the star, the dominant source
of nuclear energy is from the fusion of ${}^{1}\\!\\!H$ into ${}^{4}\\!\\!He$,
which has the largest binding energy of the lighter elements.

For most stars in the Universe, there are two main processes by which
${}^{1}\\!\\!H$ can fuse to ${}^{4}\\!\\!He$: these are the pp-chain and the
CNO-cycle $\bm{[2-5]}$. In yellow main sequence stars like the Sun, the PP-
chain is the main thermonuclear energy source. The main branches are

$\displaystyle\mathrm{I}:~{}^{1}\\!\mathrm{H}+^{1}\\!\mathrm{H}\longrightarrow^{2}\\!\\!\mathrm{H}+e^{+}+\nu_{e}+1.18Mev$
$\displaystyle\mathrm{II}:~{}^{1}\\!\mathrm{H}+^{2}\\!\mathrm{H}\longrightarrow^{3}\\!\\!\mathrm{He}+e^{+}+\gamma+5.49Mev$
$\displaystyle\mathrm{III}:^{3}\\!\mathrm{He}+^{3}\\!\mathrm{He}\longrightarrow^{4}\\!\\!\mathrm{He}+^{1}\\!\mathrm{H}+^{1}\\!\mathrm{H}+12.85Mev$
(2.27)

The pp cycle can be briefly analysed as follows:

  1. (1)

PP-I branch: The cross-section for this reaction is very small for protons
with energy $<1MeV$. This is due to the weak nuclear force which ’transmutes’
a proton into a neutron. The average lifetime of a proton in the core of any
star like the Sun is of the order of $10^{1}0$ years which is comparable to
the age of very old stars. In terms of the inner quark structure of the
proton, the heavy $W^{+}$ boson which carries the weak nuclear force,
transmutes a u-quark in the proton to a d-quark, thus creating a neutron. The
overall cross-section for this reaction is so very small because of the weak
interaction and the suppression factor ( 1/100) of the mutual proton Coulomb
barriers. This very small reaction rate is the reason why stars burn their
nuclear fuel so slowly. The quantum mechanical calculation of probability of
deuteron formation consists of two parts:(i)The penetration of the mutual
Coulomb barriers by two protons in a thermal bath at temperature
$\Theta\geq\Theta_{I}$, where $\Theta_{I}$ is the ’ignition temperature’.(ii)
The probability of inverse $\beta$ decay with positron and neutrino emission.
For this part, Bethe $\bm{[4]}$ applied the original Fermi theory of weak
interactions, which is sufficient for this low energy process.

  2. (2)

PP-II branch: This branch of the cycle is a non-resonant direct capture in a
thermal bath to the ground state of Helium-3, and involves only the strong
nuclear force. This reaction therefore proceeds very rapidly. Other possible
reactions such as $D+D\rightarrow P+^{3}H$,$D+D\rightarrow n+^{3}He$,
$D+^{3}He\rightarrow p+^{4}He$ have larger cross-sections, but the
overwhelming number density of protons in the cores of stellar infernos
ensures that the reaction of deuterium fusing with a proton dominates the $PP-
cycle$. Consequently, deuterium is virtually instantaneously destroyed in all
stars.

  3. (3)

PP-III branch: The pp-chain is complete with the burning of helium-3, creating
a helium-4 nucleus and also 2 protons which become potential new fuel within
the stellar core. Another possible reaction is ${}^{3}He+D\rightarrow^{4}He+p$
but the concentration of deuterium is so small that (PP-II) dominates.
Helium-3 can also be fused with Helium-4 to create ${}^{7}Be$ and the berylium
can in turn capture a proton to form boron-8. These reactions are important in
neutrino astronomy $\bm{[2]}$.

For heavier and hotter (bluer) stars the main nuclear energy production
process is the CNO-cycle involving carbon, nitrogen and oxygen $\bm{[2,3]}$.

$\displaystyle(i):~{}^{1}\\!\mathrm{H}+^{12}\\!\\!\mathrm{C}\longrightarrow^{13}\\!\mathrm{N}+\gamma+1.95Mev$
$\displaystyle(ii):^{13}\mathrm{N}\longrightarrow^{13}\mathrm{C}+e^{+}+\nu_{e}+1.5Mev$
$\displaystyle(iii):^{1}\\!\mathrm{H}+^{13}\\!\\!\mathrm{C}\longrightarrow^{14}\\!\mathrm{N}+\gamma+7.45Mev$
$\displaystyle(vi):^{1}\\!\mathrm{H}+^{14}\\!\\!\mathrm{N}\longrightarrow^{15}\\!\mathrm{O}+\gamma+7.35Mev$
$\displaystyle(vi):^{15}\\!\mathrm{O}\longrightarrow^{15}\\!\mathrm{N}+e^{+}+\nu_{e}+1.73Mev$
$\displaystyle(vii):^{1}\\!\mathrm{H}+^{15}\\!\\!\mathrm{N}\longrightarrow^{4}\\!\mathrm{He}+\gamma+4.96Mev$
(2.28)

The $\mathrm{C,N,O}$ nuclei were created in earlier generations of stars and
act like thermonuclear catalysts in that they are neither created or destroyed
over any complete $\mathrm{CNO}$ cycle. The energies listed here are the total
energies dumped into the stellar material, excluding the neutrino energies
since the neutrinos pass unhindered through the star, but including the
energies from the positron $m_{e}c^{2}$ and the electron with which it
annihilates.

Hydrostatic equilibrium and nuclear energy production are in positive
feedback: HE ensures that the nuclear energy production rate $\mathfrak{T}$,
and hence the luminosity or energy flux $\mathcal{L}(r)$ remains steady over
billions of years. The nuclear energy production in turn maintains the
radiation/photon pressure gradient throughout the star which maintains the
hydrostatic equilibrium. From (2.24) one has

$\displaystyle\frac{d\mathcal{L}(r)}{dr}=4\pi
r^{2}\mathfrak{T}(r)\rho(r)=\frac{d\mathcal{M}(r)}{dr}\mathfrak{T}(r)=\frac{d\mathcal{M}(r)}{dp(r)}\frac{dp(r)}{dr}\mathfrak{T}(r)$
(2.29)

and the temperature, pressure and flux gradients are stable and (to a high
approximation) constant and consistent.

######  Definition 2.5.

For a star in hydrostatic equilibrium, one can estimate the nuclear time, aka
the Einstein time, which is the time required to burn through its nuclear fuel
(hydrogen), concentrated in a central core of mass $M_{core}$. It is
$\bm{[4]}$.

$\displaystyle T_{Nuc}\sim\frac{{\epsilon}M_{core}c^{2}}{\mathcal{L}}$ (2.30)

where $\mathcal{L}$ is the luminosity and $\epsilon$ is an efficiency factor,
typically $\epsilon>0.7$.For the Sun, $\mathcal{L}\sim 3.9\times
10^{9}Js^{-1}$ so that $T_{Nuc}\sim 10^{10}$ years, which is around $10,000$
million years.

####  2.1.2. Central pressure and temperature estimates

######  Definition 2.6.

The Newtonian gravitational potential ${\Phi}(r)$ in the interior of a
spherically symmetric distribution of matter has the gradient

$\displaystyle\nabla{\Phi}(r)=\frac{d\Phi(r)(r)}{dr}=\frac{\mathscr{G}\mathcal{M}(r)}{r^{2}}$
(2.31)

Then the FEHE can be expressed as

$\displaystyle\frac{dp(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}=-\rho(r)\frac{d{\Phi}(r)}{dr}$
(2.32)

Integrating with ${\mathcal{M}}(0)=0$ and $\mathcal{M}(R)=M$ the potential
$\Phi(r)$ for $r\geq R$ is $\Phi(r)=-\frac{\mathscr{G}M}{r}$ and the boundary
potential is $\Phi(R)=-\frac{\mathscr{G}M}{R}$

######  Lemma 2.7.

The FEHE can be expressed in the Poisson equation forms

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{\rho(r)}\frac{dp(r)}{dr}\right)=-4\pi\mathscr{G}\rho(r)$
(2.33)
$\displaystyle\Delta_{r}\Phi(r)=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\Phi(r)}{dr}\right)=-4\pi{\mathscr{G}}\rho(r)$
(2.34)

where (2.32) is the Emden equation.

###### Proof.

This follows easily from (2.8) and (2.12) since

$\displaystyle\frac{r^{2}}{\rho(r)}\frac{dp(r)}{dr}=-{\mathscr{G}}\mathcal{M}(r)=-{\mathscr{G}}\int_{0}^{r}4\pi
r^{\prime 2}\rho(r^{\prime})dr^{\prime}$ (2.35)

Taking the derivative

$\displaystyle\frac{d}{dr}\left(\frac{r^{2}}{\rho(r)}\frac{dp(r)}{dr}\right)=-4\pi{\mathscr{G}}r^{2}\rho(r)$
(2.36)

and (2.32) then follows. Equation (2,31) follows immediately from substituting
$\tfrac{1}{\rho(r)}(dp(r)/dr)=-d\Phi(r)/dr$. ∎

######  Definition 2.8.

The Lagrangian form of the FEHE is derived quite simply. The FEHE is

$\displaystyle\frac{dp(r)}{dr}=\frac{dp(r)}{d\mathcal{M}(r)}\frac{d\mathcal{M}(r)}{dr}=-\frac{{\mathscr{G}}\mathcal{M}(r)\rho(r)}{r^{2}}$
(2.37)

so that

$\displaystyle\mathbf{D}_{\mathcal{M}}p(r)=\frac{dp(r)}{d\mathcal{M}(r)}=-\frac{{\mathscr{G}}\mathcal{M}(r)\rho(r)}{r^{2}}\frac{dr}{d\mathcal{M}(r)}=-\frac{\mathscr{G}G\mathcal{M}(r)\rho(r)}{r^{2}}\frac{1}{4\pi
r^{2}\rho(r)}=-\frac{{\mathscr{G}}\mathcal{M}(r)}{4\pi r^{4}}$ (2.38)

This form is often more useful in stellar structure codes and other
applications.

An immediate consequence of the hydrostatic equilibrium equation in Lagrangian
form, are strict lower bound estimates on the central pressure and central
temperature of the Sun or any main sequence star of mass M and radius R
composed primarily of hydrogen and obeying the universal gas law
$\bm{[1,2,3]}$.

######  Theorem 2.9.

Any gaseous main-sequence star like the Sun of mass $M$ and radius $R$
composed primarily of ionised hydrogen and obeying the universal/perfect gas
law, has the following lower bounds for its central pressure $p_{c}=p(0)$ and
temperature $\Theta_{c}=\Theta(0)$.

$\displaystyle p_{c}\geq\frac{\mathscr{G}M^{2}}{8\pi
R^{4}},~{}~{}~{}\Theta_{c}\geq\frac{\mu
k_{B}}{\mathscr{R}}\frac{\mathscr{G}M}{R}$ (2.39)

where $\mathscr{R}$ is the universal gas constant, $\mu$ is the molecular
weight of hydrogen, and $k_{B}$ is the Boltzmann constant.

###### Proof.

Integrating the Lagrangian form of the FEE (2.38)

$\displaystyle\int_{p(0)}^{p(R)}dp(r)=p(0)=\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{4\pi
r^{2}}=\frac{\mathscr{G}M^{2}}{8\pi R^{4}}$ (2.40)

since $p(R)=0,\mathcal{M}(R)=M$ and $\mathcal{M}(0)=0$. Since
$1/r^{4}>1/R^{4}$ we have the lower bound estimate on the central pressure.

$\displaystyle p_{c}\geq\frac{\mathscr{G}M^{2}}{8\pi R^{4}}$ (2.41)

Since the gas obeys $p(r)=\frac{\mu k_{b}}{\mathscr{R}}\Theta(r)\rho(r)$ so at
the centre of the star $p(0)=\frac{\mu k_{b}}{\mathscr{R}}\Theta(0)\rho(0)$.
Then using (2.41)

$\displaystyle p_{c}=p(0)=\frac{\mu
k_{b}}{\mathscr{R}}\Theta(0)\rho(0)\geq\frac{\mathscr{G}M^{2}}{8\pi R^{4}}$
(2.42)

In terms of the central temperature $\Theta(0)$

$\displaystyle\Theta(0)\geq\frac{\mathscr{G}M^{2}}{8\pi
R^{4}}\left|\frac{1}{\rho(0)}\right|\frac{\mu k_{b}}{\mathscr{R}}\Theta(0)$
(2.43)

The central density can be represented by the average density so that
$\rho(0)\sim\big{\langle}\rho\big{\rangle}=3M/4\pi R^{2}$ and the lower bound
on the central temperature of the star is

$\displaystyle\Theta(0)\geq\frac{3}{8}\frac{GM}{R}\left|\frac{1}{\rho(0)}\right|\frac{\mu
k_{b}}{\mathscr{R}}\Theta(0)$ (2.44)

∎

These are rough estimates/bounds but for the Sun give correct the order of
magnitude for pressure and temperature. The central core temperatures and
pressures of stars are therefore hot enough to initiate thermonuclear fusion
of hydrogen and the radiation/photon pressure which is generated as a result
and maintains equilibrium. It can be see that the central temperature and
pressure in the star will increase if R is decreased or M is increased. Note
also that in the limit as the radius tends to zero.

$\displaystyle\lim_{R\uparrow 0}p_{c})=\lim_{R\uparrow 0}\frac{\mu
k_{b}}{\mathscr{R}}\Theta(0)\rho(0)\geq\frac{\mathscr{G}M^{2}}{8\pi
R^{4}}=\infty$ (2.45)

This does not hold in general relativity since an infinite central pressure
occurs for a critical finite mass and radius of $\mathscr{G}M/R=4/9$. (See
Section 10)).

###  2.2. Hydrostatic equilibrium in higher dimensions

These results extend naturally to n dimensions $\mathbf{R}^{n}$, which is now
very briefly considered However, one can conclude that stars in equilibrium
cannot exist for $n>3$. Let $\mathbf{B}(0,R)\subset\mathbf{R}^{n}$ be an
n-dimensional ball of radius R centered at the origin. The boundary is denoted
$\partial\mathbf{B}(0,R)$. Let $\mathrm{A}(\partial\mathbf{B}(0,1))$ be the
surface are of the unit n-dimensional ball and let
$\mathrm{V}(\mathbf{B}(0,1))$ be the volume. Then the surface area of the
n-ball of radius r is

$\displaystyle\mathrm{A}(\partial\mathbf{B}(0,r))=\mathrm{A}(\partial\mathbf{B}(0,1))r^{n-2}=\frac{2\pi^{n/2}}{\Gamma(n/2)}r^{n-2}$
(2.46)

++ and the volume is

$\displaystyle\mathrm{V}(\mathbf{B}(0,r))=\mathrm{V}(\mathbf{B}(0,1))r^{n}=\frac{2\pi^{n/2}}{\Gamma(n/2+1)}r^{n}$
(2.47)

where $\Gamma$is the standard gamma function. If $\mathbf{B}(0,R)$ supports a
gaseous star of density $\rho(r)$ and pressure $p(r)$ then the FEHE is

$\displaystyle\frac{dp(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{n-1}}$
(2.48)

and the mass $\mathcal{M}(r)$ contained within a radius $r$ is

$\displaystyle\mathcal{M}(r)=\int_{0}^{r}|\mathrm{A}(\partial\mathbf{B}(0,1))|r^{n-1}dr=\int_{0}^{r}\frac{2\pi^{n/2}}{\Gamma(n/2)}r^{n-2}dr$
(2.49)

The n-dimensional generalisation is

$\displaystyle\frac{1}{r^{n-1}}\frac{d}{dr}\left(\frac{r^{n-1}}{\rho(r)}\frac{dp(r)}{dr}\right)=-\mathrm{A}(\partial\mathbf{B}(0,1))\mathscr{G}\rho(r)$
(2.50)

The Lagrangian form of (-) for $n>3$ is then

$\displaystyle\frac{dp(r)}{d\mathcal{M}(r)}=\frac{\mathscr{G}\mathcal{M}(r)}{\mathrm{A}(\mathbf{B}(0,1))r^{2(n-1)}}$
(2.51)

The central temperature and pressure lower bounds can be estimated as before.
However, the central and pressure and temperature estimates are now
drastically lowered when $n>3$ for a star of the same mass and radius.

######  Theorem 2.10.

For $n>3$, the central pressure and temperature lower bounds can be estimated
as

$\displaystyle
p_{c}\geq\frac{\mathscr{G}M^{2}}{2\mathcal{A}(\mathbf{B}(0,1)R^{2(n-1)}}=\frac{\mathscr{G}M^{2}\Gamma(n/2)}{4\pi^{n/2}R^{2(n-1)}}$
(2.52)
$\displaystyle\Theta_{c}\geq\frac{\mathscr{R}}{k_{B}\mu}\frac{\Gamma(n/2)}{\Gamma(\tfrac{n}{2}+1)}\frac{\mathscr{G}M}{R^{(n-2)}}$
(2.53)

###### Proof.

The central pressure is estimated as before via integrating (2.48) so that

$\displaystyle
p(0)\geq\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{\mathrm{A})\mathbf{B}(0,1)r^{2(n-1)}}=\frac{\mathscr{G}M^{2}}{2\mathrm{A}(\mathbf{B}(0,1))R^{2(n-1)}}$
(2.54)

To estimate the central temperature bound, apply the perfect gas law as before
then

$\displaystyle\frac{k_{B}\mu}{\mathscr{R}}\Theta(0)\rho(0)=p(0)\geq\frac{\mathscr{G}M^{2}}{2\mathrm{A}(\mathbf{B}(0,1))R^{2(n-1)}}=\frac{\mathscr{G}M^{2}\Gamma(n/2)}{4\pi^{n/2}R^{3n-4}}$
(2.55)

Using the mean density for $\rho(0)$

$\displaystyle\rho(0)\sim\big{\langle}\rho\big{\rangle}=\frac{M}{\mathrm{V}(\mathbf{B}(0,R))}=\frac{M\Gamma(\frac{n}{2}+1)}{2\pi^{n/2}R^{n}}$
(2.56)

and substituting in (2.55), the result (2.52) then follows. ∎

It can be observed that for $n>$ the central temperature and pressure in the
star are suppressed for the same mass and radius. In the limit as
$n\rightarrow\infty$ the temperature and pressure are zero. This suggests that
for $n>3$ and especially $n\gg 3$, stars are unstable and cannot achieve
equilibrium or ignition temperature for thermonuclear fusion of hydrogen.

####  2.2.1. Some basic theorems concerning hydrostatic equilibrium

Returning to $n=3$, a star in hydrostatic equilibrium also obeys a number of
basic theorems and presented in Chandrasekhar in $\bm{[1]}$ with proofs given
in $\bm{[1]}$.

######  Theorem 2.11.

If the average density
$\langle\rho(r)\rangle\sim\mathcal{M}(r)/[(4\pi/3)r^{3}]$ is monotonically
non-increasing as one moves outward in the star towards the surface, then in
any equilibrium configuration, the function $\mathcal{F}(r)$ given by

$\displaystyle\mathcal{F}(r)=p(r)+\frac{3}{8\pi}\frac{|\mathcal{M}(r)|^{2}}{r^{4}}$
(2.57)

is also non-increasing.

###### Proof.

$\displaystyle\frac{dp(r)}{dr}=\frac{d}{dr}[p(r)+\frac{3}{8\pi}\frac{|\mathcal{M}(r)|^{2}}{r^{4}}]$
(2.58)
$\displaystyle=\underbrace{\frac{dp(r)}{dr}}_{use~{}Lagrange~{}form~{}of~{}FEHE}+\frac{3}{4\pi}\mathscr{G}\frac{\mathcal{M}(r)\mathcal{M}^{\prime}(r)}{r^{4}}-\frac{3}{2\pi}\frac{|\mathcal{M}(r)|^{2}}{r^{5}}$
$\displaystyle=-\frac{\mathscr{G}\mathcal{M}^{\prime}(r)\mathcal{M}(r)}{4\pi
r^{4}}+\frac{3}{4\pi}\frac{\mathscr{G}\mathcal{M}^{\prime}(r)\mathcal{M}(r)}{4\pi
r^{4}}-\frac{3}{2\pi}\frac{|\mathcal{M}(r)|^{2}}{r^{5}}$
$\displaystyle=\frac{1}{2\pi}\frac{\mathscr{G}\mathcal{M}(r)\mathcal{M}^{\prime}(r)}{r^{4}}-\frac{3}{2\pi}\frac{|\mathcal{M}(r)|^{2}}{r^{5}}$
$\displaystyle=\frac{1}{2\pi}\frac{\mathscr{G}\mathcal{M}(r)}{r^{4}}\left(\mathcal{M}^{\prime}(r)-\frac{3\mathcal{M}(r)}{r}\right)$
$\displaystyle=\frac{1}{2\pi}\frac{\mathscr{G}\mathcal{M}(r)}{r}\frac{d}{dr}\left|\frac{\mathcal{M}(r)}{r^{3}}\right|=\frac{2}{3}\mathscr{G}\frac{\mathcal{M}(r)|}{r}\frac{d}{dr}[\langle\rho(r)\rangle]\leq
0$ (2.59)

∎

An immediate corollary is a lower bound on the central pressure in the star

######  Corollary 2.12.

$\displaystyle p(0)=p_{c}\geq
p(r)+\frac{3}{8\pi}\frac{\mathscr{G}|\mathcal{M}(r)|^{2}}{r^{4}}\geq\frac{3\mathscr{G}}{8\pi}\left|\frac{M^{2}}{R^{4}}\right|$
(2.60)

which is consistent with (2.41).

######  Theorem 2.13.

For any equilibrium configuration

$\displaystyle\mathbf{I}_{\gamma}=\int_{o}^{R}\frac{{\mathscr{G}}\mathcal{M}(r)d\mathcal{M}(r)}{r^{\xi}}=4\pi(4-\xi)\int_{0}^{R}p(r)r^{3-\xi}d\xi$
(2.61)

if $\xi<4$.

######  Theorem 2.14.

$\displaystyle\xi\pi
p_{c}R^{4}r^{4-\xi}+\frac{|4-\xi|}{8}\frac{\mathscr{G}M^{2}}{R^{3}}>\mathbf{I}_{\xi}>\frac{\mathscr{G}M^{2}}{2R^{\xi}}$
(2.62)

if $\xi<4$.

###  2.3. Gravitational and Thermal Energies of a Star in Hydrostatic
Equilibrium

The following energies can be defined for a self-gravitating star. Because a
star is a self-gravitating thermodynamic system this leads to some interesting
properties; in particular, the fact that a star has negative specific heat: as
it loses energy (via radiation) it can get hotter. Gravitational energy can
also be converted to thermal energy; for example, as a protostar radiates and
collapses in a (Kelvin-Helmholtz contraction) it gets hotter until
thermonuclear reactions are initiated in its core. The star achieves
hydrostatic equilibrium and becomes stable. The most commonly utilised
energies in Newtonian astrophysical theory are the gravitational and thermal
energies which are now defined.

######  Definition 2.15.

The total gravitational potential energy (TGPE) or binding energy
$\mathbf{E}_{G}$ of a star is the energy required to disperse the matter
comprising the star to infinity, peeling it shell by shell from the outside in
$\bm{[1,2,3]}$.

$\displaystyle\mathbf{E}_{G}=-4\pi\mathscr{G}\int_{0}^{R}r\mathcal{M}(r)\rho(r)dr\equiv\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r}$
(2.63)

The TGPE can also be written as

$\displaystyle\mathbf{E}_{G}=-3\int_{0}^{R}4\pi r^{2}p(r)dr$ (2.64)

This follows by using the FEHE and integrating by parts so that

$\displaystyle\mathrm{E}_{G}=4\pi\int_{0}^{R}\frac{dp(r)}{dr}r^{3}dr=[r^{3}p(r)]_{0}^{R}-3\int_{0}^{R}4\pi
r^{2}p(r)r^{2}dr$ (2.65)

and the first term on the RHS vanishes since $p(R)=0$.

######  Definition 2.16.

The thermal energy of the star is defined as $\bm{[2,3]}$.

$\displaystyle\mathbf{E}_{T}=\int_{0}^{R}\mathcal{E}(r)4\pi
r^{2}dr=\int_{0}^{R}\mathcal{E}(r)d\mathrm{V}(r)$ (2.66)

where $\mathcal{E}(r)$ is the density of the internal thermal energy only, not
including the gravitational energy or self gravitation and the rest-mass
energy. For example, a perfect/ideal gas of N particles at temperature
$\Theta$ contained in a domain $\mathbf{B}\subset\mathbf{R}^{3}$ of volume
$\mathrm{V}(\mathbf{B})$, the homogenous energy density is
$\mathcal{E}=\tfrac{3}{2}Nk_{B}\Theta/\mathrm{V}=E_{\Theta}/\mathrm{V}(\mathbf{B})$
so that $\mathbf{E}_{\Theta}=\mathrm{V}(\mathbf{B})\mathcal{E}=\tfrac{4}{3}\pi
R^{3}\mathcal{E}$. For a gas in spherical ’box’ the thermal energy is
$\mathbf{E}_{\Theta}=\int_{\mathbf{B}}\mathcal{E}(r)d\mathrm{V}=\int_{\mathbf{B}}\mathcal{E}(r)4\pi
r^{2}dr$. The thermal energy can also be expressed as integrals over
$d\mathcal{M}(r)$ so that

$\displaystyle\mathbf{E}_{T}=\int_{0}^{R}\mathcal{E}(r)4\pi
r^{2}dr=\int_{0}^{R}\frac{\mathcal{E}(r)}{\rho(r)}d\mathcal{M}(r)=\int_{0}^{M}\mathscr{J}(r,\rho(r))d\mathcal{M}(r)$
(2.67)

where $\mathscr{J}$ is the thermal energy density per unit mass.

######  Definition 2.17.

The total energy $E$ of the star is then the sum of the thermal and
gravitational energies

$\displaystyle\mathbf{E}=\mathbf{E}_{T}+\mathbf{E}_{G}=\int_{0}^{R}[\mathcal{E}(r)-2p(r)]4\pi^{2}dr\equiv\int_{0}^{R}[\mathcal{E}(r)-2p(r)]d\mathrm{V}(r)$
(2.68)

The star must have negative energy and is stable against thermal
dispersion/expansion of its matter if $\mathcal{E}(r)<3p(r)$.

######  Remark 2.18.

The following basic (heuristic) observations can be made regarding the total
energy.

  * •

If $\mathbf{E}<0$ with $\mathbf{E}(r)<3p(r)$ then the star is stable to
dispersion/explosion.

  * •

If $\mathbf{E}>0$ then the star will disperse/explode or collapse.

  * •

A star with $\mathbf{E}=0$ is at the point of dynamical instability.

In Section (4) it will be shown that the vanishing of the 1st variation of E
(namely $\delta\mathbf{E}=0$) leads to the FEHE.

###  2.4. The stellar virial theorem

The virial theorem, well known in classical and celestial mechanics is another
consequence of the FEHE and the gravitational potential/binding energy. The
virial theorem for self-gravitating gaseous stars follows quite simply from
the Lagrangian form of the FEHE and the gravitational potential energy
$\mathbf{E}_{G}$.

######  Theorem 2.19.

The virial theorem for a self-gravitating gaseous star is

$\displaystyle
3\int_{0}^{\mathrm{V}}p(r)d\mathrm{V}(r)=\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r}$
(2.69)

or

$\displaystyle
3\int_{0}^{\mathrm{V}}p(r)d\mathrm{V}(r)+{\mathbf{E}}_{\scriptstyle{G}}=0$
(2.70)

###### Proof.

Using the Lagrangian form of the FEHE

$\displaystyle\frac{dp(r)}{d\mathcal{M}(r)}=-\frac{\mathscr{G}\mathcal{M}(r)}{4\pi
r^{4}}$ (2.71)

multiply through by $4\pi r^{2}$ so that

$\displaystyle 4\pi
r^{3}\frac{dp(r)}{d\mathcal{M}(r)}=-\frac{\mathscr{G}\mathcal{M}(r)}{r}$
(2.72)

and

$\displaystyle 4\pi
r^{3}dp(r)=-\frac{\mathscr{G}\mathcal{M}(r)}{r}d\mathcal{M}(r)$ (2.73)

Integrating over the volume of the star with $\mathrm{V}(r)=\tfrac{4}{3}r^{3}$
and volume element $d\mathrm{V}=4\pi r^{2}dr$

$\displaystyle\underbrace{3\int_{p(0)}^{p(R)}\mathrm{V}(r)dp(r)}_{int~{}by~{}parts}=-\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)}{r}d\mathcal{M}(r)$
(2.74)
$\displaystyle|p(r)\mathrm{V}(r)|_{p(0)}^{p(R)}-3\int_{\mathrm{V}(0)}^{\mathrm{V}(R)}p(r)d\mathrm{V}(r)=-\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r}$
(2.75)

and $\mathrm{V}(0)=0$ and $p(R)=0$. Hence (2.69) or (2.70) follows. ∎

The virial theorem can also be expressed in terms of the thermal energy
$\mathbf{E}_{T}$ and gravitational energy $\mathbf{E}_{G}$. An immediate
corollary is that the specific heat of a gaseous star at fixed volume is
actually negative.

######  Theorem 2.20.

The virial theorem can be expressed in terms of the thermal and gravitational
energies as

$\displaystyle 2{\mathbf{E}}_{T}=-{\mathbf{E}}_{G}$ (2.76)

so that the total energy is minus the thermal energy

$\displaystyle\mathbf{E}=\mathbf{E}_{T}+\mathbf{E}_{G}=-\mathbf{E}_{T}$ (2.77)

The specific heat is then negative so $\mathrm{C}_{V}<0$. For a perfect gas

$\displaystyle\mathrm{C}_{V}=-\tfrac{3}{2}Nk_{b}<0$ (2.78)

###### Proof.

Since the pressure $p(r)$ and the thermal energy density $\mathcal{E}(r)$ are
related by $3p(r)=2\mathcal{E}$, the virial theorem becomes

$\displaystyle 2\int_{0}^{R}{\mathcal{E}}(r)4\pi
r^{2}dr=2{\mathbf{E}_{T}}=-{\mathbf{E}_{G}}$ (2.79)

The total energy $\mathbf{E}$ is then

$\displaystyle{\mathbf{E}}={\mathbf{E}_{T}}+{\mathbf{E}_{G}}=-\mathbf{E}_{T}$
(2.80)

The specific heat in hydrostatic equilibrium at fixed $\mathrm{V}$ is then the
derivative of $\mathbf{E}$ with respect to temperature $\Theta$ so that

$\displaystyle{\mathrm{C}}_{\mathrm{V}}=\frac{d\mathbf{E}}{d\Theta}=-\frac{d\mathbf{E}_{T}}{d\Theta}<0$
(2.81)

For an ideal gas at temperature $\Theta$ one has
$\mathbf{E}_{T}=\frac{3}{2}Nk_{b}\Theta$ so
that$\mathrm{C}_{\mathrm{V}}=-\frac{3}{2}Nk_{b}<0$. ∎

Physically, this means that as a self-gravitating system loses energy by
radiation it gravitationally contracts since it cannot be supported in
hydrostatic equilibrium. Therefore it heats up. During the lifetime of the
star, the radiation pressure is maintained by thermonuclear reactions in the
core. But once this reservoir of fuel becomes exhausted the star will shift
out of hydrostatic equilibrium and contract, and gravitational energy will be
converted to thermal energy. In a main sequence star like the Sun, the
increased core temperature initiates the next phase of thermonuclear burning
whereby helium is fused into carbon and oxygen. The temperature of ordinary
matter with positive specific heat, decreases as the matter radiates energy;
for example, a red-hot piece of iron will cool off. The statistical mechanics
and thermodynamics of self gravitating systems of particles has been discussed
in various works(refs).

###  2.5. The isothermal sphere for a perfect gas

The equilibrium of an isothermal gas sphere (IGS) is now considered for a
perfect gas. This is a self-gravitating gaseous sphere of uniform temperature
$\bm{[1,6]}$. Real stars have thermal gradients so this is not a realistic
model of any normal star. However, once a star has exhausted its nuclear fuel,
the central core of the star–consisting primarily of helium for main sequence
stars–will become isothermal, as a prelude to becoming a red giant. It can
also describe globular clusters. Thus, the isothermal core is still astro-
physically relevant. Also, the isothermal sphere has a number of interesting
thermodynamic and statistical mechanical properties. $\bm{[1,6,7,8]}$.

######  Lemma 2.21.

Let $\mathbf{B}(0,R)\subset\mathbf{R}^{3}$ be a ball/domain of radius $R$
centred at the origin and supporting a perfect gas at uniform temperature
$\Theta$ consisting of N particles of mass $m$. The total mass is then $M=mN$.
The pressure and density at any radius $r<R$ are then related by
$p(r)=\mathscr{R}\rho(r)=\frac{k_{B}\Theta}{m}\rho(r)$. The equilibrium
equations for the isothermal gas sphere are then the Poisson and Poisson-
Boltzmann equations

$\displaystyle\Delta(\log\rho(r))=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d}{dr}(\log\rho(r))\right)=-\frac{4\pi\mathscr{G}\rho(r)}{\mathcal{R}}$
(2.82)
$\displaystyle\Delta\Phi(r)=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\Phi(r)}{dr}\right)=4\pi\mathscr{G}\rho_{c}\exp\left(-\frac{\Phi(r)}{\mathscr{R}}\right)$
(2.83)

###### Proof.

Differentiating the perfect gas law

$\displaystyle\frac{dp(r)}{dr}=\mathscr{R}\frac{d\rho(r)}{dr}$ (2.84)

and applying the FEHE $dp(r)/dr=-\frac{\mathscr{G}\mathcal{M}(r)}{r^{2}}$
gives

$\displaystyle\frac{dp(r)}{dr}=\mathscr{R}\frac{d\rho(r)}{dr}=-\rho(r)\frac{d\Phi(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}$
(2.85)

Integrating

$\displaystyle\mathscr{R}\int_{\rho(0)}^{\rho(r)}\frac{d\overline{\rho(r)}}{\overline{\rho(r)}}=-\int_{\Phi(0)}^{\Phi(r)}d\overline{\Phi(r)}=-\Phi(r)$
(2.86)

since $\Phi(0)=0$. Then

$\displaystyle
p(r)=\log\left|\frac{\rho(r)}{\rho(0)}|\right|=-\frac{\Phi(r)}{\mathscr{R}}$
(2.87)

so that

$\displaystyle\rho(r)=\rho_{c}\exp\left(-\frac{\Phi(r)}{\mathscr{R}}\right)$
(2.88)

To derive (2.82), begin with the Emden equation

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{\rho(r)}\frac{dp(r)}{dr}\right)=-4\pi\mathscr{G}\rho(r)$
(2.89)

Then from the FEHE and (2.87)

$\displaystyle\frac{1}{\rho(r)}\frac{dp(r)}{dr}=\frac{-\mathscr{G}\mathcal{M}(r)}{r^{2}}=-\frac{d\Phi(r)}{dr}=\mathscr{R}\frac{d}{dr}\log\left|\frac{\rho(r)}{\rho(0)}\right|$
(2.90)
$\displaystyle=\mathscr{R}\frac{d}{dr}|\log\rho(r)-\log\rho(0)|=\mathscr{R}\frac{d}{dr}\log\rho(r)$
(2.91)

so that the Emden equation becomes

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d}{dr}(\log\rho(r))\right)=-\frac{4\pi\mathscr{G}\rho(r)}{\mathcal{R}}$
(2.92)

∎

######  Corollary 2.22.

The IGE has a singular solution

$\displaystyle\Phi(r)=\frac{1}{2\pi\mathscr{R}\mathscr{G}r^{2}}$ (2.93)

Since $\rho(r)\sim r^{-2}$ the total mass is $M=\int_{0}^{\infty}4\pi
r^{2}dr=\int_{0}^{\infty}4\pi dr=\infty$

####  2.5.1. Relevance of the isothermal gas sphere to red giant evolution

Once the core of a main-sequence star has exhausted its supply of hydrogen
fuel, it gradually ceases to produce nuclear energy and, in the limit of
thermal equilibrium, becomes isothermal. To see this, consider equations
(2.24) and (2.25). Once thermonuclear reactions begin to cease within the core
region then the thermonuclear energy production rate per mass diminishes so
that $\mathfrak{N}\rightarrow 0$ and

$\displaystyle\frac{d\mathcal{L}(r)}{dr}=4\pi r^{2}\rho(r)\mathfrak{T}=0$
(2.94)

Equation (-) for the temperature gradient then becomes

$\displaystyle\frac{d\Theta(r)}{dr}\sim
C\frac{1}{|\Theta(r)|^{3}}\frac{\mathcal{L}(r)}{4\pi r^{2}}=0$ (2.95)

so that $\Theta\rightarrow const.$ and the core region becomes isothermal;
essentially an isothermal gas sphere of $He^{4}$. Schonberg and Chandrasekhar
(1942) showed that, if the envelope is polytropic with index $n=3$, then there
is a maximum fractional mass that the core can achieve. If the core is less
massive, it can remain isothermal while nuclear reactions continue in a
surrounding shell, and the isothermal core temperature is set by the
temperature of this shell. If the SC mass is exceeded then the core contracts
until it is supported by electron degeneracy pressure or helium begins to
burn/fuse at the centre (the triple-alpha reaction) to create carbon and
oxygen. The idealised result is sufficiently accurate that it has become a
well-established element of the theory of the evolution of stars once they
move off the main sequence. It is referred to simply as the
Schonberg–Chandrasekhar (SC) limit $\bm{[9]}$.

######  Definition 2.23.

The Schonberg-Chandrasekhar limit is the maximum fraction of a star’s mass
that can exist as an isothermal core and still support the outer mass envelope
in a state of hydrostatic equilibrium. It is

$\displaystyle\frac{M_{CORE}}{M}\sim 0.37|\frac{\mu_{ENV}}{\mu_{CORE}}$ (2.96)

where $\mu_{ENV}$ and $\mu_{CORE}$ are the molecular weights in the envelope
and core regions respectively.

Similar limits have been computed for other polytropic solutions. Beech (1988)
calculated the corresponding limit for an isothermal core surrounded by an
envelope with n = 1. Eggleton, Faulkner and Cannon (1998) discovered that,
when n = 1 in the envelope and n = 5 in the core, a fractional mass limit
exists if the density decreases discontinuously at the core-envelope boundary
by a factor exceeding 3. They went further to propose conditions on the
polytropic indices of the core and envelope that lead to fractional mass
limits. We refer to all these limits, including the original result of
Schonberg and Chandrasekhar (1942) as S-C-like limits.

##  3\. Deviations From Hydrostatic Equilibrium and Stellar Timescales

So far, only stars in hydrostatic equilibrium have been considered. Suppose
now the star deviates from hydrostastic equilibrium via radial perturbations
or via a decrease/increase in the interior pressure. If gravitation dominates
the star will collapse; if internal/thermal forces dominate the star will
dissipate or explode. This deviation from equilibrium occurs only towards the
end of the star’s life when it has depleted its nuclear fuel. A more formal
treatment of perturbations away from equilibrium is given in Section (4).
Here, a heuristic result can be derived from the FEHE $\bm{[3]}$.

$\displaystyle\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}+\frac{1}{\rho(r)}\frac{dp(r)}{dr}=0$
(3.1)

If the star deviates slightly from hydrostatic equilibrium then

$\displaystyle\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}+\frac{1}{\rho(r)}\frac{dp(r)}{dr}~{}{\gtrless}~{}~{}0$
(3.2)

This induces an acceleration of the mass $\mathcal{M}(r)$ either inward or
outward so that

$\displaystyle\frac{d^{2}\lambda(t)}{dt^{2}}=\frac{\mathscr{G}\mathcal{M}(\lambda(t))\rho(r)}{\lambda(t)^{2}}+\frac{1}{\rho(r)}\frac{dp(r)}{d\lambda(t)}~{}{\gtrless}~{}~{}0$
(3.3)

where $\lambda(t)$ is now a scale factor or time-dependent radius and
$\mathcal{M}(\lambda(t))$ is the mass contained within a (moving) radius
$\lambda(t)$ at time $t$. This can be written as

$\displaystyle\frac{d^{2}\lambda(t)}{dt^{2}}=\frac{r\mathcal{M}(\lambda(t))\rho(r)}{|\lambda(t)|^{2}}\left(1+\frac{\frac{1}{\rho(r)}\frac{dp(r)}{d\lambda(t)}}{\frac{\mathscr{G}\mathcal{M}(\lambda(t))\rho(r)}{\lambda(t)^{2}}}\right)~{}{\gtrless}~{}~{}0\
$ (3.4)

###  3.1. Zero-pressure gravitational collapse and the hydrodynamical time

Suppose now the pressure is ’switched off’ entirely so that $p(r)=0$ and the
pressure gradient vanishes so that $dp(r)/dr=0$. Then (3.4) reduces to

$\displaystyle\frac{d^{2}\lambda(t)}{dt^{2}}=\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{|\lambda(t)|^{2}}<0,~{}~{}\lambda(t)\leq
R$ (3.5)

Introducing a velocity $\mathrm{U}(r(t))=d\lambda(t)/dt$

$\displaystyle\frac{d\mathrm{U}(t)}{dt}=\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{\lambda(t)^{2}}<0$
(3.6)

Pressureless matter or matter with $p\ll 1$ is often referred to as ’dust’ in
the general relativity literature. Within the Newtonian theory, the time for
the sphere of pressureless dust to collapse to zero, is the free fall time or
the hydrodynamic time. With $p(r)=0$ within the star, it undergoes
gravitational collapse to zero size or infinite density within a time $T_{H}$,
the hydrodynamical time.

######  Lemma 3.1.

Let $\mathbf{B}\subset\mathbf{R}^{3}$ be a domain supporting a gaseous star of
radius R and mass M at time t=0. The initial density is then $\rho(0)=3M/4\pi
R^{3}$ and $\mathrm{V}(\mathbf{B}(0,R))=\tfrac{4}{3}\pi R^{3}$. For $t>0$ the
star is no longer in hydrostatic equilibrium so that the support
$\mathbf{B}(0,\lambda(t))$ evolves (shrinks) for increasing t with
$\mathbf{B}(0,\lambda(t))\subset\mathbf{B}(0,R)$ and
$\mathrm{V}(0,\lambda(t))=\tfrac{4}{3}\pi|\lambda(t)|^{3}$. Since the star is
collapsing for $t>0$ then at a sequence of times
$t_{1}<t_{2}<...<t_{n-1},t_{n}$.

$\displaystyle\mathbf{B}(0,\lambda(t_{n})\subset\mathbf{B}(0,\lambda(t_{n-2})\subset...\subset\mathbf{B}(0,\lambda(t_{2}))\subset\mathbf{B}(0,\lambda(t_{1}))$

Then the star collapses to zero size at the hydrodynamical time $T_{H}$ which
is

$\displaystyle T_{H}=\left(\frac{3\pi}{32\mathscr{G}\rho(0)}\right)^{1/2}$
(3.7)

Then $\mathrm{V}(\mathbf{B}(\lambda(T_{H}))=0$ and
$\rho(0,\lambda(t_{H}))=3M/4\pi|\lambda(t)|^{3}=\infty$

###### Proof.

The ODE (3.6) is expressible as

$\displaystyle\frac{d^{2}\lambda(t)}{dt^{2}}=\frac{d}{dt}\left(\frac{d\lambda(t)}{dt}\right)=\frac{d}{dt}\mathrm{U}(\lambda(t))=\frac{d\lambda(t)}{dt}\frac{d}{d\lambda(t)}\mathrm{U}(\lambda(t))$
$\displaystyle=\mathrm{U}(\lambda(t))\frac{d}{d\lambda(t)}\mathrm{U}(\lambda(t))=\frac{1}{2}\frac{d}{d\lambda(t)}|\mathrm{U}(\lambda(t))|^{2}=\frac{\mathscr{G}\mathcal{M}(r)}{\lambda(t)^{2}}=-\frac{\mathscr{G}M^{2}}{|\lambda(t)|^{2}}$
(3.8)

Using the boundary conditions
$\mathcal{M}(0)=0,\mathcal{M}(R)=M,\rho(R)=0,\mathrm{U}(\lambda(0))=0$ and
integrating gives

$\displaystyle\int_{0}^{\mathrm{U}(\lambda(t))}d|\mathrm{U}(\lambda(t))|^{2}=-2\int_{R}^{r}\frac{\mathscr{G}M}{\overline{r}^{2}}d\overline{r}=2\mathscr{G}M\left(\frac{1}{\lambda(t)}-\frac{1}{R}\right)$
(3.9)

where the radial integral is inward since the sphere is collapsing. Then

$\displaystyle\mathrm{U}(\lambda(t))=\frac{d\lambda(t)}{dt}=\pm
2\mathscr{G}M\left(\frac{1}{\lambda(t)}-\frac{1}{R}\right)^{1/2}$ (3.10)

Taking the negative sign and integrating

$\displaystyle\int_{R}^{\lambda(t)}\frac{d\overline{\lambda(t)}}{\left(\frac{1}{\overline{\lambda(t)}}-\frac{1}{R}\right)^{1/2}}=2\mathscr{G}M\int_{0}^{t}dt=-2\mathscr{G}Mt$
(3.11)

so that‘

$\displaystyle\left|R\left(\frac{1}{\overline{\lambda(t)}}-\frac{1}{R}\right)^{1/2}r+R^{3/2}\tan^{-1}\left(\sqrt{R}\left(\frac{1}{\overline{\lambda(t)}}-\frac{1}{R}\right)\right)\right|_{R}^{\lambda(t)}=\left(2\mathscr{G}M\right)^{1/2}t$
(3.12)

The lhs vanishes for $\lambda(t)=\lambda(0)=R$ so that

$\displaystyle
R\left(\frac{1}{\overline{\lambda(t)}}-\frac{1}{R}\right)^{1/2}r+R^{3/2}\tan^{-1}\left(\sqrt{R}\left(\frac{1}{\overline{\lambda(t)}}-\frac{1}{R}\right)\right)=\left(2\mathscr{G}M\right)^{1/2}t$
(3.13)

Since $\tan^{-1}(0)=0,\tan^{-1}(\infty)=\frac{\pi}{2}$ then at
$r=\lambda(t)=\lambda(T_{H})=0$

$\displaystyle R^{3/2}\frac{\pi}{2}=(2\mathscr{G}M)^{1/2}T_{H}$ (3.14)

and squaring

$\displaystyle R^{3}\frac{\pi}{2}=(2\mathscr{G}M)T_{H}$ (3.15)

Since the initial density is $\rho(0)=3M/4\pi R^{3}$ then

$\displaystyle\frac{8\pi\mathscr{G}\rho_{o}T_{H}^{2}}{3}=\frac{\pi^{2}}{4}$
(3.16)

Solving for $T_{H}$ then gives (3.7) as required. ∎

The exact same result also follows from the Tolman-Oppenheimer-Syner models
which apply an Einstein-matter system to pressureless gravitational
collapse.(See $\bm{[10]}$. The star collapses to zero size within a finite
comoving proper time and essentially forms a black hole. However, from the
Newtonian theory one can still deduce that a sufficiently heavy star with
negligible pressure could collapse out of existence within a finite time, once
it is no longer in HE. However, main-sequence Newtonian stars like the Sun and
stars of comparable mass, will collapse to form white dwarfs since the
pressure is always finite and never exactly zero.

###  3.2. The Kelvin-Helmholtz time scale

A self-gravitating gaseous body can convert gravitational energy into internal
thermal energy by a process of contraction $\bm{[2,3]}$. From (2.39) it was
seen that the central pressure and temperature increase as the radius
decreases. This is also consistent with the fact that self-gravitating gases
have negative specific heats so the star will heat up internally as it
radiates away energy and contracts. The total energy of the star is the sum of
the thermal and gravitational energies so that
$\mathbf{E}=\mathbf{E}_{T}+\mathbf{E}_{G}<0$. The Kelvin-Helmhotz time
$T_{KH}$ for a self-gravitating gaseous body of mass M and initial radius R,
to convert its gravitational potential energy to thermal energy, at a constant
luminosity $\mathcal{L}$. Suppose a star is powered only by the process of
conversion of gravitational to thermal energy then $T_{KH}$ is given by

$\displaystyle
T_{KH}=-\frac{\mathbf{E}_{G}}{\mathcal{L}}=\frac{1}{\mathcal{L}}\int_{\mathcal{M}(0)}^{\mathcal{M}(R)}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r}=\frac{\mathscr{G}M^{2}}{2R\mathcal{L}}$
(3.17)

For the Sun, $\mathcal{L}\sim 3.9\times 10^{26}W$ so that $T_{KH}\sim 30\times
10^{6}$ years.
Historical remark
For the Sun, $T_{KH}\sim 30\times 10^{6}$ years, which is much less than the
age of the Sun or the Earth. Also, the Sun would have greatly contracted in
size over this timescale if this was its only energy source. In the 19th
century this theory was proposed by Lord Kelvin (William Thompson) and
Helmholtz (and also Sitter) as an explanation as to how the Sun and stars
could produce their energy. At this time, enough was known about gravitation,
gas dynamics, hydrodynamics and thermodynamics that physically viable models
of the Sun and stars could be formulated. But even at the time of Kelvin,
there was geological and fossil evidence to support the age of the Earth being
far greater than the estimate $T_{KH}$ of Kelvin. Two years earlier, Darwin
had also published his work _Origin of Species_ which suggested that life on
Earth took billions of years of evolve, so the Sun had to be at least as old.
Kelvin also abandoned his theory because he learned that Alexander the Great
had observed a total solar eclipse when he crossed the River Oxus in 329 BC.
This put an upper limit on the size of the Sun at that date, suggesting that
the Sun was not contracting, or at least not contracting fast enough, to
provide the required energy output. It would be many decades until the
discovery of radioactivity and also quantum mechanics which would (in the 20th
century) lead to a full understanding of thermonuclear reactions and
processes; the actual source of the immense heat of the stars that maintains
them in perpetual equilibrium with constant luminosity. However, Kelvin-
Helmholtz contractions still describes collapsing protostars which form in gas
clouds and which radiate thermally by this process in the phase before nuclear
reactions are initiated in their core.

##  4\. Derivation Of The Hydrostatic Equilibrium Equation Via Variational
Methods

The hydrostatic equilibrium criteria are now examined from a more formal
perspective via variational methods. The HEE for the isothermal gaseous sphere
can also be derived via a variation method from a statistical mechanical
distribution function. It is seen in both cases that the equilibria correspond
to the critical points. First, the definitions of 1st and 2nd variations of a
function are briefly reviewed.

######  Definition 4.1.

Let $\Theta=\Theta(f(r))$ be a functional of a function $f(r)$ with
$f:\mathbf{R}\rightarrow\mathbf{R}$ and
$\Theta:\mathbf{R}\rightarrow\mathbf{R}$. For spherical symmetry $r\leq R$.
The 1st variation with respect to $\delta f(r)$ is then

$\displaystyle\delta\Lambda(f(r))=\frac{\partial\Lambda(f(r))}{\partial
f(r)}\delta f(r)=\frac{\partial\Lambda(f(r))}{\partial
f(r)}\frac{df(r)}{dr}\delta r$ (4.1)

The second variation is

$\displaystyle\delta^{2}\Lambda(f(r))=\frac{\partial\Lambda(f(r))}{\partial
f(r)}\delta f(r)+\frac{1}{2}\frac{\partial^{2}\Lambda(f(r))}{|\partial
f(r)|^{2}}|\delta f(r)|^{2}$ (4.2)

For a functional $\Lambda=\Lambda(f(r),h(r))$, the 1st variation is

$\displaystyle\delta\Lambda(f(r))=\frac{\partial\Lambda(f(r),h(r))}{\partial
f(r)}\delta f(r)+\frac{\partial\Lambda(f(r),h(r))}{\partial h(r)}\delta h(r)$
(4.3)

and the second variation is

$\displaystyle\delta\Lambda(f(r))=\frac{\partial\Lambda(f(r),h(r))}{\partial
f(r)}\delta f(r)+\frac{\partial\Lambda(f(r),h(r))}{\partial h(r)}\delta h(r)$
$\displaystyle+\frac{1}{2}\frac{\partial^{2}\Lambda(f(r),h(r))}{|\partial
f(r)|^{2}}|\delta
f(r)|^{2}+\frac{1}{2}\frac{\partial^{2}\Lambda(f(r),h(r))}{|\partial
h(r)|^{2}}|\delta h(r)|^{2}$ (4.4)

The first theorem derives the FEHE via the 1st variation of the total energy
$\bm{[11]}$.

######  Theorem 4.2.

Let a Newtonian gaseous star of mass $M$ and radius $R$ in hydrostatic
equilibrium have total energy $\mathbf{E}=\mathbb{E}_{T}+\mathbb{E}_{G}$,
occupying a ball $\mathbf{B}_{R}(0)\subset\mathbf{R}^{3}$. Let $\mathscr{S}$
be the entropy per nucleon or mass element and $(\mathcal{M}(r),p(r),\rho(r))$
have the usual definitions. Then the vanishing of the 1st variation of the
total energy is

$\displaystyle\delta\mathbf{E}=\frac{\partial\mathbf{E}}{\partial\rho(r)}\delta\rho(r)+\frac{\partial\mathbf{E}}{\partial
r}\delta r=0$ (4.5)

yields the fundamental hydrostatic equilibrium equation. The equilibrium state
is therefore a minimizer or extremum of $\mathscr{H}$.

###### Proof.

The gravitational and thermal energies of the star are given by

$\displaystyle\mathbf{E}_{G}=-4\pi\mathscr{G}\int_{0}^{R}\frac{\mathcal{M}(r)\rho(r)}{r}=-\int_{0}^{M}\frac{\mathcal{M}(r)d{\mathcal{M}(r)}}{r}$
(4.6) $\displaystyle\mathbf{E}_{T}=\int_{0}^{R}\mathcal{E}(r,\rho(r))4\pi
r^{2}dr=\int_{0}^{M}\mathscr{J}(r,\rho(r))d\mathcal{M}(r)$ (4.7)

where $\mathcal{E}$ is the internal thermal energy density and $\mathscr{J}$
is the energy per unit mass. The variation of the total energy is then

$\displaystyle\delta\mathbf{E}=\frac{\partial\mathbf{E}}{\partial\rho(r)}\delta\rho(r)+\frac{\partial\mathbf{E}}{\partial
r}\delta
r=\frac{\partial\mathbf{E}_{T}}{\partial\rho(r)}\delta\rho(r)+\frac{\partial\mathbf{E}_{G}}{\partial
r}\ \delta r=0$ (4.8)
$\displaystyle\delta\mathbf{E}=\frac{\partial}{\partial\rho(r)}\int_{0}^{M}\frac{\mathcal{E}(r,\rho(r))}{\rho(r)}d\mathcal{M}(r)\delta\rho(r)-\frac{\partial}{\partial
r}\int_{0}^{M}\frac{\mathcal{M}(r)d\mathcal{M}(r)}{r}\delta r$
$\displaystyle=\int_{0}^{M}\frac{\partial}{\partial\rho(r)}{\mathscr{J}(r,\rho(r))}d\mathcal{M}(r)\delta\rho(r)-\int_{0}^{M}\frac{\partial}{\partial
r}\left|\frac{\mathcal{M}(r)d\mathcal{M}(r)}{r}\right|\delta r$ (4.9)

Using the thermodynamic relationship

$\displaystyle
p(r)=\frac{\partial\mathcal{E}}{\partial(1/\rho(r))}=|\rho(r)|^{2}\frac{\partial\mathscr{J}}{\partial\rho(r)}$
(4.10)

for the pressure at constant entropy as well as
$\delta(\tfrac{1}{\rho(r)})=\delta(4\pi r^{2}\tfrac{dr}{d\mathcal{M}(r)})$ and
the variational calculus result
$\delta(\tfrac{dr}{df(r)})g(r)=\tfrac{dg(r)}{df(r)}\delta r$, for functions
$f(r)$ and $g(r)$, the 1st variation $\delta\mathbf{E}$ becomes

$\displaystyle\delta\mathbf{E}$
$\displaystyle=\int_{0}^{M}\frac{p(r)}{|\rho(r)|^{2}}\delta\rho(r)d\mathcal{M}(r)-{\mathscr{G}}\int_{0}^{M}\frac{\partial}{\partial
r}\frac{\mathcal{M}(r)d\mathcal{M}(r)}{r}\delta r$
$\displaystyle=\int_{0}^{M}\frac{p(r)}{|\rho(r)|^{2}}\delta\rho(r)d\mathcal{M}(r)+{\mathscr{G}}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$
$\displaystyle=-\int_{0}^{M}p(r)\delta(\frac{1}{\rho(r)})d\mathcal{M}(r)+{\mathscr{G}}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$ $\displaystyle=-\int_{0}^{M}p(r)\delta\left|4\pi
r^{2}\frac{dr}{d\mathcal{M}(r)}\right|d\mathcal{M}(r)+{\mathscr{G}}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$ $\displaystyle=-\int_{0}^{M}p(r)\delta(4\pi
r^{2})\frac{dr}{d\mathcal{M}(r)}d\mathcal{M}(r)-\int_{0}^{M}p(r)4\pi
r^{2}\delta(\frac{dr}{d\mathcal{M}(r)})d\mathcal{M}(r)+{\mathscr{G}}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$ $\displaystyle=-\int_{0}^{M}8\pi
P(r)r\frac{dr}{d{\mathcal{M}}(r)}d\mathcal{M}(r)\delta
r-\int_{0}^{M}4\pi\delta(\frac{dr}{d\mathcal{M}(r)})p(r)r^{2}d\mathcal{M}(r)+{\mathscr{G}}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$ $\displaystyle=-\int_{0}^{M}8\pi
p(r)r\frac{dr}{d\mathcal{M}(r)}\mathcal{M}(r)\delta
r-\int_{0}^{M}4\pi\frac{d}{d\mathcal{M}(r)}(p(r)r^{2})\delta
rd\mathcal{M}(r)+\mathscr{G}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$ $\displaystyle=-\int_{0}^{M}8\pi
p(r)r\frac{dr}{d\mathcal{M}(r)}d\mathcal{M}(r)\delta r-\int_{0}^{M}8\pi
p(r)r\frac{dr}{d\mathcal{M}(r)}d\mathcal{M}(r)\delta r$ (4.11)
$\displaystyle+\int_{0}^{M}4\pi r^{2}\frac{dp(r)}{d\mathcal{M}(r)}(p(r)\delta
rd\mathcal{M}(r)+\mathscr{G}\int_{0}^{M}\frac{\mathcal{M}(r)}{r^{2}}d\mathcal{M}(r)\delta
r$
$\displaystyle=\int_{0}^{M}\frac{1}{\rho(r)}\frac{dp(r)}{dr}+\int_{0}^{M}\frac{\bm{\mathscr{G}}\mathcal{M}(r)}{r^{2}}=\int_{0}^{M}\left[\frac{1}{\rho(r)}\frac{dp(r)}{dr}+\frac{\bm{\mathscr{G}}\mathcal{M}(r)}{r^{2}}\right]d\mathcal{M}(r)\delta
r=0$ (4.12)

Hence, the condition for hydrostatic equilibrium follows

$\displaystyle\underbrace{\frac{dp(r)}{dr}=-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}}$
(4.13)

∎

###  4.1. 2nd variation of the thermodynamic Massieu function for an
isothermal self-gravitating sphere of perfect gas

The condition for hydrostastic equilibrium can be established via the
isothermal gas sphere via the 2nd variation of the thermodynamic Massieu
function for a self-gravitating perfect gas. Equilibrium is the critical point
for te Massieu function. The statistical mechanics of a self-gravitating
system of particles raises many new issues and conundrums, initiated with the
work of Antonov, and which have been discussed extensively in a number of
works. (See $\bm{[7,8,9]}$] and reference their in). However, one can still
(carefully) apply statistical mechanics to a Newtonian gas of particles in a
finite domain. In what follows the canonical ensemble is utilised, whereby the
gas is at fixed temperature $\Theta$ and energy can fluctuate.

Let $\mathbf{B}(0,R)\subset\mathbf{R}^{3}$ be a spherical domain of radius of
radius $R$ and $r\leq R$ and volume
$\mathrm{V}(\mathbf{B}(0,R))=\tfrac{4}{3}\pi r^{3}$. Let the domain support an
isothermal mono-atomic gas of N particles of mass m so that the total mass of
the gas is $M=Nm$ and the temperature $\Theta$ is fixed. The system is in non-
rotating and is not expanding/collapsing. This can describe a star of N
hydrogen atoms or a globular cluster of N stars with $N\gg 0$. For a self-
gravitating gas there is now a long-range interaction described by the
Newtonian potential $\Phi(r)$. Let $\mathfrak{F}(r,\bm{u},t)$ be a
distribution function for the velocities and positions. Let $\mathbf{C}$ be a
’cell’ of the phase space $\\{x_{1},x_{2},x_{3},u_{1},u_{2},u_{3}\\}$ so that

$\displaystyle\mathbf{C}=\\{((x_{1},x_{1}+dx_{1}),(u_{2}+du_{2}),(x_{3},u_{3}+dx_{3}),(u_{1}++du_{1}),(u_{2},u_{2}+dx_{1}),(u_{3},u_{3}+dx_{1})\\}$
(4.14)

where $r=\|\bm{r}\|=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$ and
$\|\bm{u}\|=\sqrt{u_{1}^{2}+u_{2}^{2}+u_{3}^{2}}$. The volume elements is
$d\mathrm{V}=dx_{1}dx_{2}dx_{3}=4\pi r^{2}dr$ and
$d\mathcal{U}=du_{1}du_{2}du_{3}$. The distribution function
$\mathfrak{F}(\bm{r},\bm{u},t)d\mathrm{V}d\mathcal{U}$ then gives the number
of particles in the cell at time with positions between $x_{1}$ and
$x_{1}+dx_{1}$ etc. and velocities $u_{1}+du_{1}$ etc. The spatial density and
total mass are then

$\displaystyle\rho(r,t)=\int_{\mathbf{B}(0,R)}\mathfrak{F}(\bm{r},\bm{u},t)d\mathrm{V}$
(4.15) $\displaystyle M=mN=\int_{\mathbf{B}(0,R)}\rho(r)d\mathrm{V}=const.$
(4.16)

The total energy $\mathbf{E}$ is now

$\displaystyle\mathbf{E}=\frac{1}{2}\int_{\mathbf{B}(0,R)}\mathfrak{F}|\bm{u}|^{2}d\mathrm{V}+\frac{1}{2}\int_{\mathbf{B}(0,R)}\rho(r)\Phi(r)d\mathrm{V}=\mathbf{E}_{K}+\mathbf{E}_{G}=\mathbf{E}_{K}+\mathbf{E}_{G}$
(4.17)

where $\Phi(r)$ is the usual Newtonian potential. $\mathbf{E}_{K}$ is the
kinetic energy of the particles and this can be related to thermal energy for
a perfect gas to reproduce since $t\frac{1}{2}m|\bm{u}|^{2}\sim k_{b}\Theta$.
The Newtonian potential in the gas is

$\displaystyle\Phi(r)=-4\pi\mathscr{G}\int_{\mathbf{B}(0,R)}\frac{\rho(r,t)d\mathrm{V}(r)}{|\bm{r}-\bm{r}^{\prime}|}$
(4.18)

The Boltzmann entropy is the standard formula

$\displaystyle\mathcal{S}=-k_{B}\int\\!\\!\\!\int\left|\frac{\mathfrak{F}}{m}\right|\log\left|\frac{\mathfrak{F}}{m}\right|d\mathcal{U}d\mathrm{V}$
(4.19)

In a canonical ensemble in which the temperature $\Theta$ is fixed, the energy
is allowed to fluctuate. The relevant thermodynamic potential is then the
Massiue function $\mathcal{J}$ related to the Helmholtz free energy
$\mathcal{F}$. They are defined as

$\displaystyle\mathcal{F}=\mathbf{E}-\Theta\mathcal{S}$ (4.20)
$\displaystyle\mathcal{J}=-\frac{\mathcal{F}}{\Theta}=\mathcal{S}-\frac{1}{\Theta}\mathbf{E}$
(4.21)

The energy density for a perfect isothermal gas is

$\displaystyle\mathcal{E}=\frac{\frac{3}{2}Nk_{B}\Theta}{|V(\mathbf{B}(0,R))|}$
(4.22)

For a perfect isothermal gas, the total energy $\mathbf{E}$ and entropy
$\mathcal{S}$ become

$\displaystyle\mathbf{E}$
$\displaystyle=\int_{\mathrm{V}(\mathbf{B}(0,R))}\frac{\frac{3}{2}Nk_{B}\Theta}{\mathrm{V}(\mathbf{B}(0,R))}d\mathrm{V}+\frac{1}{2}\int\rho(r)\phi(r)d\mathrm{V}$
$\displaystyle=\frac{\frac{3}{2}Nk_{B}\Theta}{\mathrm{V}(\mathbf{B}(0,R))}\int_{\mathrm{V}(\mathbf{B}(0,R))}d\mathrm{V}+\frac{1}{2}\int\rho(r)\phi(r)d\mathrm{V}$
(4.23)
$\displaystyle=\frac{\frac{3}{2}Nk_{B}\Theta}{\mathrm{V}(\mathbf{B}(0,R))}\mathrm{V}(\mathbf{B}(0,R))+\frac{1}{2}\int_{\mathrm{V}(\mathbf{B}(0,R))}\rho(r)\phi(r)d\mathrm{V}$
$\displaystyle=\frac{3}{2}Nk_{B}\Theta+\frac{1}{2}\int_{\mathrm{V}(\mathbf{B}(0,R))}\rho(r)\phi(r)d\mathrm{V}\equiv\int_{\mathrm{V}(\mathbf{B}(0,R))}\mathscr{E}d\mathrm{V}+\mathbf{E}_{G}=\mathbf{E}_{T}+\mathbf{E}_{G}$
(4.24)

The entropy (4.19) for the perfect gas is

$\displaystyle\mathcal{S}=\frac{3k_{B}N}{2}+\frac{3k_{B}N}{2}\log\left(\frac{2\pi
k_{B}\Theta}{m}\right)-k_{B}\int_{\mathbf{B}(0,R)}\left|\frac{\rho(r)}{m}\right|\log\left|\frac{\rho(r)}{m}\right|d\mathrm{V}$
(4.25)

The main theorem for the isothermal gas sphere is now stated.

######  Theorem 4.3.

Let a spherical domain $\mathbf{B}(0,R)$ of radius R support a gas obeying the
perfect gas law at fixed temperature so that $d\Theta(r)/dr=0$. Then the
entropy and energy are given by (4.23) and (4.24). The Massieu thermodynamic
potential $\mathcal{J}=\mathcal{S}-\frac{\mathcal{E}}{\Theta}$ is then

$\displaystyle\mathcal{J}$
$\displaystyle=\frac{3k_{B}N}{2}+\frac{3k_{B}N}{2}\log\left(\frac{2\pi
k_{B}\Theta}{m}\right)-k_{B}\int_{\mathbf{B}(0,R)}\left|\frac{\rho(r)}{m}\right|\log\left|\frac{\rho(r)}{m}\right|d\mathrm{V}$
$\displaystyle-\frac{3}{2}Nk_{B}-\frac{1}{2\Theta}\int_{(\mathbf{B}(0,R))}\rho(r)\phi(r)d\mathrm{V}$
(4.26)

Let the 2nd variation of $\mathcal{J}$ with respect to density $\rho(r)$ be
$\bm{\delta}^{2}\mathcal{J}$. Let the first variation of the total mass with
respect to density be $\delta M=0$ so that the total mass M is an invariant.
Then the constrained optimisation condition with Lagrange multiplier $\chi$

$\displaystyle\delta^{2}\mathcal{J}-\chi\delta
M=\frac{1}{2}\frac{\partial^{2}\mathcal{J}}{\partial\rho(r)^{2}}|\delta\rho(r)|^{2}+\frac{\partial\mathcal{J}}{\partial\rho(r)}\delta\rho(r)-\chi\frac{\partial
M}{\partial\rho(r)}(\delta\rho(r))$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}\mathcal{S}}{\partial\rho(r)^{2}}|\delta\rho(r)|^{2}+\frac{\partial\mathcal{S}}{\partial\rho(r)}\delta\rho(r)$
$\displaystyle+\frac{1}{2}\frac{1}{\Theta}\frac{\partial^{2}\mathbf{E}}{\partial\rho(r)^{2}}|\delta\rho(r)|^{2}+\frac{\partial\mathbf{E}}{\partial\rho(r)}\delta\rho(r)-\chi\frac{\partial
M}{\partial\rho(r)}(\delta\rho(r))=0$ (4.27)

yields a Boltzmann distribution for the density as a critical point of the
free energy so that

$\displaystyle\rho(r)=\rho_{c}\exp\left(-\frac{mk_{B}\Phi(r)}{k_{N}\Theta}\right)=\rho_{c}\exp\left(-\frac{\Phi(r)}{\mathscr{R}}\right)$
(4.28)

The equations of hydrostatic equilibrium (2.82) and (2.83) for an isothermal
perfect-gas sphere then follow.

###### Proof.

Let $\mathbf{D}_{\rho}=\frac{\partial}{\partial\rho(r)}$ and
$\mathbf{D}^{2}_{\rho}=\mathbf{D}_{\rho}\mathbf{D}_{\rho}=\frac{\partial^{2}}{|\partial\rho(r)|^{2}}$.
The constrained optimizational problem to solve is then

$\displaystyle\delta^{(2)}\mathcal{J}-\chi\delta M$
$\displaystyle=-k_{B}\mathbf{D}_{\rho}\int_{\mathbf{B}(0,R)}\left|\frac{\rho(r)}{m}\right|\log\left|\frac{\rho(r)}{m}\right|d\mathrm{V}\delta\rho(r)$
$\displaystyle-\frac{1}{2}k_{B}\mathbf{D}_{\rho}\mathbf{D}_{\rho}\int_{\mathbf{B}(0,R)}\left|\frac{\rho(r)}{m}\right|\log\left|\frac{\rho(r)}{m}\right|d\mathrm{V}|\delta\rho(r)|^{2}$
$\displaystyle-\frac{1}{\Theta}\int_{\mathbf{B}(0,R)}\rho(r)\Phi(r)d\mathrm{V}\delta\rho(r)-\frac{1}{2\Theta}\int_{\mathbf{B}0,R)}\rho(r)\Phi(r)d\mathrm{V}\delta\rho(r)$
$\displaystyle-\frac{1}{4\Theta}\mathbf{D}^{2}_{\rho}\int_{\mathbf{B}0,R)}\rho(r0\Phi(r)d\mathrm{V}|\delta\rho(r)|^{2}-\chi\mathbf{D}_{\rho}\int_{\mathbf{B}0,R)}\rho(r)d\mathrm{V}\delta\rho(r)$
$\displaystyle=-k_{B}\int_{\mathbf{B}(0,R)}\left(\frac{\rho(r)}{m}\mathbf{D}_{\rho}(\log\rho(r)-\log
m)+\frac{1}{m}\log\left|\frac{\rho(r)}{m}\right|\right)d\mathrm{V}\delta\rho(r)$
$\displaystyle-
k_{B}\int_{\mathbf{B}(0,R)}\mathbf{D}_{\rho}\left(\frac{\rho(r)}{m}\mathbf{D}_{\rho}(\log\rho(r)-\log
m)+\frac{1}{m}\log\left|\frac{\rho(r)}{m}\right|\right)d\mathrm{V}\delta\rho(r)$
$\displaystyle-\frac{1}{\Theta}\int_{\mathbf{B}(0,R)}\Phi(r)d\mathrm{V}\delta\rho(r)-\frac{1}{4\Theta}\int_{\mathbf{B}(0,R)}\delta\Phi(r)d\mathrm{V}\delta\rho(r)-\chi\int_{\mathbf{B}(0,R)}d\mathrm{V}\delta\rho(r)$
$\displaystyle=-\frac{K_{b}}{m}\int_{\mathbf{B}(0,R)}\left(1+\log\left|\frac{\rho(r)}{m}\right|\right)d\mathrm{V}\delta\rho(r)-\frac{1}{2}\frac{k_{B}}{m}\int_{\mathbf{B}(0,R)}\frac{|\delta\rho(r)|^{2}}{\rho(r)}d\mathrm{V}$
$\displaystyle-\frac{1}{\Theta}\int_{\mathbf{B}(0,R)}\phi(r)\delta\rho(r)d\mathrm{V}-\frac{1}{2\Theta}\int_{\mathbf{B}(0,R)}\delta\rho(r)\delta\phi(r)d\mathrm{V}-\chi\int_{\mathbf{B}(0,R)}d\mathrm{V}\delta\rho(r)$
$\displaystyle=-\int_{\mathbf{B}(0,R)}d\mathrm{V}\delta\rho(r)\left\\{1+\log\left|\frac{\rho(r)}{m}\right|+\frac{\Phi(r)}{\Theta}\right\\}\delta\rho(r)d\mathrm{V}$
$\displaystyle-\frac{1}{2}\frac{K_{B}}{m}\int_{\mathbf{B}(0,R)}\frac{|\delta\rho(r)|^{2}}{\rho(r)}d\mathrm{V}-\frac{1}{2\Theta}\int_{\mathbf{B}(0,R)}\delta\rho(r)\delta\phi(r)d\mathrm{V}-\chi\int_{\mathbf{B}(0,R)}d\mathrm{V}\delta\rho(r)$
(4.29)

Then to first order

$\displaystyle\delta^{(2)}\mathcal{J}-\chi\delta
M=-\int_{\mathbf{B}(0,R)}d\mathrm{V}\delta\rho(r)\left\\{1+\log\left|\frac{\rho(r)}{m}\right|+\frac{\Phi(r)}{\Theta}+\chi\right\\}\delta\rho(r)d\mathrm{V}=0$
(4.30)

This is satisfied if

$\displaystyle\chi=-\frac{\Phi(r)}{\Theta}-\frac{k_{B}}{m}\left(\log\left|\frac{\rho(r)}{m}\right|+1\right)=C$
(4.31)

Since the Lagrange multiplier must be constant then $d\chi/dr=0$ so that

$\displaystyle\frac{d\chi}{dr}=-\frac{1}{\Theta}\frac{d\Phi(r)}{dr}-\frac{k_{B}}{m}\frac{d}{dr}\log\left|\frac{\rho(r)}{m}\right|=0$
(4.32)

Multiplying through by $dr$

$\displaystyle-\frac{1}{\Theta}d\Phi(r)\frac{k_{B}}{m}-d\log\rho(r)$ (4.33)

then integrating

$\displaystyle\frac{1}{\Theta}\int_{\Phi(0)}^{\Phi(r)}d\Phi(r)=-\frac{k_{B}}{m}\int_{\log\rho(0)}^{\log\rho(r)}d\overline{\log\rho(r)}=-\frac{k_{B}}{m}\log\left|\frac{\rho(r)}{\rho(0)}\right|$
(4.34)

so that

$\displaystyle\frac{\Phi(r)}{\Theta}=-\frac{k_{B}}{m}\log\left|\frac{\rho(r)}{\rho(0)}\right|$
(4.35)

since $\Phi(0)=0$. Solving for the density profile then gives an equilibrium
or stationary Boltzmann density distribution as the critical point of the
Massieu thermodynamic potential

$\displaystyle\rho(r)=\rho(0)\exp\left(-\frac{\Phi(r)m}{k_{B}\Theta}\right)=\rho_{c}\exp(\left(-\frac{\Phi(r)}{\mathscr{R}}\right)$
(4.36)

This must then be associated with a Newtonian gravitational potential
$\Phi(r)$ that satisfies a Boltzmann-Poisson equation

$\displaystyle\Delta\Phi(r)=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\Phi(r)}{dr}\right)=4\pi\mathscr{G}\rho_{c}\exp\left(-\frac{\Phi(r)}{\mathscr{R}}\right)$
(4.37)

describing an equilibrium or stationary configuration. Comparing with (2.83)
it can be seen that this is an isothermal perfect gas sphere. Since the gas
obeys the idea gas law then the Boltzmann density of the gas obeys the
differential equation (2.83)

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d}{dr}\log\rho(r)\right)=\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{\rho(r)}\frac{d\rho(r)}{dr}\right)=\frac{4\pi\mathscr{G}\rho(r)}{\mathscr{R}}$
(4.38)

We now essentially take the reverse of Lemma (2.21). From the isothermal
perfect gas law $p(r)=\mathscr{R}\rho(r)$ it follows that

$\displaystyle\frac{d\rho(r)}{dr}=\frac{1}{\mathscr{R}}\frac{dp(r)}{dr}$
(4.39)
$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{\mathscr{R}\rho(r)}\frac{dp(r)}{dr}\right)=\frac{4\pi\mathscr{G}\rho(r)}{\mathscr{R}}$
(4.40)

Integrating

$\displaystyle\frac{r^{2}}{\rho(r)}\frac{dp(r)}{dr}=-\int_{0}^{r}4\pi\mathscr{G}\overline{r}^{2}\rho(\overline{r})d\overline{r}$
(4.41)

then

$\displaystyle\underbrace{\frac{dp(r)}{dr}=-\frac{\left|\int_{0}^{r}4\pi\mathscr{G}\overline{r}^{2}\rho(\overline{r})d\overline{r}\right|\rho(r)}{r^{2}}\equiv-\frac{\mathscr{G}\mathcal{M}(r)\rho(r)}{r^{2}}}$
(4.42)

which is again the fundamental equation of hydrostatic equilibrium ∎

##  5\. Polytropic Gaseous Stars And The Lane-Emden Equation:Derivation,
Existence And Uniqueness Of Solutions, Analytical Solutions And Stability

There are classes of stars for which the pressure is proportional to the
density raised to a power; that is,a polytropic gas law. Such polytropic
gaseous stars (PGS) have been extensively studied beginning over a century
ago. A detailed study is still the classic 1939 work by Chandrasekhar
$\bm{[1]}$ and more recently $\bm{[12]}$. When combined with the FEHE this
leads to the Lane-Emden stellar structure equation which has both analytical
and numerical solutions. While approximations to real stars, these models have
been efficacious and useful in stellar structure theory. However, it should be
noted that there is still no complete mathematically rigorous treatment of the
nonlinear stability of polytropic gas stars.

######  Definition 5.1.

Let p be a function of $\rho>0$ such that $p\geq 0$ and $dp(\rho)/d\rho>0$ for
$\rho>0$. Let $\Xi$ be a smooth function on $\mathbf{R}$ such that $\Xi(0)=0$.
The most general polytropic gas law has the form

$\displaystyle p=\mathsf{K}|\rho|^{\lambda}(1+\Xi(\rho^{\gamma-1}))$ (5.1)

where $\mathsf{K}$ is a constant and the polytropic index is $\gamma\in(1,2)$.
Usually $\Xi=0$ so that

$\displaystyle p=\mathsf{K}|\rho|^{\lambda}=\mathsf{K}|\rho|^{1+\frac{1}{n}}$
(5.2)

One can define an enthalpy variable $\mathfrak{E}(\rho)$

$\displaystyle\mathfrak{E}=\int_{0}^{\rho}\frac{dp}{\rho}$ (5.3)

so that

$\displaystyle\frac{d\mathfrak{E}(\rho)}{d\rho}=\frac{1}{\rho}\frac{dp}{d\rho}=\frac{1}{\rho(r)}\frac{dp(r)}{dr}\frac{dr}{d\rho(r)}$
(5.4)

and the function $\psi(\rho)$

$\displaystyle\psi(\rho)=\int_{0}^{\rho}\mathfrak{E}(\rho)d\rho,~{}~{}\frac{d\psi}{d\rho}=\mathfrak{E}(\rho)$
(5.5)

For a polytropic gas

$\displaystyle\mathfrak{E}(\rho)=\frac{\mathsf{K}\gamma}{\gamma-1}|\rho|^{\gamma-1}$
(5.6)

######  Lemma 5.2.

Let a domain $\mathbf{B}\subset\mathbf{R}^{n}$ with volume
$\mathrm{V}=\mathrm{V}(\mathbf{B})$ be filled with a perfect mono-atomic gas
at temperature $\Theta$; for example, ionized hydrogen. The domain is
adiabatic or ’insulated’ so no heat enters or leaves the domain under
variations $\delta\mathrm{V}(\mathbf{B})$. For a perfect (non-relativistic)
gas comprised of $N$ particles $p=Nk_{b}\Theta\mathrm{V}^{-1}$ and the
internal energy is E=$\tfrac{3}{2}Nk_{B}\Theta$ and total (fixed) mass is
$M=mN$. This gas is then a polytropic gas with $\gamma=5/3$ so that

$\displaystyle p=\mathsf{K}|\rho|^{5/3}$ (5.7)

###### Proof.

The energy of the perfect gas is $E=\tfrac{3}{2}p\mathrm{V}$ since
$pV=Nk_{B}\Theta$ and the differential form of this

$\displaystyle d\mathbf{E}=\frac{3}{2}(\mathrm{V}dp+pd\mathrm{V})$ (5.8)

Now the 1st law of thermodynamics states that the differential $dE$ due to
$d\mathrm{V}$ at fixed pressure is $dE=-pd\mathrm{V}$ so that (5.8) becomes

$\displaystyle-pd\mathrm{V}=\tfrac{3}{2}\mathrm{V}dp+\tfrac{3}{2}pd\mathrm{V}$
(5.9)

so that upon integrating

$\displaystyle\int\frac{dp}{p}=-\frac{5}{3}\int\frac{d\mathrm{V}}{\mathrm{V}}$
(5.10)

Hence $\log p=-\frac{5}{3}\ln\mathrm{V}=\ln\mathrm{V}^{-5/3}$ and so
$p\sim\mathrm{V}^{-5/3}$. Since the mass of the gas is left unchanged and
$\rho=M/\mathrm{V}$ it follows that

$\displaystyle p\sim|\rho|^{5/3}$ (5.11)

and so an adiabatic perfect gas is also a polytropic gas with
$p=\mathsf{K}|\rho|^{5/3}$. ∎

######  Remark 5.3.

The polytropic gas model can be applied to the following types of stars:

  * •

Stars with efficient convective energy transport obey a polytropic law
$p=\mathsf{K}|\rho|^{5/3}$ so that $\gamma=5/3$ and $\mathsf{K}$ depends on
boundary conditions such as $p(0),\rho(0)$.

  * •

Light white dwarf stars with $\gamma\sim 5/3$ and very heavy white dwarf stars
with $\gamma\sim 4/3$. Then $\mathsf{K}$ depends on chemical composition and
fundamental constants, and the equation of state can expressed in polytropic
form.

  * •

Supermassive stars supported predominantly by radiation pressure with
$p=\mathsf{K}|\rho|^{4/3}$, and $\mathsf{K}$ depends on the ratio of matter
and radiation pressure.

The following formulas for $\mathbf{E}_{G}$ and $\mathbf{E}_{T}$ for a
polytropic gas star are well known and applies the FEHE $\bm{[1,10,12]}$.

######  Theorem 5.4.

($\underline{Betti~{}and~{}Ritter}$. The gravitational potential energy
$\mathbf{E}_{G}$ and the thermal energy $\mathbf{E}_{T}$ for a polytropic
gaseous star with polytropic equation of state
$p(r)=\mathsf{K}|\rho(r)|^{\gamma}=\mathsf{K}|\rho(r)|^{\frac{n+1}{n}}$ is

$\displaystyle\mathbf{E}_{G}=-\int_{0}^{M}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r}=\frac{3}{5-n}\frac{\mathscr{G}M^{2}}{R}\equiv-\frac{3(\gamma-1)}{(5\gamma-6)}\frac{\mathscr{G}M^{2}}{R}$
(5.12) $\displaystyle\mathbf{E}_{T}=\int_{0}^{R}\mathcal{E}4\pi
r^{2}dr=\frac{1}{(5\gamma-1)}\frac{\mathscr{G}M^{2}}{R}=-\frac{1}{(3\gamma-1)}\mathbf{E}_{G}$
(5.13)

The total energy is then

$\displaystyle\mathbf{E}=\mathbf{E}_{T}+\mathbf{E}_{G}=-\frac{(3\gamma-4)}{(5\gamma-6)}\frac{\mathscr{G}M^{2}}{R}$
(5.14)

###### Proof.

The proof utilises the FEHE and the polytropic gas equation. From the
polytropic state equation

$\displaystyle\frac{dp(r)}{\rho(r)}=\mathsf{K}\left(\frac{n+1}{n}\right)|\rho(r)|^{\frac{1}{n-1}}d\rho(r)=(n+1)d\left(\frac{p(r)}{\rho(r)}\right)$
(5.15)

Since there are several changes of variable, it is convenient to label the
integrals throughout as $\int_{C}^{S}$, where C refers to any quantity at the
centre and S refers to the surface. We take $p(R)=0,\rho(R)=0$ and
$\mathcal{M}(0)=0,\mathcal{M}(R)=M$.

$\displaystyle\mathbb{E}$
$\displaystyle=-\int_{0}^{M}\mathscr{G}\frac{\mathcal{M}(r)d\mathcal{M}(r)}{r}=-\underbrace{\frac{1}{2}\int_{C}^{S}\frac{\mathscr{G}\mathcal{M}(r)d|\mathcal{M}(r)|^{2}}{r}}_{int~{}by~{}parts}$
$\displaystyle=\left|\frac{\mathscr{G}{m}(r)}{2r}\right|_{C}^{S}-\frac{1}{2}\int_{C}^{S}\frac{\mathscr{G}|\mathcal{M}(r)|^{2}}{r^{2}}dr$
$\displaystyle=\left|\frac{\mathscr{G}\mathcal{M}(r)}{2r}\right|_{C}^{S}-\frac{1}{2}\int_{C}^{S}\mathcal{M}(r)\underbrace{\left(\frac{\mathscr{G}\mathcal{M}(r)}{r^{2}}dr\right)}_{use~{}FEHE}$
$\displaystyle=\frac{\mathscr{G}M^{2}}{2R}+\frac{1}{2}\int_{C}^{S}\mathcal{M}(r)\left(\frac{dp(r)}{\rho(r)}\right)$
$\displaystyle=\frac{\mathscr{G}M^{2}}{2R}+\underbrace{\frac{1}{2}(n+1)\int_{C}^{S}\mathcal{M}(r)\left(\frac{dp(r)}{\rho(r)}\right)}_{int~{}by~{}parts}$
$\displaystyle=-\frac{\mathscr{G}M^{2}}{2R}+\left|\frac{1}{2}(n+1)\mathcal{M}(r)\frac{p(r)}{\rho(r)}\right|_{C}^{S}-\frac{1}{2}(n+1)\int_{C}^{S}\frac{p(r)}{\rho(r)}d\mathcal{M}(r)$
$\displaystyle=-\frac{\mathscr{G}M^{2}}{2R}-\frac{1}{2}(n+1)\int_{C}^{S}\frac{[p(r)}{\rho(r)}d\mathcal{M}(r)$
$\displaystyle-\frac{\mathscr{G}M^{2}}{2R}-\frac{1}{2}(n+1)\int_{C}^{S}\frac{p(r)}{\rho(r)}4\pi
r^{2}\rho(r)dr$
$\displaystyle=-\frac{\mathscr{G}M^{2}}{2R}-\underbrace{\frac{1}{2}(n+1)\int_{C}^{S}p(r)\frac{4\pi}{3}d|r^{3}|}_{int~{}by~{}parts}$
$\displaystyle-\frac{\mathscr{G}M^{2}}{2R}-\left|\frac{1}{2}(n+1)\frac{4\pi}{3}p(r)r^{3}\right|_{C}^{S}+\frac{1}{6}(n+1)\int_{C}^{S}\underbrace{4\pi
r^{2}dp(r)}_{use~{}FEHE}$
$\displaystyle=-\frac{\mathscr{G}M^{2}}{2R}-\frac{1}{6}(n+1)\int_{C}^{S}4\pi
r^{3}\frac{\mathscr{G}\mathcal{M}(r)}{r^{2}}\rho(r)dr$
$\displaystyle=-\frac{\mathscr{G}M^{2}}{2R}-\frac{1}{6}(n+1)\int_{C}^{S}\frac{\mathscr{G}\mathcal{M}(r)d\mathcal{M}(r)}{r}$
$\displaystyle=-\frac{\mathscr{G}M^{2}}{2R}-\frac{1}{6}(n+1)\mathbf{E}_{G}$
(5.16)

so

$\displaystyle\mathbf{E}_{G}=-\frac{GM^{2}}{2R}+\frac{1}{6}(n+1)\mathbf{E}_{G}$
(5.17)

which then gives (5.12) as required. To prove (5.12), first note that for a
polytropic gas the energy density $\mathcal{E}$ is related to the pressure by
$\mathcal{E}(r)=(\gamma-1)^{-1}p(r)$. Then using the FEHE

$\displaystyle\mathbf{E}_{G}=4\pi\int_{0}^{R}r^{3}\frac{dp(r)}{dr}dr=-12\pi\int_{0}^{R}r^{2}p(r)dr=-12\pi\int_{0}^{R}\mathcal{E}(\gamma-1)r^{2}dr$
(5.18)

Hence

$\displaystyle\mathbf{E}_{G}=-(\gamma-1)\mathbf{E}_{T}$ (5.19)

and so (5.13) immediately follows using (5.12). The total energy $\mathbf{E}$
also immediately follows.

######  Corollary 5.5.

It can be see that $\mathbf{E}=0$ if $\gamma=\tfrac{4}{3}$ so such polytropic
stars are teetering on the brink of catastrophic instability. Stability at
least requires $\mathbf{E}<0$.

∎

###  5.1. The Lane-Emden equations

Since the polytropic gas law is temperature independent, then the structure of
such a star can be described by the FEHE in the Poisson form

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{\rho(r)}\frac{dp(r)}{dr}\right)=-4\pi\mathscr{G}\rho$
(5.20)

Using $p(r)=\mathsf{K}|\rho(r)|^{\gamma}$ this becomes

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{\rho(r)}\frac{d}{dr}|\rho(r)|^{\gamma}\right)+\frac{4\pi\mathscr{G}}{\mathsf{K}}r^{2}\rho=0$
(5.21)

with the boundary conditions that $\rho(0)$ has some prescribed value and
$\rho^{\prime}(0)=0,\rho(R)=0$. We can state the following theorem for a
polytropic gas equation of state $\bm{[13]}$.

######  Lemma 5.6.

Let $\rho(r)$ be a solution for $r_{o}\leq r<R$ and let $[r_{o},R]$ be a right
maximal interval of existence for $\rho>0,R<\infty,|dp(r)/dr|_{r=r_{o}}<0$.
Then $\exists C>0$ such that

$\displaystyle\rho(r)=C|R-r|^{\frac{1}{\gamma-1}}\left(1+\left\llbracket\left|\frac{R-r}{r}\right|,\overline{C}\left(\frac{R-r}{R}\right)^{\frac{\gamma}{\gamma-1}}\right\rrbracket\right)$
(5.22)

where $\overline{C}=QR^{\frac{\gamma}{\gamma-1}}C^{2-\gamma}$ and
$Q=4\pi\mathscr{G}(\gamma-1)/\mathsf{K}\gamma$.
Here,$\big{\llbracket}X\rrbracket$ is a notation denoting a power series of
the form $\llbracket X\rrbracket_{q}=\sum_{j\geq q}a_{j}X^{j}$ and so

$\displaystyle\big{\llbracket}X,Y\big{\rrbracket}_{q}=\sum_{j+k\geq
q}q_{jk}X^{j}X^{k}$ (5.23)

####  5.1.1. Non-dimensional reduction

The Emden equation can be reduced to a dimensionless form, the Lane-Emden
equation. This is a very standard derivation in the literature $\bm{1,2]}$.
Defining a variable $\theta$ by $\theta=|\rho(r)/\rho(0)|^{\gamma-1}$ then

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\theta(r)}{dr}\right)+\frac{4\pi\mathscr{G}(\gamma-1)}{\mathsf{K}\gamma}|\rho(0)|^{2-\gamma}|\theta(r)|^{1/(\gamma-1)}$
(5.24)

Defining a second variable $\xi(r)$ by

$\displaystyle\xi(r)=\left|\frac{4\pi\mathscr{G}(\gamma-1)}{\mathsf{K}\gamma}\right|^{1/2}|\rho(0)|^{(2-\gamma)/2}r$
(5.25)

the differential equation becomes

$\displaystyle\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d}{d\xi}\theta(\xi)\right)+|\theta(\xi)|^{1/(\gamma-1)}=0$
(5.26)

with boundary conditions $\theta(0)=1,\theta^{\prime}(0)=0$ and
$\theta(\xi_{1})=0$ on the surface. Since $n=\frac{1}{\gamma}-1$ the LE
equation is identically.

$\displaystyle\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d}{d\xi}\theta(\xi)\right)=-|\theta(\xi)|^{n}$
(5.27)

This is then the Lane-Emden equation, which can be solved analytically or
numerically depending on the polytropic index $n$. The following lemma
establishes that if $\theta(\xi)=0$ then one must have
$[d\theta(\xi)/d\theta]_{\xi=0}=\theta^{\prime}(0)=0$; that is, the derivative
must vanish at $\xi=0$ for bounded solutions that are nonsingular at $\xi=0$.
The function $\theta(\xi)$ is analytic at $\xi=0$.

######  Lemma 5.7.

Given the Lane-Emden equation

$\displaystyle\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta(\xi)}{d\theta}\right)\equiv\frac{1}{\xi^{2}}(\xi^{2}\theta^{\prime}(\xi))^{\prime}=f(\theta(\xi))$
(5.28)

where $f(\theta(\xi))=-|\theta(\xi)|^{n}$. Let $f\in C(\mathbf{R})$ and
$\theta\in C^{2}(0,Q)$ for some $Q>0$ and $\theta\in[0,Q]$ and
$\frac{1}{\xi^{2}}(\xi^{2}\theta^{\prime}(\xi))^{\prime}=f(\theta)$. Then
$\lim_{\xi\uparrow\infty}\theta^{\prime}(\xi)=0$ if $\theta(0)=0$.

###### Proof.

Let $\xi\in[0,Q]$ and express the LE equation as

$\displaystyle(\xi^{2}\theta^{\prime}(\xi))^{\prime}=\xi^{2}f(\theta(\xi))$
(5.29)

then integrating

$\displaystyle\xi^{2}\theta^{\prime}(\xi)=\int_{\xi_{o}}^{\xi}\beta^{2}f(\theta(\beta))d\beta+\xi_{o}^{2}\theta(\xi_{o})$
(5.30)

and

$\displaystyle
0=\int_{\xi_{o}}^{0}\beta^{2}f(\theta(\beta))d\beta+\xi_{o}^{2}\theta(\xi_{o})$
(5.31)

so that

$\displaystyle\int_{\xi_{o}}^{0}\beta^{2}f(\theta(\beta))d\beta=-\xi_{o}^{2}\theta(\xi_{o})$
(5.32)

Then (5.29) becomes

$\displaystyle\xi^{2}\theta^{\prime}(\xi)$
$\displaystyle=\int_{\xi_{o}}^{\xi}\beta^{2}f(\theta(\beta))d\beta+\xi_{o}^{2}\theta(\xi_{o})$
$\displaystyle=\int_{\xi_{o}}^{0}\beta^{2}f(\theta(\beta)d\beta+\int_{0}^{\xi}\beta^{2}f(\theta(\beta))d\beta+\xi_{o}^{2}\theta(\xi_{o})$
$\displaystyle=-\int_{\xi}^{0}\beta^{2}f(\theta(\eta)d\eta+\int_{0}^{\xi}\beta^{2}f(\theta(\beta))d\beta+\xi_{o}^{2}\theta(\xi_{o})$
$\displaystyle=-\xi_{o}^{2}\theta(\xi_{o})+\xi_{o}^{2}\theta(\xi)+\int_{0}^{\xi}\beta^{2}f(\theta(\beta))d\beta=\int_{0}^{\xi}\beta^{2}f(\theta(\beta))d\beta$
(5.33)

Now let $\mathcal{L}$ be the limit as$\xi\rightarrow\infty$

$\displaystyle\lim_{\xi\uparrow 0}\xi^{2}\theta^{\prime}(\xi)=\mathcal{L}$
(5.34)

and consider two cases:

  1. (1)

Assume that the limit is strictly greater than zero so that $\mathcal{L}>0$.
Then

$\displaystyle\lim_{\xi\uparrow
o}\theta^{\prime}(\xi)=\frac{\mathcal{L}}{\xi^{2}}$ (5.35)

If some $\alpha\geq 0$ then

$\displaystyle\lim_{\xi\uparrow
0}\theta^{\prime}(\xi)\geq\frac{\mathcal{L}}{\alpha\xi^{2}}$ (5.36)

and let $\xi\in[0,\epsilon]$. Integrating

$\displaystyle(\theta(\epsilon)-\theta(0))=\lim_{t\uparrow\infty}(\theta(\epsilon)-\theta(\xi))\geq\lim_{t\uparrow\infty}\int_{\xi}^{\epsilon}\frac{\mathcal{L}}{\alpha\beta^{2}}d\beta=\lim_{t\uparrow\infty}\frac{\mathcal{L}}{\alpha\xi}-\frac{\mathcal{L}}{\alpha\epsilon}=-\infty$
(5.37)

This then contradicts the requirement that $\theta(0)=0$. The derivative at
$\xi=0$ would then be singular and so does not exist.

  2. (2)

The case $\mathcal{L}=0$. Now $\lim_{\xi\uparrow
0}\theta^{\prime}(\xi)=\mathcal{L}=0$. Using (5.30) it is then required that

$\displaystyle\lim_{\xi\uparrow
0}\xi^{2}\theta^{\prime}(\xi)=\lim_{\xi\uparrow
0}\int_{0}^{\xi}\beta^{2}f(\theta(\beta))d\eta=0$ (5.38)

or

$\displaystyle\lim_{\xi\uparrow
0}\theta^{\prime}(\xi)=\frac{1}{\xi^{2}}\lim_{\xi\uparrow
0}\int_{0}^{\xi}\beta^{2}f(\theta(\beta))d\beta=0$ (5.39)

In terms of the supremum of $f(\theta)$

$\displaystyle\lim_{\xi\uparrow
0}\theta^{\prime}(\xi)\leq\frac{1}{\xi^{2}}\lim_{\xi\uparrow
0}\int_{0}^{\xi}\beta^{2}\sup|f(\theta(\beta))|d\beta=\frac{1}{\xi^{2}}\lim_{\xi\uparrow
0}\int_{0}^{\xi}\beta^{2}C|d\beta=\frac{1}{\xi^{2}}\int_{0}^{\xi}\beta^{2}C=\lim_{\xi\uparrow
0}\frac{1}{3}C\xi=0$ (5.40)

so that $\lim_{\xi\uparrow 0}\theta^{\prime}(\xi)=0$ as required.

∎

Within the LE formalism, the total energy can also be expressed as

$\displaystyle\mathbb{E}=\frac{4-3\gamma}{\gamma-1}\int_{0}^{R}p(r)4\pi
r^{2}dr$ (5.41)
$\displaystyle=\frac{4-3\gamma}{\gamma-1}4\pi\mathsf{K}\left(\frac{\mathsf{K}\gamma}{4\pi\mathscr{G}(\gamma-1)}\right)^{3/2}|\rho(0)|^{\frac{3\gamma-6}{2}}\int_{0}^{\xi_{1}}|\theta(\xi)|^{n+1}\xi^{2}d\xi$
(5.42)

For $\gamma>6/5$, the solution vanishes at $\xi=\xi_{1}$ so the radius of a
polytropic star is

$\displaystyle
R=\left(\frac{4\pi\mathscr{G}(\gamma-1)}{\mathsf{K}\gamma)^{-1/2}}\right)|\rho(0)|^{(2-\gamma)/2}\xi$
(5.43)

The mass is

$\displaystyle M=\int_{0}^{R}4\pi r^{2}\rho(r)dr$
$\displaystyle=|4\pi\rho(0)|^{(3\gamma-4)/2}\left(\frac{\mathsf{K}\gamma}{4\pi\mathscr{G}(\gamma-1)}\right)^{3/2}\int_{0}^{\xi_{1}}\xi^{2}|\theta(\xi)|^{1/\gamma-1}d\xi$
$\displaystyle=|4\pi\rho(0)|^{(3\gamma-4)/2}\left(\frac{\mathsf{K}\gamma}{4\pi\mathscr{G}(\gamma-1)}\right)^{3/2}\xi_{1}^{2}\left|\frac{d\theta(\xi)}{d\xi}\right|_{\xi=\xi_{1}}$
(5.44)

These formulas are useful in determining the masses and radii of white dwarf
stars.

###  5.2. Existence and uniqueness

Having formulated the Lane-Emden equations, it is necessary to establish that
a solution exists and is unique. The following general theorem appears in
$\bm{[14]}$.

######  Theorem 5.8.

Let $u(x)$ be a positive solution to the general Lane-Emden equation

$\displaystyle\Delta u+u^{n}=0,~{}~{}x\in\mathbf{R}^{N},N>2$ (5.45)

If

$\displaystyle 1\leq n<\frac{N+2}{N-2}$ (5.46)

then $u(x)=0$. Otherwise there is a non-trivial positive solution.

The case N=2 is relevant to stellar structure which is our concern here, but
for $n=(N+2)/(N-2)$ and for $N=4$ it is relevant to Yang Mills theory, while
the case $N=2$ is relevant to differential geometry. The coupled Lane-Emden
system

$\displaystyle\Delta u+v^{n}=0$ $\displaystyle\Delta v+u^{m}=0$

also arises in many mathematical, biological, chemical and physical problems.
$\bm{[15-20]}$

We are really concerned with the existence and uniqueness of solutions of the
LE equation in the spherically symmetric ’astrophysical form’ so that

$\displaystyle\Delta_{\xi}\theta(\xi)=\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta(\xi)}{d\xi}\right)\equiv\theta^{\prime\prime}(\xi)+\frac{2}{\xi}\theta^{\prime}(\xi)+|\theta(\xi)|^{n}=0,n>0$
(5.47) $\displaystyle\theta(0)=\alpha,~{}~{}~{}\theta^{\prime}(0)=0$ (5.48)

This equation can be considered as a special case of a general equation of the
form

$\displaystyle\theta^{\prime\prime}(\xi)+\mathcal{P}(\xi)\theta^{\prime}(\xi)+\mathcal{Q}(\theta(\xi),\xi)=0,~{}~{}\xi>0$
(5.49) $\displaystyle\theta(0)=\alpha,~{}~{}\theta^{\prime}(0)=0$ (5.50)

where $\mathcal{P}:\mathbf{R}\rightarrow\mathbf{R}$,
$\mathcal{Q}:\mathbf{R}\times\mathbf{R}\rightarrow\mathbf{R}$, $\alpha>0$ and
$p(\xi)$ can be singular at $\xi=0$ $\bm{[21]}$.

######  Proposition 5.9.

The following general conditions are imposed on $\mathcal{P}$ and
$\mathcal{Q}$

  1. (1)

$\mathcal{P}$ is measurable on $[0,1]$,

  2. (2)

$\mathcal{P}\geq 0$ on $(0,1]$.

  3. (3)

$\int_{0}^{1}s\mathcal{P}(s)ds<\infty$

  4. (4)

$\mathcal{Q}$ satisfies the general Caratheodry conditions and is bounded so
that:

     * •

$\exists,a,b$ with $a<\alpha<\beta$ and $C>0$ such that for all
$t\in(0,1],\mathcal{Q}(\xi,\bullet)$ is continuous on $(a,b)$.

     * •

For all $x\in(a,b)$, $\mathcal{Q}(\xi,\bullet)$ is measurable on $(0,1]$.

     * •

$\sup_{(t,x)\in[0,1]\times[a,b]}|\mathcal{Q}(\xi,x)|<C.$

######  Definition 5.10.

$\theta(\xi)$ is considered a solution to the LE equation iff $\exists~{}T>0$
such that:

  1. (1)

$\theta$ and $\theta^{\prime}$ are absolutely continuous on $(0,T]$

  2. (2)

$\theta$ satisfies (5.47) almost everywhere on $(0,T]$.

  3. (3)

$\theta$ satisfies the conditions of (5.48).

######  Lemma 5.11.

Given $\mathcal{P}(\xi)$, define a function $\psi(\xi)$ by the integral

$\displaystyle\psi(\xi)=\exp\left(\int_{1}^{\xi}\mathcal{P}(s)ds\right),~{}\xi>0$
(5.51)

where $\psi(\xi)$ is a nonnegative, nondecreasing and continuous function
which is bounded then the ODE (5.47) can be expressed in the form

$\displaystyle\frac{1}{\psi(\xi)}\frac{d}{d\xi}(\psi(\xi)\theta^{\prime}(\xi))=-\mathcal{Q}(\theta(\xi),\xi)$
(5.52)

###### Proof.

Since

$\displaystyle\mathcal{P}(\xi)=\frac{\psi^{\prime}(\xi)}{\psi(\xi)}$ (5.53)

then (-) becomes

$\displaystyle\theta^{\prime\prime}(\xi)+\frac{\psi^{\prime}(\xi)}{\psi(\xi)}\theta^{\prime}(\xi)+\mathcal{Q}(\theta(\xi),\xi)=0,~{}~{}\xi>0$
(5.54) $\displaystyle\theta(0)=\alpha,~{}~{}\theta^{\prime}(0)$
$\displaystyle=\frac{1}{\psi(t)}\big{(}\psi(\xi)\theta^{\prime\prime}(\xi)+\psi^{\prime}(\xi)\theta^{\prime}(\xi)\big{)}+\mathcal{Q}(\theta(\xi),\xi)=0$
(5.55)

which is (5.52) ∎

The standard Lane-Emden equation is the case $\psi(\xi)=\xi^{-2}$.

######  Theorem 5.12.

Let the conditions of Prop 4.9 be satisfied. Then a solution to the boundary
value problem (5.47) exists.

###### Proof.

The full proof is quite detailed and is given in $\bm{[21]}$. ∎

We will however, consider the uniqueness problem in detail. First two
preliminary lemmas are required.

######  Lemma 5.13.

The solution of the ODE

$\displaystyle\frac{1}{\psi(\xi)}\frac{d}{d\xi}(\psi(\xi)\theta^{\prime}_{i}(\xi))=-\mathcal{Q}(\theta_{i}(\xi),\xi),~{}i=1,2$
(5.56)

has the integral representation

$\displaystyle\theta_{i}(\xi)=\alpha-\int_{0}^{\xi}\left(\int_{s}^{\xi}\frac{\psi(s)}{\psi(\xi)}d\xi\right)\mathcal{Q}(\theta_{i}(s),s)ds$
(5.57)

###### Proof.

From (-)

$\displaystyle\psi(\xi)\frac{d\theta(\xi)}{d\xi}=-\int_{0}^{\xi}\psi(\xi)\mathcal{Q}(\xi,\theta_{i}(\xi))ds$
(5.58)

so that

$\displaystyle\frac{d\theta(\xi)}{d\xi}=-\frac{1}{\psi(\xi)}\int_{0}^{\xi}\psi(\xi)\mathcal{Q}(\xi,\theta_{i}(\xi))ds$
(5.59)

Integrating once more for $\xi\in[0,T]$ gives

$\displaystyle\theta(\xi)=\alpha-\int_{0}^{t}\frac{1}{\psi(\xi)}\left(\int_{0}^{\xi}\psi(s)\mathcal{Q}(\xi,\theta(\xi))ds\right)d\xi$
$\displaystyle=\alpha-\int_{0}^{\xi}\left(\int_{s}^{\xi}\frac{\psi(s)}{\psi(\xi)}d\xi\right)\mathcal{Q}(\theta_{i}(s),s)ds$
(5.60)

where the integration order has been changed. ∎

The next lemma is the basic Gronwall lemma

######  Lemma 5.14.

Let $f:\mathbf{R}\rightarrow\mathbf{R}$ and let $C_{1},C_{2}>0$. Then if

$\displaystyle f(x)\leq C_{1}+C_{2}\int_{I}f(\bar{x})d\bar{x}$ (5.61)

where $I=[x_{0},x]$ is an interval then one has the estimate or bound

$\displaystyle f(x)\leq
C_{1}\exp\left(C_{2}\int_{I}dx\right)=C_{1}\exp(C_{2}|x-x_{0}|)$ (5.62)

and $f(x)\leq 0$ if $C_{1}=0$. The result holds with $>$ replacing $\leq$.

The main theorem on uniqueness of solutions can now be stated and proved.

######  Theorem 5.15.

Let the conditions of definition (5.10) hold and let $\mathcal{Q}$ be
Lipschitz on $[a,b]$. Then the nonlinear elliptic boundary value problem
(5.47) has a unique solution.

###### Proof.

If $\mathcal{Q}$ is Lipschitz on $\theta(\xi)$ on an interval $[a,b]$ then
$\exists$ constant $\Lambda>0$ such that

$\displaystyle|\mathcal{Q}(\xi,\theta_{1}(\xi)-\mathcal{Q}(\xi_{2},\theta_{2}(\xi)|\leq\Lambda|\theta_{1}(\xi)-\theta_{2}(\xi)|$
(5.63)

for all $\theta_{1},\theta_{2}\in(a,b)$ and $\xi\in[0,1]$, where
$\theta_{1},\theta_{2}$ are any two solutions to (5.47) on $[0,T]$. If
$s\leq\xi\leq\zeta$ then $\psi(s)/\psi(\xi)|\ll 1$. Now using the integral
representation (5.57)

$\displaystyle|\theta_{2}(\zeta)-\theta_{1}(\zeta)|$
$\displaystyle=\left|\int_{0}^{\zeta}\left(\int_{s}^{\zeta}\frac{\psi(s)}{\psi(\xi)}d\xi\right)[\mathcal{Q}(\theta_{2}(s),s)-\mathcal{Q}(s,\theta(s))ds\right|$
$\displaystyle=\int_{0}^{\zeta}\left(\int_{s}^{\zeta}d\xi\right)\big{|}\mathcal{Q}(\theta_{2}(s),s)-\mathcal{Q}(\theta_{1}(s),s)\big{|}ds$
$\displaystyle\leq\int_{0}^{\zeta}\left(|\zeta-s|\right)\big{|}\mathcal{Q}(\theta_{2}(s),s)-\mathcal{Q}(\theta_{1}(s),s)\big{|}ds$
$\displaystyle\leq\Lambda\int_{0}^{\zeta}\left(|\zeta-s|\right)\big{|}\theta_{2}(s)-\theta_{1}(s)\big{|}ds$
$\displaystyle\leq\Lambda\int_{0}^{\zeta}B\big{|}\theta_{2}(s)-\theta_{1}(s)\big{|}ds$
(5.64)

Applying the Gronwall Lemma to

$\displaystyle|\theta_{2}(\zeta)-\theta_{1}(\zeta)|\leq\int_{0}^{\zeta}C_{2}\big{|}\theta_{2}(s)-\theta_{1}(s)\big{|}ds$
(5.65)

where $C_{2}=\lambda B$ and $C_{1}=0$ gives
$|\theta_{2}(\zeta)-\theta_{1}(\zeta)|\leq 0$ so that one must have
$\theta_{1}(\zeta)=\theta_{2}(\zeta)$ for any two solutions. Hence, the
solution to the Lane-Emden nonlinear elliptic boundary value problem is
unique. ∎

###  5.3. Transformations of the LE equations and a dynamical systems
perspective

The LE equations can be transformed to alternative forms which are useful for
their solution $\bm{[1,12]}$. They can be recast into an autonomous 1st-order
system of ODEs much like a dynamical system.

######  Lemma 5.16.

The transformation $\theta=\chi/\xi$ reduces the LE equations to

$\displaystyle\frac{d^{2}\xi}{d\xi^{2}}=-\frac{\chi^{n}}{\xi^{n-1}}$ (5.66)

The Kelvin transform $x=1/\xi$ reduces the LE equation to

$\displaystyle x^{4}\frac{d{2}\theta}{dx^{2}}=-|\theta|^{n}$ (5.67)

with the solution

$\displaystyle\theta=\alpha
x^{\eta},~{}~{}\eta=\frac{2}{n-1},~{}~{}\alpha=\left(\frac{2(n-3)}{(n-1)}\right)^{\frac{1}{n-1}},~{}~{}n>3$
(5.68)

###### Proof.

It is easy to show that (5.66) is a solution of (5.65). ∎

######  Lemma 5.17.

The substitution

$\displaystyle\theta(x)=\beta x^{n}w(x)$ (5.69)

transforms $x^{4}\frac{d{2}\theta}{dx^{2}}=-|\theta|^{n}$ to the form

$\displaystyle x^{2}\frac{d^{2}w}{dx^{2}}+2\eta
x\frac{dw}{dx}+\eta(\eta-1)w+|w|^{n}=0$ (5.70)

###### Proof.

The derivatives of (5.67) are

$\displaystyle\frac{d\theta}{dx}=\beta x^{\eta}\frac{dw}{dx}+\eta\beta
x^{\eta-1}w(x)$ (5.71)
$\displaystyle\frac{d^{2}\theta}{dx^{2}}=\beta(x^{\eta}\frac{d^{2}w}{dx^{2}}+2\eta
x^{\eta-1}\frac{dw}{dx}+\eta(\eta-1)x^{\eta-2}w$ (5.72)

Substituting into (5.65) then gives (5.68) ∎

######  Lemma 5.18.

Let

$\displaystyle x=\frac{1}{\xi}=exp(w),~{}~{}w=\log x=-\log\xi,~{}~{}\lambda=1$
(5.73)

then (5.68) becomes

$\displaystyle\frac{d^{2}w}{dx^{2}}+\left(\frac{5-n}{n-1}\right)\frac{dw}{dz}-\frac{2(n-3)}{(n-1)^{2}}w+|w|^{n}=0$
(5.74)

or

$\displaystyle\frac{d^{2}w}{dx^{2}}=-2\Xi\left(w,\frac{dw}{dz}\right)$ (5.75)

where

$\displaystyle\Xi\left(w,\frac{dw}{dz}\right)=\frac{1}{2}\left(-\frac{n-5}{n-1}\right)\frac{dw}{dz}+\frac{2(3-n)}{(n-1)^{2}}w+|w|^{n}$
(5.76)

######  Corollary 5.19.

Equation (5.72) can be recast as a first-order system of autonomous ODES, much
like a dynamical system so that

$\displaystyle\frac{dw}{dz}=\mathcal{Q}$ (5.77)
$\displaystyle\frac{d\mathcal{Q}}{dz}=-2\Xi\left(w,\frac{dw}{dz}\right)\equiv-2\Xi\left(w,\mathcal{Q}\right)$
(5.78)

and the critical points then lie along

$\displaystyle\Xi(w,0)=0$ (5.79)

These can be used to study stability. (REF)

###  5.4. explicit solutions

The Lane-Emden system $LE_{n}$ can be solved analytically $\bm{[1,2,12]}$ for
three values of $\gamma=1+\frac{1}{n}$, namely $n=0,1,5$.

  * •

If $n=0$ or $\gamma=\infty$ the LE equation is linear and homogenous with
solution $\theta(\xi)=-\tfrac{1}{6}\xi^{2}$ plus any linear combination of
$\xi$ and $1$. Then $LE_{0}$ is

$\displaystyle\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta(\xi)}{d\xi}\right)+1=0$
(5.80)

Rearranging and integrating once

$\displaystyle\xi^{2}\frac{d\theta(\xi)}{d\xi}=c_{1}-\frac{1}{3}\xi^{3}$
(5.81)

Diving through by $\xi^{2}$ and integrating once more

$\displaystyle\theta(\xi)=c_{o}-\frac{c_{1}}{\xi}-\frac{1}{6}\xi^{2}$ (5.82)

The condition $\theta(0)=1$ fixes the solution to be
$\theta(\xi)=1-\tfrac{1}{6}\xi^{2}$ with $\theta(0)=\sqrt{6}$ and
$\xi_{1}^{2}|d\theta(\xi)/d\xi|_{\xi=\xi_{1}}=-3\sqrt{6}$. So
$c_{0}=1,c_{1}=0$.

  * •

If $\gamma=2$ then $n=1$ so that

$\displaystyle\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta(\xi)}{d\xi}\right)+\theta=0$
(5.83)

Using a power series ansatz $\theta(\xi)=\sum_{n=0}^{\infty}a_{n}\xi^{n}$
leads to recurrence relation $a_{n+2}=-a_{n}/(n+3)(n+2)$. The solution is

$\displaystyle\theta(\xi)=\frac{a}{\xi}\sin\xi+\tfrac{b}{\xi}\cos\xi$ (5.84)

The condition $\theta(0)=1$ fixes the solution to be
$\theta(\xi)=(\sin\xi)/\xi$. with $\xi_{1}=\pi$ and
$\xi_{1}^{2}|d\theta(\xi)/d\xi|_{\xi=\xi_{1}}=-\pi$

  * •

If $\gamma=6/5$ or $n=5$ the solution is
$\theta(\xi)=(1+\tfrac{1}{3}\xi^{2})^{-1/2}$ with $\theta(\xi)\rightarrow 0$
as $\xi\rightarrow\infty$ so $\xi_{1}=\infty$ and
$\xi_{1}^{2}|d\theta(\xi)/d\xi|_{\xi=\xi_{1}}\rightarrow-\sqrt{3}$ for
$t\rightarrow\infty$. From (5.41) and (5.42) such polytropic stars have
infinite radius but finite mass. For $n=5$, equation (-) reduces to

$\displaystyle\frac{d^{2}w}{dz^{2}}=\frac{1}{4}w(1-w^{4})$ (5.85)

Multiply both sides by $dw/dz$ so that

$\displaystyle\frac{1}{2}\frac{d}{dz}\left(\left(\frac{dw}{dz}\right)^{2}\right)=\frac{1}{4}w(1-w^{4})\frac{dw}{dz}$
(5.86)

Integrating

$\displaystyle\left(\frac{dw}{dz}\right)^{2}=\frac{1}{2}\int(w+w^{5})dw=\frac{1}{4}w^{2}-\frac{1}{12}w^{6}+C$
(5.87)

then (choosing the negative so that $w\rightarrow\infty$)

$\displaystyle\frac{dw}{w(1-\frac{1}{3}w^{4})^{1/2}}=-\frac{1}{2}dt$ (5.88)

Now make the substitution

$\displaystyle\frac{1}{3}w^{4}=\sin^{2}\zeta$ (5.89)

First integrate it so that

$\displaystyle\frac{4}{3}w^{3}dw=2\sin\zeta\cos\zeta d\zeta$ (5.90)

then it follows that

$\displaystyle\frac{dw}{w}=\frac{1}{2}\frac{\cos\zeta}{\sin\zeta}d\zeta$
(5.91)

from (5.88) and (5.89) we obtain

$\displaystyle cosec\zeta d\zeta=-dz$ (5.92)

Integrating

$\displaystyle\int cosec\zeta
d\zeta=\log[\sin(\tfrac{1}{2}\zeta)]-\log[\cos(\tfrac{1}{2}\zeta)]=\log[\tan(\frac{1}{2}\zeta]=-\int
dz$ (5.93)

so

$\displaystyle\tan(\frac{1}{2}\zeta)=C\exp(-z)$ (5.94)

Then

$\displaystyle\frac{1}{3}w^{4}=\frac{4\tan^{2}(\frac{1}{2}\zeta)}{1+\tan^{2}(\frac{1}{2}\zeta))^{2}}$
(5.95)

which is

$\displaystyle
w=\left(\frac{12\tan^{2}(\frac{1}{2}\zeta)}{1+\tan^{2}(\frac{1}{2}\zeta))^{2}}\right)^{1/4}$
(5.96)

Using (5.92)

$\displaystyle
w=\left(\frac{12C^{2}\exp(-2z)}{1+C^{2}\exp(-2z))^{2}}\right)^{1/4}$ (5.97)

then

$\displaystyle\theta(\xi)=\left(\frac{3C^{2}}{1+C^{2}\xi^{2})^{2}}\right)^{1/4}$
(5.98)

Since $\theta(-)=1$ it follows that $C^{2}=1/3$ so the solution is

$\displaystyle\theta(\theta)=\left(\frac{1}{1+\frac{1}{3}\xi^{2}}\right)^{1/4}$
(5.99)

  * •

Other values of n or $\gamma$, also the most interesting and astro-physically
relevant, require numerical solutions. The case $\tfrac{1}{2}<\gamma<1$ or
$3<\gamma<2$ can model neutron stars. The case $\gamma=5/3$ or $n=3/2$
describes stars with convective cores such as red giants and brown dwarfs and
stars with perfect gas compositions.

  * •

The case $\gamma=1$ or $n=\infty$ corresponds to an isothermal gas sphere and
will be considered separately.

###  5.5. The isothermal gas sphere again

The isothermal gas sphere was previously discussed. Here, it is reprised again
for the case of a polytropic gas and can be described by a Lane-Emden type
equation known as the Chandrasekhar equation $\bm{[1]}$.

######  Lemma 5.20.

Let a spherical region $\mathbf{B}(0,R)$ of radius R, centre zero, support a
polytropic gas with pressure $p(r)$ and density $\rho(r)$ and polytropic index
$\gamma=1$ or $n=\infty$. Then
$p=\mathsf{K}|\rho|^{\gamma}=\mathsf{K}|\rho|^{1+\frac{1}{n}}$. The
(Poisson)equation of hydrostatic equilibrium of the isothermal sphere is then
given by (2.80) and (2.81) as

$\displaystyle\Delta\log(\rho)=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d}{dr}(\log\rho)\right)=-\frac{4\pi\mathscr{G}\rho(r)}{\mathsf{K}}$
(5.100)

or equivalently

$\displaystyle\Delta\Phi(r)=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\Phi(r)}{dr}\right)=4\pi\mathscr{G}\rho_{c}\exp\left(-\frac{\Phi(r)}{\mathsf{K}}\right)$
(5.101)

and there is a singular solution

$\displaystyle\rho(r)=\frac{A}{2\pi\mathscr{G}r^{2}}$ (5.102)

Using the substitutions

$\displaystyle\rho=\lambda e^{-\psi},~{}~{}r=a\xi=\left|\frac{A}{4\pi
G\lambda}\right|^{1/2}\xi$ (5.103)

where $\lambda\in\mathbf{R}$, then the ODE (5.98) reduces to the dimensionless
Chandrasekhar equation

$\displaystyle\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\psi(\xi)}{d\xi}\right)=e^{-\psi(\xi)}$
(5.104)

which is the Lane-Emden equation for $\gamma=1$ or $n=\infty$. One can choose
$\lambda=\rho_{c}$ with $\lim_{\xi\rightarrow 0}\frac{d\psi(0)}{d\xi}=0$ if
$\psi(\xi)$ bounded. The boundary conditions are
$\psi(0)=0,\tfrac{d\psi(\xi)}{d\xi}=0,\xi=0$. The analysis of the CE is then
similar to that of the LE equation. The CE can be transformed in various ways
$\bm{[1]}$ using

$\displaystyle X=\frac{1}{\lambda}$ (5.105)

then the CE becomes

$\displaystyle X^{4}\frac{d^{2}\psi(X)}{dX^{2}}=\exp(-\psi(X))$ (5.106)

which has a solution

$\displaystyle\psi(X)=-\log 2-2\log X$ (5.107)

The Emden transform then introduces a variable Z

$\displaystyle Z=-\psi-2\log X$ (5.108)

so that

$\displaystyle\frac{d\psi(X)}{dX}=-\frac{2}{X}-\frac{dZ}{dX}$ (5.109)
$\displaystyle\frac{d^{2}\psi(X)}{dX^{2}}=-\frac{2}{X}^{2}-\frac{d^{2}Z}{dX^{2}}$
(5.110) $\displaystyle X^{2}\frac{d^{2}Z}{dX^{2}}=e^{Z}-2=0$ (5.111)

To eliminate $X$ set $X=1/\xi=e^{\beta}$ then

$\displaystyle\frac{dZ}{d\beta}=e^{-\beta}\frac{dZ}{d\beta}$ (5.112)
$\displaystyle\frac{d^{2}Z}{d\beta^{2}}=e^{-2\beta}\left(\frac{d^{2}Z}{d\beta^{2}}-\frac{dZ}{d\beta}\right)$
(5.113)

giving

$\displaystyle\frac{d^{2}Z}{d\beta^{2}}-\frac{dZ}{d\beta}+e^{Z}-2=0$ (5.114)

##  6\. The Fundamental Hydrostatic Equilibrium Equation (FHEE) And The Euler-
Poisson Equations For A Self-Gravitating Gaseous Star

The previous derivation and applications of the FEHE were mostly heuristic. In
this section, the FEHE is now considered and derived from a more formal
mathematical perspective. The fundamental description of a self-gravitating
gaseous star are the Euler-Poisson system of nonlinear partial differential
equations for a self-gravitating gas/fluid. These are essentially the Euler
equations of fluid dynamics (with viscosity set to zero) coupled to a
Newtonian gravitational potential. The Euler-Poisson system has been studied
extensively in many works. (For example $\bm{[22-26]}$ and references their
in.)

######  Definition 6.1.

The scenario describing a self-gravitating gas/fluid described by the Euler-
Poisson equations is defined as follows.

  1. (1)

Let $\mathbf{D}\subset\mathbf{R}^{3}$. Generally, we use the functions
$\rho:[0,T)\times\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq
0},p:[0,T)\times\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq 0}$ and
$\Phi:[0,T)\times\mathbf{R}^{3}\rightarrow\mathbf{R}$ for some temporal
interval $[0,T)$, to describe the spatio-temporal distributions of mass
density, pressure and Newtonian potential of the gas/fluid. Then
$\rho=\rho(\mathbf{x},t),p=p(\mathbf{x},t),\Phi=\Phi(\mathbf{x},t)$.

  2. (2)

The self-gravitating gas is supported within a spherical region
$\mathbf{D}(\lambda(t))\equiv\mathbf{D}(t)$ defined by

$\displaystyle\mathbf{D}(\lambda(t))=\mathbf{D}(t)=supp~{}\rho(t,\bullet)=\\{(\mathbf{x}\in\mathbf{R}^{3}|\rho(\mathbf{x},t)>0,\lambda(t)\gtrless\lambda(0)=R\\}\subset\mathbf{R}^{3}$

with varying radius $\lambda(t)\gtrless\lambda(0)=R$ depending on whether the
self-gravitating gaseous sphere is collapsing or expanding, or
$\lambda(t)=\lambda(0)=R$ in hydrostatic equilibrium. Then
$\mathrm{V}(\mathbf{D}(t))=\tfrac{4}{3}\pi|\lambda(t)|^{3}$ is the volume at
time t.

  3. (3)

The (moving) boundary is $\partial\mathbf{D}(t)\subset\mathbf{D}(t)$ with
surface area $\mathrm{A}(\mathbf{D}(t))=4\pi|\lambda(t)|^{2}$. The boundary
conditions are that the pressure and density vanish here so that
$\rho(\mathbf{x},t)=p(\mathbf{x},t)=0$ for all $x\in\partial\mathbf{D}(t)$. If
$\mathbf{D}(t)\bigcup\mathbf{D}^{c}(t)=\mathbf{R}^{3}$ then the vacuum region
is defined by $\rho(\mathbf{x},t)=0$ in $\mathbf{D}^{c}(t)$ and the matter or
hydrodynamic region by $\rho(\mathbf{x},t)>0,x\in\mathbf{D}(t)$. The fluid
velocity $\mathbf{u}$ is defined only on $\mathbf{D}(t)$ by $\mathbf{u}\equiv
u^{i}:[0,T)\times\mathbf{D}(t)\rightarrow\mathbf{R}^{3}$ and there is no
definition in $\mathbf{D}^{c}(t)$. The space $\mathbf{D}(0)$ is taken to be
precompact. A diffuse boundary condition can also be defined as (ref)

$\displaystyle\lim_{(\mathbf{x},t)\rightarrow(\overline{\mathbf{x}},\overline{t})}\rho(\mathbf{x},t)\nabla_{i}p(\mathbf{x},t))=0$
(6.1)

where $(x,t)\in[0,T)\times\mathbf{D}(t)$ and
$(\overline{x},\overline{t})\in[0,T)\times\partial\mathbf{D}(t)$, where again
$\partial\mathbf{D}(t)$ is the boundary of the volume $\mathbf{D}(t)$, that
is, the moving interface between the fluid/gas matter and the exterior vacuum.
This boundary condition characterizes the gradually vanishing fluid pressure
per mass near the boundary.

  4. (4)

The Newtonian potential $\Phi$ exists on
$\mathbf{R}^{3}=\mathbf{D}(t)\bigcup\mathbf{D}^{c}(t)$, in both the matter and
vacuum regions and obeys the Poisson equation

$\displaystyle\Delta\Phi=4\pi\mathscr{G}\rho,~{}~{}~{}in~{}~{}\mathbf{R}^{3}$
(6.2)

with the general solution

$\displaystyle\Phi(\mathbf{x},t)=-\mathscr{G}{\int}_{\mathbf{D}}\frac{\rho(\mathbf{y},t)d^{3}\mathbf{y}}{|\mathbf{x}-\mathbf{y}|},~{}~{}~{}~{}for~{}(t,\mathbf{x},\mathbf{y})\in[0,T)\times\mathbf{R}^{3}$
(6.3)

  5. (5)

It is generally assumed that the gas is _isentropic_ having constant entropy
or entropy per nucleon, and obeys a barytropic pressure law of the form
$p=p(\rho)$, typically a polytropic gas law

$\displaystyle
p=\mathsf{K}|\rho|^{\gamma},~{}~{}\mathsf{x}\in\mathbf{D}(t)\in\mathbf{R}^{3}$
(6.4)

(Polytropic stars will be considered in detail in Section (5) The set
$\\{\rho(x,t),p(x,t),\Phi(x,t)\\}$ is also a solution to the Euler-Poisson
system of equations but this will be considered formally in Section (6). For a
gas in hydrostatic equilibrium $\mathbf{u}$=0 and
$\rho:\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq
0},p:\mathbf{R}^{3}\rightarrow\mathbf{R}_{\geq 0}$ and
$\Phi:\mathbf{R}^{3}\rightarrow\mathbf{R}$ so that
$\rho=\rho(\mathbf{x}),p=p(\mathbf{x}),\Phi=\Phi(\mathbf{x})$. The volume of
the support is now fixed so that
$\mathrm{V}(\mathbf{D}(t))\mathrm{V}(\mathbf{D}(0))=\tfrac{4}{3}\pi
R^{3}\equiv\tfrac{4}{3}\pi|\lambda(0)|^{3}$.

  6. (6)

For self-gravitating fluid-gas systems like stars in equilibrium at radius
$R$, it is also natural to assume spherical symmetry and coordinates so that
$\rho=\rho(r),p=p(r),\Phi=\Phi(r)$ with central values $\rho(0),p(0)$ and
boundary values $\rho(R)=p(R)=0$. The domain
$\mathbf{D}(\lambda(t))=\mathbf{D}(\lambda(0))\equiv\mathbf{D}(0)=\mathbf{B}(\mathbf{x}_{o},R)$
is now a static ball with centre $\mathbf{x}_{o}$ an radius $R$. Then

$\displaystyle\mathbf{D}(\lambda(0)=R)=\mathbf{B}(0,R)=\\{\mathbf{x},\mathbf{y}\in\mathbf{R}^{3}|d(x,y)\leq
R\\}$ (6.5)

The mass contained within a radius $r$ is or ball
$\mathbf{B}(0,r)\subset\mathbf{B}(0,R)$ is then $\mathcal{M}(r)$ and
$\mathcal{M}(R)=M$, the total mass contained within $\mathbf{B}(0,R)$ It is
defined as

$\displaystyle\mathcal{M}(r)=\int_{\mathbf{B}(0,r)}4\pi\widehat{r}^{2}\rho(\widehat{r})d{r}\equiv\int_{0}^{r}4\pi\widehat{r}^{2}\rho(\widehat{r})d\overline{r}$
(6.6)

######  Definition 6.2.

If the fluid/gas is described by a fundamental Boltzmann equation then the
continuity of mass and Navier-Stokes equations can be established as a
consequence. In the absence of gravitation they are

$\displaystyle\partial_{t}\rho+\nabla.(\rho\mathbf{u})=0~{}~{}~{}~{}x\in\mathbf{D}(t),t>0$
(6.7) $\displaystyle\rho(\partial_{t}\mathbf{u}+(\bm{u}.\nabla)\bm{u})+\nabla
p=\eta\Delta\bm{u}~{}~{}~{}~{}x\in\mathbf{D}(t),t>0$ (6.8)

and with the following boundary and initial conditions.

$\displaystyle P=\rho=0~{}on~{}\partial\mathbf{D}$ (6.9)
$\displaystyle\mathcal{U}_{\partial\mathbf{D}}=\bm{u}.\bm{n}(t)~{}on~{}\partial\mathbf{D}$
(6.10)
$\displaystyle(\rho(0,\bullet),\bm{u}(0,\bullet))=(\rho_{o},u_{o}),~{}in~{}\mathbf{D}(t)$
(6.11)

where $\mathcal{U}_{\partial\mathbf{D}}$ is the normal velocity of the moving
domain boundary $\partial\mathbf{D}(t)$ with unit normal $\bm{n}(t)$. If the
gas is now coupled to Newtonian gravitation and becomes a self-gravitating gas
then one obtains the Euler-Poisson system system of nonlinear PDEs.

$\displaystyle\partial_{t}\rho+\nabla.(\rho\mathbf{u})=0,~{}~{}x\in\mathbf{D}(t),t>0$
(6.12) $\displaystyle\rho(\partial_{t}\mathbf{u}+(\bm{u}.\nabla)\bm{u})+\nabla
p=-\rho\nabla{\Phi},~{}~{}~{}~{}x\in\mathbf{D}(t),t>0$ (6.13)
$\displaystyle\Delta{\Phi}=4\pi\mathscr{G}\rho,~{}~{}~{}~{}x\in\mathbf{D}(t),t>0$
(6.14) $\displaystyle p=\rho=0~{}on~{}\partial\mathbf{D}$ (6.15)
$\displaystyle\mathcal{U}_{\partial\mathbf{D}}=\bm{u}.\bm{n}(t)~{}on~{}\partial\mathbf{D}$
(6.16)
$\displaystyle(\rho(0,\bullet),\bm{u}(0,\bullet))=(\rho_{o},u_{o}),~{}in~{}\mathbf{D}(t)$
(6.17)

and the gas-vacuum boundary is $\mathbf{R}^{3}\char
92\relax\partial\mathbf{D}$. One can also introduce a mass function
$\mathcal{M}(r,t)$

$\displaystyle\mathcal{M}(\rho)=\int_{\mathbf{D}(t)}\rho(x,t)d\mathrm{V}(x)$
(6.18)

where
$\int_{\mathbf{B}(t)}d\mathrm{V}(x)=\int_{\mathbf{B}(t)}d^{n}x=\mathrm{V}(\mathbf{B}(t))$.

The standard notations apply so that

$\displaystyle\nabla.(\rho\bm{u})=\sum_{k=1}^{3}\frac{\partial}{\partial
x^{k}}(\rho u^{k})\equiv\sum_{k=1}^{3}\nabla_{k}(\rho u^{k}),j,k=1,2,3$ (6.19)
$\displaystyle(\bm{u}.\nabla)\bm{u}=\sum_{k=1}^{3}u^{k}\frac{\partial
u^{j}}{\partial x^{k}}\equiv\sum_{k=1}^{3}u^{k}\nabla_{k}u^{j}$ (6.20)
$\displaystyle\nabla
Z=(\nabla_{1}Z,\nabla_{2}Z,\nabla_{3}Z),~{}Z=p,\mathscr{B}$ (6.21)
$\displaystyle\Delta\Phi=\sum_{k=1}^{3}\frac{\partial^{2}}{\partial
x^{k}\partial
x^{k}}=\sum_{k=1}^{3}\partial_{k}\partial^{k}\Phi\equiv\nabla^{2}$ (6.22)

The EL system can also be expressed in various coordinate systems $\bm{[22]}$.
One can also use Lagrangian derivatives $D/Dt$ whereby

$\displaystyle\frac{D}{Dt}=\partial_{t}+\bm{u}.\nabla\equiv\partial_{t}+\sum_{k=1}^{3}u^{k}\nabla_{k}$
(6.23)

so that (6.14) and (6.15) become

$\displaystyle\frac{D\rho}{Dt}+\rho\nabla.\bm{u}=0$ (6.24)
$\displaystyle\rho\frac{Dv}{Dt}+\nabla p=-\rho\nabla\Phi$ (6.25)

The EL system of PDEs can be closed by an equation of state $p=p(\rho)$ with
$dp/d\rho>0$ and the polytropic equation can be utilised so that
$p=\mathsf{K}|\rho|^{\gamma}$ with $\gamma\in[1,2]$. We also reprise the
enthalpy variable $\mathfrak{E}(\rho)$ and the functional $\psi(\rho)$ so that
$\mathfrak{E}(\rho)=\int_{0}^{\rho}dp/\rho$ and
$\psi(\rho)=\int_{0}^{\rho}\mathfrak{E}(\rho)d\rho$. The notation
$EP_{\gamma}$ will denote an Euler-Poisson system with a polytropic gas with
index $\gamma$.

######  Definition 6.3.

A solution $\rho=\rho(x,t),\phi=\phi(x,t)$ is a compactly supported classical
solution (CSCS) iff $\forall(\rho,u,\phi)\in
C^{1}([0,T]\times\mathbf{R}^{2}),\phi(t,\bullet)\in
C^{2}(\mathbf{R}^{2}),\rho>0$,everywhere and the support of $\rho(t,\bullet)$
is compact in $\mathbf{D}(t)$ for all $t\in[0,T]$.

For any CSCS the solution of the Poisson equation (6.14) is the Newtonian
potential

$\displaystyle\phi(x,t)=-\mathscr{G}\int_{\mathbf{D}}\rho\frac{(x^{\prime},t)}{|x-x^{\prime}|}d\mathbf{V}(x^{\prime})$
(6.26)

The EP equations have been studied quite extensively and there are expanding
and collapsing solutions.

######  Definition 6.4.

The total mass $\mathcal{M}(\rho)$ was defined. The total energy at any time
$t>0$ is defined as

$\displaystyle\mathbf{E}(\rho,\bm{u})=$
$\displaystyle\int_{\mathbf{D}(t)}\frac{1}{2}\rho|\bm{u}|^{2}d\mathrm{V}(x)+\int_{\mathbf{D}(t)}\psi(\rho)d\mathrm{V}(x)+\int_{\mathbf{D}(t)}\frac{1}{2}\rho\phi
d\mathrm{V}(x)$
$\displaystyle=\int_{\mathbf{D}(t)}\left(\frac{1}{2}\rho|\bm{u}|^{2}d\mathrm{V}(x)+\psi(\rho)\right)d\mathrm{V}(x)-\frac{1}{2}\mathscr{G}\int\int\frac{\rho(x,t)\rho(x^{\prime},t)}{|x-x^{\prime}|}d\mathrm{V}(x)d\mathrm{V}(x^{\prime}$
(6.27)

For a polytropic gas

$\displaystyle\mathbf{E}(\rho,\bm{u})=$
$\displaystyle\int_{\mathbf{D}(t)}\frac{1}{2}\rho|\bm{u}|^{2}d\mathrm{V}(x)+\int_{\mathbf{D}(t)}\frac{\mathsf{K}}{\gamma-1}\rho^{\gamma}d\mathrm{V}(x)+\int_{\mathbf{D}(t)}\frac{1}{2}\rho\phi
d\mathrm{V}(x)$
$\displaystyle=\int_{\mathbf{D}(t)}\left(\frac{1}{2}\rho|\bm{u}|^{2}d\mathrm{V}(x)+\frac{\mathsf{K}}{\gamma-1}\rho^{\gamma}\right)d\mathrm{V}(x)-\frac{\mathscr{G}}{2}\int\int\frac{\rho(x,t)\rho(x^{\prime},t)}{|x-x^{\prime}|}d\mathrm{V}(x)d\mu(x^{\prime}$
$\displaystyle=\int_{\mathbf{D}(t)}\left(\frac{1}{2}\rho|\bm{u}|^{2}d\mu(x)+\frac{p}{\gamma-1}\right)d\mathrm{V}(x)-\frac{\mathscr{G}}{2}\int\\!\\!\\!\int\frac{\rho(x,t)\rho(x^{\prime},t)}{|x-x^{\prime}|}d\mu(x)d\mathrm{V}(x^{\prime})$
(6.28)

######  Theorem 6.5.

The total energy $\mathbf{E}(\rho,\bm{u})$ is conserved for any solution
$(\rho,\bm{u})$ so that

$\displaystyle\frac{d}{dt}\mathbf{E}(\rho,\bm{u})=0$ (6.29)

###### Proof.

$\displaystyle\frac{d}{dt}\mathbf{E}(\rho,\bm{u})$
$\displaystyle=\frac{d}{dt}\int\frac{1}{2}\rho\bm{u}|^{2}d\mathrm{V}=\frac{1}{2}\int\partial_{t}(\rho\sum_{k}|u_{k}|^{2}d\mathrm{V}$
(6.30)
$\displaystyle=\int\rho\sum_{k}u^{k}\partial_{t}u^{k}d\mathrm{V}+\int\frac{1}{2}\partial_{t}\rho\sum_{k}|u_{k}|^{2}d\mathrm{V}$
$\displaystyle=-\int(\frac{1}{2}\rho\sum_{j,k}u^{j}\nabla)_{j}|u^{k}|^{2}d\mathrm{V}+\int\sum_{k}u^{k}\nabla_{k}pd\mathrm{V}$
$\displaystyle+\int\sum_{k}\rho u^{k}\nabla_{k}\Phi
d\mathrm{V}+\int\frac{1}{2}\nabla_{j}|\rho
u^{j}|\sum_{k}|u_{k}|^{2}d\mathrm{V}$
$\displaystyle=-\int(\frac{1}{2}\rho\sum_{j,k}|u^{j}\nabla_{j}|u^{k}|^{2}d\mathrm{V}+\int\sum_{k}u^{k}\nabla_{k}pd\mathrm{V}$
$\displaystyle+\int\sum_{k}\rho
u^{k}\nabla_{k}\Phi+\int\frac{1}{2}\nabla_{j}|\rho
u^{j})\sum_{k}|u^{k}|^{2}d\mathrm{V}$
$\displaystyle=-\underbrace{\int\frac{1}{2}\sum_{j,k}\nabla_{j}(\rho
u^{j}|u^{k}|^{2})-\int\sum_{k}u^{k}(\nabla_{k}p+\rho\nabla_{k}\Phi)d\mathrm{V}}$
$\displaystyle=-\int\sum_{k}\rho
u^{k}\nabla_{k}\mathfrak{E}d\mathrm{V}+\int\sum_{k}\rho u^{k}\nabla_{k}\Phi
d\mathrm{V}$ $\displaystyle=\int\frac{\partial\Psi}{\partial
t}d\mathrm{V}+\int\underbrace{\sum_{k}\nabla_{k}(\rho
u^{k})}_{use~{}EP~{}eqn~{}again}\Phi d\mathrm{V}$
$\displaystyle=-\frac{d}{dt}\int\Psi
d\mathrm{V}-\int\frac{\partial\rho}{\partial t}\Phi d\mathrm{V}$ (6.31)

The last term can be re-expressed as

$\displaystyle\int\partial_{t}\rho\Phi d\mathrm{V}$
$\displaystyle=-\mathscr{G}\int\\!\\!\\!\int\frac{\partial_{t}\rho(x,t)\rho(x,t)}{|x-x^{\prime}|}d\mathrm{V}(x)d\mathrm{V}(x^{\prime})$
$\displaystyle=-\frac{1}{2}\frac{d}{dt}\int\\!\\!\\!\int\frac{\rho(x,t)\rho(x,t)}{|x-x^{\prime}|}d\mathrm{V}(x)d\mathrm{V}(x^{\prime})=\frac{1}{2}\frac{d}{dt}\int\rho\Phi
d\mathrm{V}$ (6.32)

Hence

$\displaystyle\frac{d}{dt}\int\frac{1}{2}\rho|u|^{2}d\mathrm{V}+\frac{d}{dt}\int\Psi
d\mathrm{V}+\frac{d}{dt}\int\frac{1}{2}\rho\Phi d\mathrm{V}=0$ (6.33)

and so the total energy is conserved along t for any compact supported
solution $(\rho,\bm{u})$, ∎

Solutions of an $EP_{\gamma}$ system have a self-similar re-scaling property
under $\overline{t}=t/\lambda^{1/(2-\gamma)}$ and $\overline{x}=x/\lambda$. If
$(\rho,\bm{u})$ is a solution of an $EP_{\gamma}$ system then so is the pair
$(\overline{\rho},\overline{\bm{u}})$ given by $\bm{[32,36]}$.

$\displaystyle{\rho(x,t)}=\lambda^{-\frac{2}{2-\gamma}}\overline{\rho}(\frac{t}{\lambda^{\frac{1}{2-\gamma}}},\frac{x}{\lambda})$
(6.34)
$\displaystyle{\bm{u}(x,t)}=\lambda^{-\frac{\gamma-1}{2-\gamma}}\overline{\bm{u}}(\frac{t}{\lambda^{\frac{1}{2-\gamma}}},\frac{x}{\lambda})$
(6.35)

Also, given a solution $(p,\phi)$ for the pressure and Newtonian potential,
the rescaled pressure and potential are

$\displaystyle{p(x,t)}=\lambda^{-\frac{2\gamma}{2-\gamma}}\overline{o}(\frac{t}{\lambda^{\frac{1}{2-\gamma}}},\frac{x}{\lambda})$
(6.36)
$\displaystyle{\phi(x,t)}=\lambda^{-\frac{2\gamma-2}{2-\gamma}}\overline{\phi}(\frac{t}{\lambda^{\frac{1}{2-\gamma}}},\frac{x}{\lambda})$
(6.37)

The mass $M(\rho)$ and the energy $E(\rho,u)$ are rescaled as

$\displaystyle
M(\rho)=\lambda^{\frac{4-3\gamma}{2-\gamma}}\overline{M(\overline{\rho})}$
(6.38) $\displaystyle
E(\rho,\bm{u})=\lambda^{\frac{6-5\gamma}{2-\gamma}}\overline{E(\overline{\rho},\overline{\bm{u}})}$
(6.39)

The mass is invariant under a re-scaling when $\gamma=4/3$ and the energy when
$\gamma=5/6$. so the re-scalings are respectively mass-critical and energy
critical.

The $EP_{\gamma}$ system applied to a spherical gaseous star is naturally best
expressed in spherical symmetry then $r=\|x\|$. Define a ball