Fluid statics of a self-gravitating perfect-gas isothermal sphere

TL;DR

This paper investigates the fluid static equilibrium of an isothermal, self-gravitating perfect gas sphere using a novel boundary condition approach, revealing multiple solutions and thermodynamic stability considerations.

Contribution

It introduces a new method for solving the isothermal sphere problem by fixing the gravitational field at the shell boundary, differing from traditional central density prescriptions.

Findings

Identification of a gravitational number upper bound.

Existence of multiple equilibrium solutions.

Preliminary thermodynamic stability analysis.

Abstract

We open the paper with introductory considerations describing the motivations of our long-term research plan targeting gravitomagnetism, illustrating the fluid-dynamics numerical test case selected for that purpose, that is, a perfect-gas sphere contained in a solid shell located in empty space sufficiently away from other masses, and defining the main objective of this study: the determination of the gravitofluid-static field required as initial field ($t=0$) in forthcoming fluid-dynamics calculations. The determination of the gravitofluid-static field requires the solution of the isothermal-sphere Lane-Emden equation. We do not follow the habitual approach of the literature based on the prescription of the central density as boundary condition; we impose the gravitational field at the solid-shell internal wall. As the discourse develops, we point out differences and similarities between the literature's and our approach. We show that the nondimensional formulation of the problem hinges on a unique physical characteristic number that we call gravitational number because it gauges the self-gravity effects on the gas' fluid statics. We illustrate and discuss numerical results; some peculiarities, such as gravitational-number upper bound and multiple solutions, lead us to investigate the thermodynamics of the physical system, particularly entropy and energy, and preliminarily explore whether or not thermodynamic-stability reasons could provide justification for either selection or exclusion of multiple solutions. We close the paper with a summary of the present study in which we draw conclusions and describe future work.

Figures

Tables

Table 1. Table 1: physical meaning of partial derivatives

	gas density (logarithmic)	$N\times \left(\text{<span> <span>gravitational</span> <span>field</span> </span>}\right)$
$\frac{\partial \ln \text{x}}{\partial \text{h}}$	slope
$\frac{{\partial }^2\ln \text{x}}{\partial {\text{h}}^2}$	curvature	slope
$\frac{{\partial }^3\ln \text{x}}{\partial {\text{h}}^3}$		curvature

Equations

$$\nabla^{\mbox{\tiny$2$}}\mbox{{\tengr y}}=\begin{cases}\;4\mbox{{\tengr p}}G\mbox{{\tengr r}}\quad&\text{gas}\\[10.76385pt] \;4\mbox{{\tengr p}}G\mbox{{\tengr r}}_{s}&\text{shell}\\[10.76385pt] \;0&\text{empty space}\end{cases}$$

$$\mbox{\boldmath g}=-\mbox{\boldmath$\nabla$}\mbox{{\tengr y}}$$

$$\frac{\partial \text{mbox}\text{tengr}r}{\partial t}$$

$$\displaystyle\mbox{\boldmath$\nabla$}p-\mbox{{\tengr r}}\mbox{\boldmath g}$$

$$\displaystyle\mbox{{\tengr r}}c_{v}\frac{\partial T}{\partial t}-\mbox{\boldmath$\nabla$}\cdot\left(\mbox{\textgr{l}}\mbox{\boldmath$\nabla$}T\right)$$

$$p=\text{mbox}\text{tengr}rRT$$

$$\mbox{\boldmath J}_{\mbox{\tiny$\mbox{U}$}}=-\mbox{\textgr{l}}\mbox{\boldmath$\nabla$}T$$

$$\mbox{\boldmath$\nabla$}\cdot\left(\mbox{\textgr{l}}\mbox{\boldmath$\nabla$}T\right)=0$$

$$\nabla^{\mbox{\tiny$2$}}=\frac{1}{r^{2}}\frac{\partial}{\partial r}\left(r^{2}\frac{\partial}{\partial r}\right)=\frac{\partial^{2}}{\partial r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}$$

$${\frac{\partial }{\partial r}|}_{r=0}=0$$

$$\underset{r\to 0}{\lim }\frac{1}{r}\frac{\partial }{\partial r}={\frac{{\partial }^2}{\partial r^2}|}_{r=0}$$

$$\lim_{r\rightarrow 0}\;\nabla^{\mbox{\tiny$2$}}=3{\left.\frac{\partial^{2}\ }{\partial r^{2}}\right|}_{r=0}$$

$$\frac{{\partial }^2\text{mbox}\text{tengr}y}{\partial r^2}+\frac{2}{r}\frac{\partial \text{mbox}\text{tengr}y}{\partial r}=4\text{mbox}\text{tengr}pG\text{mbox}\text{tengr}r$$

$${\frac{\partial \text{mbox}\text{tengr}y}{\partial r}|}_{r=0}=-g(0)=0$$

$$g(a^{\text{mbox}-})=g(a^{\text{mbox}+})$$

$$g[(a+\text{mbox}\text{tengr}d{)}^{\text{mbox}-}]=g[(a+\text{mbox}\text{tengr}d{)}^{\text{mbox}+}]$$

$$g(\infty )=0$$

$$\begin{gathered} g(r)=\{\begin{array}{ll}-\frac{G}{r^2}\left(4\text{mbox}\text{tengr}p{\int }_0^r\text{mbox}\text{tengr}r(x)x^{2}dx\right)& \text{gas}\\ -\frac{G{m}_g}{r^2}\left(1+\frac{m_s}{m_g}\frac{r^3-a^3}{(a+\text{mbox}\text{tengr}d{)}^3-a^3}\right)& \text{shell}\\ -\frac{G{m}_g}{r^2}\left(1+\frac{m_s}{m_g}\right)& \text{parbox}20.00003pt\text{centering}empty \\ space\text{@add@centering}\end{array} \end{gathered}$$

$$m_g$$

$$m_s$$

$$\overline{\text{mbox}\text{tengr}r}=\frac{m_g}{V}=\frac{m_g}{\frac{4}{3}\text{mbox}\text{tengr}p{a}^3}$$

$$\bar{p}=\overline{\text{mbox}\text{tengr}r}RT$$

$$g(0)=-\underset{r\to 0}{\lim }\frac{G}{r^2}\cdot \left[4\text{mbox}\text{tengr}p{\int }_0^r\text{mbox}\text{tengr}r(x)x^{2}dx\right]=0$$

$$g(a)=-\frac{G{m}_g}{a^2}$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(\text{mbox}\text{textgr}l{r}^2\frac{\partial T}{\partial r}\right)=0$$

$$\text{mbox}\text{textgr}l{r}^2\frac{\partial T}{\partial r}={\left(\text{mbox}\text{textgr}l{r}^2\frac{\partial T}{\partial r}\right)}_{r=0}=0$$

$$\frac{\partial T}{\partial r}=0$$

$$T(a)=T_s$$

$$T(r)=T_s$$

$$p=\text{mbox}\text{tengr}rR{T}_s$$

$$\mbox{\boldmath g}=RT\mbox{\boldmath$\nabla$}\ln\mbox{{\tengr r}}$$

$$g=RT\frac{\partial \ln \text{mbox}\text{tengr}r}{\partial r}$$

$$\mbox{\boldmath$\nabla$}\left(\mbox{{\tengr y}}+RT\ln\mbox{{\tengr r}}\right)=0$$

$$\text{mbox}\text{tengr}y=A-RT\ln \text{mbox}\text{tengr}r$$

Full text

1 Introduction

The motivation and the inspiration for the research activities that we begin to describe in this paper originate in a somewhat peculiar manner from aerothermodynamics, a physical/engineering discipline rather distant from gravitation. In 2001, the aerothermodynamics section of the European Space Agency initiated a series of studies driven by interest in investigating the influence of electromagnetic fields on the heat transfer to a body in hypersonic flow. The studies addressed theoretical, numerical and experimental aspects of the complex physics encompassed in electromagnetic fluid dynamics; the most challenging part was undoubtedly the one aiming at the development of numerical algorithms for the solution of the fully coupled Maxwell equations and Navier-Stokes equations, the latter with due inclusion of the physics characterizing flows in strong thermochemical nonequilibrium. The coupling is conceptually smooth and straightforward from a theoretical point of view; on the other hand, coupling the two sets of equations from a numerical point of view is a different story in consequence of the huge disparity in orders of magnitude between the speeds of light, on the electromagnetism side, and sound, on the fluid-dynamics side. As a matter of fact, the numerical challenge turned out to be formidable and it is still unbeaten nowadays, many assaults by practitioners of computational fluid dynamics notwithstanding. The core difficulty lies mainly in the diffusion and relaxation processes taking place in thermochemical nonequilibrium that add numerical heaviness to the field equations with terms of such different orders of magnitude as to make their interplay practically uncontrollable, at least until now, by all attempted numerical-algorithm machineries.

How about gravitational fields? Terrestrial fluid dynamicists are prevalently occupied with applications on the surface of the planet; for them, the gravitational field is externally imposed, uniform and, above all, known because the contribution to the field from the fluid mass under consideration and from the other masses in the universe is absolutely negligible with respect to the planet’s contribution. Astrophysical fluid dynamicists dealing with non-relativistic applications resort to Newton’s theory of gravity to study the dynamics of (self-)gravitating fluid masses. Yet, Newton’s theory of gravity is an action-at-distance theory; its incompatibility with the relativity theory has been overwhelmingly detected, addressed and resolved by well-known eminent physicists. More humbly, a fluid dynamicist could add that Newton’s theory of gravity smoothly matches situations belonging to fluid statics and steady-state fluid dynamics but presents conceptual rugosities when confronted with time-dependent fluid-dynamics circumstances. Straightforward analysis suggests that a theory of gravity based on a single-vector field, like Newton’s, is not sufficient when the non-relativistic dynamics of a (self-)gravitating fluid mass is unsteady; two vector fields are necessary and they are governed by differential equations whose mathematical structure coincide with the Maxwell equations in which densities of matter’s mass and linear momentum play the role of field sources [1, 2, 3, 4, 5]. Interestingly enough, the bottom-up conclusions achieved in a fluid-dynamics context are confirmed by a top-down approach that starts from Einstein’s theory of gravity and, in the weak but time-dependent field approximation, leads into the realm of the fiercely debated subject known as gravitomagnetism [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44], the battleground in which scientists struggle to chase experimentally the proof of existence of the Lorentz-type force for gravity.

Gravitational fluid dynamics offers, therefore, the potential opportunity to skirmish with the numerical difficulties of the Maxwell-Navier-Stokes system of differential equations but without the numerical heaviness deriving from the complex physics of flows in thermochemical nonequilibrium; a perfect-gas model in standard conditions of pressure and temperature will do!

This is the challenge we settled to deal with and the test case we selected for this purpose is illustrated in Fig. 1. A perfect gas is contained in a spherical solid shell of radius $a$ and thickness d placed in empty space far away from other masses in the universe. The gas is self-gravitating and we imagine to induce non-relativistic motion in it by adequately changing the internal-wall temperature distribution, in an axisymmetric manner to conveniently simplify the flow pattern. We intend to calculate the ensuing flow field with both theories of gravity: Newton’s and gravitomagnetic Maxwell’s.

The main objective of this paper is the determination of the initial fluid and gravitational fields ($t=0$) required by the time-dependent calculation of the dynamics ($t>0$). We assume the gas initially quiescent [$\mbox{\boldmath v}(\mbox{\boldmath r},0)=0$]; then, the pursuit of the initial fields leads straightforwardly to the kind of Lane-Emden equations, a well-known family of differential equations recurrent in astrophysics [45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65] and in terrestrial applications such as fuel ignition [66] and thermal explosions [67, 68], that describes the so-called isothermal-gas sphere. The main difference between our manner to deal with such a differential equation and the approach followed, more or less systematically, in the astrophysical community resides in the boundary conditions. In our problem, gas confinement rules and the gravitational side of the physics provides the required boundary condition at the shell internal wall ($r=a$); there is no need to, and in fact we do not, impose the gas density either at the center or at any other radial position in the gas sphere. In the Lane-Emden problem, there is no gas confinement and a general consensus exists supporting the idea that the imposition of the density at the center of the gas sphere constitutes a legitimate boundary condition. We have grappled with such an idea for a long time but have not been able to reconcile its presumed legitimacy with the arbitrariness descending from the undeniable lack of any independent physical information fixing the density at the center of the gas sphere.

The paper is logically organized in two parts. In the first part (Sec. 2), we describe the solution procedure to obtain the initial fields for our test case; we also briefly review the salient aspects of the Lane-Emden solution for the isothermal gas sphere and emphasize differences/equivalences with our approach (Sec. 2.5). In the second part (Sec. 3), we concentrate on the thermodynamics of the physical system and discuss entropy, energy and aspects of thermodynamic stability.

2 Gravitofluid-static fields

2.1 Governing equations and boundary conditions

Let us consider a prescribed mass $m_{g}$ of perfect gas in mechanical equilibrium [$\mbox{\boldmath v}(\mbox{\boldmath r},t)=0$] contained in the spherical solid shell illustrated in Fig. 1. The shell has internal radius $a$, contains a volume $V=4/3\,\mbox{{\tengr p}}a^{3}$ and has thickness d. The gravitational field occupies all space and is governed by

[TABLE]

according to Newton’s theory of gravity. The fluid-static field is governed by the static versions of the traditional balance equations of mass, momentum and energy

[TABLE]

together with the thermodynamic state equation

[TABLE]

In Eqs. (1)–(3),

[TABLE]

The continuity equation [Eq. (2a)] enforces the time independence of the gas density. The mechanical equilibrium exists also in the shell and, therefore, its density must also be time independent $(\partial\mbox{{\tengr r}}_{s}/\partial t=0)$; moreover, we consider negligible the deformation field and assume the shell density uniform and prescribed. Consequently, time independence is sequentially passed on to gravitational potential and field by Eqs. (1a) and (1b), to pressure by the momentum equation [Eq. (2b)], and to temperature by the state equation [Eq. (3)]. The temperature time derivative in the energy equation [Eq. (2c)] vanishes identically (we left it in there only for formal consistency and clarity) and the heat flux, assumed according to Fourier law

[TABLE]

becomes necessarily solenoidal

[TABLE]

The mechanical equilibrium of the gas implies temperature uniformity in the containing shell; otherwise, temperature gradients would settle in at the gas/solid wall and would induce motion in the gas through the pressure gradient in the momentum equation. An important consequence ensuing from the shell-temperature uniformity is the spherical symmetry of the physical system; hence, we conveniently select standard polar coordinates $r,\mbox{{\tengr j}},\mbox{{\tengr f}}$ with origin coincident with the center of the spherical geometry (Fig. 1). Before further processing the governing equations, it is convenient to pay a bit of attention to the spherically-symmetric Laplace operator

[TABLE]

because it may become unbounded when $r\rightarrow 0$. Clearly, we must assume

[TABLE]

in order to avoid the physically unacceptable discontinuity; then, the second term on the rightmost-hand side of Eq. (6) becomes an indeterminate form that can be cured by de l’Hôpital’s theorem

[TABLE]

and the spherically-symmetric Laplace operator [Eq. (6)] becomes discontinuity-free

[TABLE]

The differential equations [Eqs. (1)–(2)] governing the statics of the physical system, composed by gas, shell and gravitational field, must be supplemented with the required boundary conditions. Let us begin with the gravitational field which, of course, has only the radial component. The expansion of Eq. (1a) (gas)

[TABLE]

provides, according to Eq. (7), the boundary condition

[TABLE]

The spherical symmetry of the physical system allows to shortcut the integration of Eqs. (1) by taking advantage of Gauss’ theorem, whose application to concentric spherical surfaces inside the gas ($r<a$), inside ($a<r<a+\mbox{{\tengr d}}$) and outside ($a+\mbox{{\tengr d}}<r$) the shell yields respectively

[TABLE]

The term in parentheses in Eq. (12) (gas) is the mass of gas contained in the sphere of radius $r\;(<a)$; the variable $x$ represents a dummy integration variable. The masses of gas and shell in Eq. (12) (shell & empty space) are expressible in terms of the corresponding densities

[TABLE]

Equation (13a) should be read from right to left: it is a constraint to which the density distribution must comply because the mass on its left-hand side is prescribed. It is also convenient at this point to introduce average density

[TABLE]

of the gas in view of the forthcoming nondimensional analysis of Sec. 2.2. Compliance of the gravitational-field solution [Eq. (12)] with the gravitational boundary conditions [Eq. (11)] is easily verified. The verification of Eq. (11a) calls for a bit of attention because the limit of Eq. (12) (gas)

[TABLE]

requires repeated application of de l’Hôpital’s theorem but verification of field-continuity and asymptotic-vanishing boundary conditions [Eqs. (11b)–(11d)] is straightforward. Of course, Eq. (12) (gas) becomes operational only when the density distribution is known; however, it nails down the gravitational field at the shell internal wall to the immutable value

[TABLE]

regardless of the specific density distribution that settles in the gas. Equation (16) plays an important role for the determination of the fluid-static field to which we move on now.

The energy equation [Eq. (5)] reduces to the spherically-symmetric form

[TABLE]

that can be easily integrated starting from $r=0$

[TABLE]

to yield

[TABLE]

According to Eq. (19), the perfect gas cannot sustain thermal gradients in a spherically-symmetric steady state and must necessarily be isothermal. The thermal boundary condition must be prescribed at the shell internal wall because, in compliance with shell-temperature uniformity, the gas temperature there must coincide with the shell’s

[TABLE]

regardless of the angular position. From the imposition of Eq. (20) we obtain

[TABLE]

and the state equation [Eq. (3)] simplifies to the isothermal form

[TABLE]

Hereinafter, we will drop the subscript ‘$s$’ from the temperature symbol to simplify the notation.

The next step consists in introducing the isothermal state equation [Eq. (22)] into the momentum equation [Eq. (2b)] and solving the latter for the gravitational field

[TABLE]

Given the spherical symmetry, Eq. (23) is equivalent to the scalar form

[TABLE]

By comparing Eq. (12) (gas) and Eq. (24), the attentive reader may wonder whether or not those differently looking expressions of the gravitational field are equivalent; we provide the answer at the end of this section.

The comparison of Eq. (23) with Eq. (1b) leads to

[TABLE]

from which we deduce the gravitational potential

[TABLE]

in terms of the gas density, save for an arbitrary and inessential constant $A$. Then, the substitution of Eq. (26) into Eq. (1a) (gas) leads to the isothermal Lane-Emden equation

[TABLE]

This differential equation must be integrated with the gravitational boundary conditions [Eqs. (11a) and (16)] reformulated, by account of Eq. (24), in terms of the gas-density radial gradient

[TABLE]

The mathematical problem [Eqs. (27b)–(28)] is thus well posed: we have a second-order differential equation and its two required boundary conditions whose physical meaning is clear and unambiguous. It is not necessary to think the density as assigned either at the sphere center or at any other radial position; as a matter of fact, the value $\mbox{{\tengr r}}(0)$ can, actually must, be obtained as a result after that the density distribution has been determined from the integration of Eq. (27b). With that distribution in hand, the pressure distribution follows straightforwardly from the isothermal state equation [Eq. (22)].

Before proceeding to cast the mathematical problem in nondimensional form, we believe a few brief considerations regarding physical consistency are in order. Multiplication of Eq. (27b) by $r^{2}$, subsequent integration from the sphere center up to the generic radial position $r$, and slight rearrangement of constant factors give

[TABLE]

Equation (29) shows the equivalence between the integral expression [Eq. (12) (gas)] originated from Gauss’ theorem and the differential expression [Eq. (24)] of the gravitational field that follows from the fluid-statics momentum equation [Eq. (2b)] of the isothermal gas. This is not surprising: the equivalence was established when we substituted Eq. (26) into Eq. (1a) (gas). Equation (29) also clearly indicates that $\partial\ln\mbox{{\tengr r}}/\partial r<0$ as a consequence of the attractive nature of the gravitational field; therefore, density and pressure [Eq. (22)] are monotonically decreasing from the sphere center ($r=0$) to the shell internal wall ($r=a$)

[TABLE]

If calculated at the shell internal wall, Eq. (29) gives

[TABLE]

Substitution of Eq. (28b) into Eq. (31) and simplification of constant factors yield

[TABLE]

Equation (32) indicates that the gas-mass constraint [Eq. (13a)] is identically satisfied; in other words, the imposition of the second boundary condition [Eq. (28b)] ensures that whatever density distribution be extracted from Eq. (27b), it complies automatically with the prescribed mass of the gas, a somewhat comforting result.

2.2 Nondimensional formulation

According to standard practice, we introduce nondimensional radial coordinate and density together with corresponding dimensional scale factors marked with a tilde

[TABLE]

in order to cast the mathematical problem [Eqs. (27b)–(28)] in nondimensional form

[TABLE]

The set of Eqs. (34) features three characteristic numbers

[TABLE]

In general, the scale factors are chosen according to, hopefully judicious, criteria of convenience and the P-type characteristic numbers are deduced therefrom. A more preferable approach is the other way around: it is more convenient to resolve the scale factors

[TABLE]

in terms of the P-type characteristic numbers and to set the latter to convenient values. The second approach has two indisputable advantages. The first is that the dependent P-type characteristic numbers are disclosed with clear evidence; indeed, Eq. (36) reveals that $\mbox{{\tengr P}}_{1}$ is unrestricted but $\mbox{{\tengr P}}_{2}$ and $\mbox{{\tengr P}}_{4}$ cannot be assigned independently. True, this could have been detected easily also from attentive scrutiny of Eq. (35) but just because we are dealing here only with two scale factors and three P-type characteristic numbers. If the set of differential equations and boundary conditions becomes overcrowded with scale factors and characteristic numbers then processing the analog of Eq. (35) for the purpose of setting the scale factors to convenient values turns into an unmanageable task prone to unavoidable mistakes. The second advantage of the second approach is that it brings to light the physical characteristic numbers that control the problem at hand. In this case, Eq. (36) displays only (the one we call) the gravitational number

[TABLE]

With the scale factors enforced by Eq. (36), the nondimensional mathematical problem [Eqs. (34)] can be rephrased accordingly

[TABLE]

and the mass constraint [Eq. (13a)] turns into the normalization condition

[TABLE]

Remembering the considerations at the end of Sec. 2.1 [Eqs. (31) and (32)], we should expect Eq. (39) to be identically satisfied due to the imposition of Eq. (38c).

The next step consists in choosing the two independent P-type characteristic numbers. In view of the necessarily numerical solution of Eq. (38a), we set $\mbox{{\tengr P}}_{4}=1$ to keep the computational domain for h fixed to the interval [0,1], a desirable feature when doing numerical calculations; we then set $\mbox{{\tengr P}}_{1}=3N$ to remove the dependence of the density scale factor on the gravitational number. With these choices, Eq. (36) simplifies to

[TABLE]

and Eqs. (33) yield the nondimensional variables

[TABLE]

that were also selected by Darwin [50]; the mathematical problem [Eqs. (38)] attains its final form

[TABLE]

conveniently predisposed for numerical processing and the accompanying normalization condition [Eq. (39)] becomes

[TABLE]

Obviously, every solution of Eqs. (42) depends parametrically on the gravitational number

[TABLE]

but we will explicitly indicate that dependence only if required by the context. We wish to point out what we believe to be fine features of the nondimensional problem [Eqs. (41) and Eqs. (42)] based on Eq. (40): it includes boundary conditions whose physical legitimization is unambiguous; it is governed by one single, clearly identified physical characteristic number, the gravitational number $N$; it involves nondimensional variables scaled with factors constructed with variables characterizing the physical system, known a priori and as precisely controllable as the gravitational number is.

To complete the nondimensional formulation, we scale the gravitational field [Eqs. (12) and (24)] with the value [Eq. (16)] attained at the shell internal wall

[TABLE]

[TABLE]

precisely as Ritter [47] and Darwin [50] also did, and the pressure with its average value [Eq. (22)]

[TABLE]

2.3 Some analytical results

It is somewhat interesting to analyze the qualitative behavior of the solutions of Eqs. (42) before proceeding to discuss numerical operations and results.

Given the selection $\mbox{{\tengr P}}_{1}=3N,\mbox{{\tengr P}}_{4}=1$, the inequalities shown in Eqs. (30) merge into

[TABLE]

Multiplication of Eq. (48) by $3\,\mbox{{\tengr h}}^{2}$ and integration on the interval [0,1], with due account of the normalization condition [Eq. (43)], leads to the interesting result

[TABLE]

Thus, the density is always greater or lower than its average value at, respectively, the sphere center or the shell internal wall; therefore, a radial position $\bar{\mbox{{\tengr h}}}$ must exists at which the density attains its average value [Eq. (14a)]

[TABLE]

and the conceptual boundary between a denser (than average) core ($0\leq\mbox{{\tengr h}}<\bar{\mbox{{\tengr h}}}$) and a less dense peripheral layer ($\bar{\mbox{{\tengr h}}}<\mbox{{\tengr h}}\leq 1$) can be set.

When $N\rightarrow\epsilon\ll 1$, Eqs. (42) provide infinite solutions

[TABLE]

The normalization condition [Eq. (43)] fixes the value of the constant to $C=1$ and singles out the only physical solution

[TABLE]

Therefore, density and pressure are uniform throughout the gas sphere and the self-gravitating condition becomes negligible, as expected. The gravitational field in the gas sphere is easily obtained from Eq. (45) (gas) with $\mbox{{\tengr x}}(u)\simeq 1$ and turns out to decrease linearly from the sphere center to the shell internal wall

[TABLE]

consistently with the existence of a uniform density. Taking into account Eq. (53), we also deduce from the alternative equivalent expression [Eq. (46)] of the gravitational field that

[TABLE]

in view of forthcoming considerations.

Considering $N\neq\epsilon\ll 1$ now, more qualitative information can be gathered from the evaluation of the logarithmic-density derivatives, up to third order, evaluated at the boundary points of the interval [0,1]. The physical meaning of these derivatives is summarized for convenience in Table 1; the third-order derivative

[TABLE]

is easily obtained from further derivation of Eq. (42a). The situation at the sphere center

[TABLE]

First and second derivatives [Eqs. (56a) and (56b)] inform that the density attains a maximum there in consequence of the vanishing of the gravitational field. The latter’s profile has negative slope and infinite curvature [Eq. (56c); evaluation requires application of de l’Hôpital’s theorem to first term on right-hand side of Eq. (55)]. The situation at the shell internal wall

[TABLE]

is more informative. According to Eqs. (56b) and (57b), the function $\ln\mbox{{\tengr x}}$ preserves its negative curvature on the whole interval [0,1] up to a value $N^{\prime}$ of the gravitational number at which

[TABLE]

If $N=N^{\prime}$ then the second derivative [Eq. (57b)] vanishes; hence, an inflection point and a minimum appear for, respectively, logarithmic density and gravitational-field profiles. The third derivative [Eq. (57c)] takes on the value

[TABLE]

and its necessary positivity (gravitational-field minimum) leads to $N^{\prime}>1$. If $N^{\prime}<N$ then the positions of both inflection point and minimum shift leftward inside the interval [0,1]. Another noticeable value $N^{\prime\prime}$ of the gravitational number is the one at which

[TABLE]

If $N=N^{\prime\prime}$ then the third derivative [Eq. (57c)] vanishes and an inflection point appears for the gravitational-field profile.

2.4 Numerical solution and results

As repeatedly mentioned and well emphasized in the literature, the differential equation [Eq. (42a)] governing the mathematical problem is not amenable to analytical integration and has to be dealt with numerically. In this regard, we have used three different numerical algorithms (labels F, FP, D in the forthcoming figures) independently coded in three different programming languages (R, matlab and fortran). The first algorithm (F) is based on the standard approach of recasting the second-order differential equation [Eq. (42a)] as a first-order system and then applying a general-purpose code for boundary-value problems [70, 71, 72, 73, 74]. The second algorithm (FP) [75, 76] is based on high-order (up to eight) finite-difference schemes to solve general second-order ordinary differential equations subjected to Neumann, Dirichlet or mixed boundary conditions; the ideas underlying this approach have also been adapted to solve Sturm-Liouville problems with one and two parameters [77, 78]. The third algorithm (D) is based on a finite-difference scheme whose nonlinear algebraic system is treated with a multidimensional Newton-Raphson iterative method. All the algorithms need initial guesses of the solution from which the iterative procedure starts and proceeds to the converged solution.

The agreement among the results produced by the different algorithms is very satisfactory, as evidenced in Figs. 2 and 3; unless otherwise specified, the data-representation convention displayed in the legends (solid line for algorithms F & FP; solid circles for algorithm D) applies systematically in all figures.

Figure 2 shows the gas-density profiles for selected values of the gravitational number, in logarithmic scale to illustrate more clearly the situation in the vicinity of the shell internal wall. In view of the discussion about the Lane-Emden solution (Sec. 2.5), we call the reader’s attention to an important detail: the density values attained at sphere center and shell internal wall are unequivocally determined by the integration of Eqs. (42); no knowledge of those values is required before integration. Figure 3 shows the corresponding profiles of the gravitational field inside the gas sphere. All the analytical characteristics described in Sec. 2.3 are exhibited by the numerical results illustrated in Figs. 2 and 3. The density profiles feature the expected monotonic behavior [Eq. (48)] from a denser central core to a less dense peripheral layer in the neighborhood of the shell internal wall; the average-density boundary location $\bar{\mbox{{\tengr h}}}$ [Eq. (50)] is clearly detectable and, expectedly, it decreases with increasing $N$. The calculation with $N=10^{-5}\ll 1$ served to crosscheck whether or not the numerical algorithms would reproduce the analytical solution [Eqs. (52) and (53)]; they clearly do. The density inflection point and gravitational-field minimum appear at $N^{\prime}\simeq 1.5558$; the gravitational-field inflection point appears at $N^{\prime\prime}\simeq 2.2829$. Their locations shift leftward and the minimum becomes deeper as the gravitational number increases.

The rightward continuation of Fig. 3 is presented in Fig. 4 which shows the analytical profile of the gravitational field inside the shell and in empty space up to two sphere radii produced by Eq. (45) (shell, empty space) with $\mbox{{\tengr d}}/a=0.1$ and $m_{s}/m_{g}=0.21$; the latter values were selected for the convenience of keeping the gravitational field inside the shell practically constant. The profile in Fig. 4 does not depend on the gravitational number and is unique on account of the spherical symmetry of the physical system, regardless of how matter is radially distributed inside the sphere.

Similarly to what already found out in the literature, particularly Darwin [50], Lynden-Bell and Wood [59] and Padmanabhan [61], with regard to the Lane-Emden solution, although detected and described within different perspectives, we also discovered two peculiar aspects of the physical problem: the existence of gravitational-number upper bound and of multiple solutions. They are clearly illustrated in the diagram of the density at the shell internal wall, peripheral density hereinafter, as function of the gravitational number shown in Fig. 5.

The profile practically follows a linear trend

[TABLE]

for $N\leq 1$, then it bends downward until it reaches vertical slope and marks the presence of an upper bound $N_{m}$ in correspondence to which the peripheral density reduces to 1/3 [50]; after that the profile spirals around the point $\,N=2,\mbox{{\tengr x}}(1,2)=1/3\,$ as a consequence of the existence of multiple solutions when $N>1.84$. At $N\simeq 2.20$, the profile goes through the absolute minimum $\mbox{{\tengr x}}(1,2.20)\simeq 0.264$ which represents a 74% drop below the average density; there is no solution of Eqs. (42) less dense than this one at the shell internal wall. Another interesting diagram, shown in Fig. 6, is the one that illustrates the density at the sphere center, central density hereinafter, as function of the gravitational number. The profile reveals no upper bound and goes through an infinite series of oscillations with decreasing amplitude about the vertical line $N=2$ as a consequence of the existence of multiple solutions. The value $N=2$ is rather peculiar: it presupposes the existence of infinite solutions, the most extreme of which has infinite central density.

However, we do not share the commonly accepted opinion that this asymptotic solution coincides with the singular isothermal sphere [64]

[TABLE]

This particular solution of Eq. (42a) matches the boundary condition at the shell internal wall [Eq. (42c)] only if $N=2$, in that case even verifies smoothly the normalization condition [Eq. (43)] and attains the correct value $\mbox{{\tengr x}}^{\ast}(1,2)=1/3$ of the peripheral density, but does not match at all the central boundary condition [Eq. (42b)]

[TABLE]

a detail correctly emphasized also by Saslaw [60]. So, the singular isothermal sphere described by Eq. (62) and the asymptotic isothermal sphere to which the profile of Fig. 6 tends to for $N=2$ have central gravitational fields that are, respectively, either infinitely great (and unphysical) or vanishing and there will always be a neighborhood of $\mbox{{\tengr h}}=0$, however small, in which they will differ. The density contrast, widely preferred and utilized in the literature, can be obtained as

[TABLE]

and presents a profile similar to the central density’s (Fig. 6).

No solution of Eqs. (42) exists beyond the upper bound $N_{m}$ and, therefore, no fluid-static field can exist in the gas if $N>N_{m}$; with numerical accuracy preset to eight significant digits, we obtained $N_{m}=2.51755148$ (rounded off to 2.5175 in the figures) as the greatest value for which it was possible to obtain a solution of Eqs. (42). Characteristic values corresponding to $N_{m}$ can be read off from Figs. 2–3 and Figs. 5–6. The central density rises to \,\mbox{{\tengr x}}(0,N_{m})\simeq{\color[rgb]{0,0,0}10.8}\, and the peripheral one lowers to $\,\mbox{{\tengr x}}(1,N_{m})=1/3\,$, so that the density contrast turns out to be \,\mbox{{\tengr x}}(0,N_{m})/\mbox{{\tengr x}}(1,N_{m})\simeq{\color[rgb]{0,0,0}32.4}. The average density [Eq. (50)] is attained at $\bar{\mbox{{\tengr h}}}\simeq 0.640$. The gravitational-field minimum is located at $\mbox{{\tengr h}}\simeq 0.335$ and its intensity is about 85% stronger than the value attained at the shell internal wall. From the gravitational-number definition [Eq. (37a)], we can extract critical values of radius

[TABLE]

temperature

[TABLE]

and mass

[TABLE]

and assert that no fluid-static field can exist in the gas sphere if $a<a_{c}$ or $T<T_{c}$ or $m_{g}>m_{g,c}$. To the best of our understanding, the existence of the maximum value $N_{m}$ is somewhat surprising because there does not seem to appear any sign either in Eqs. (42) or in the profiles of Figs. 2 and 3 that explicitly hints at the occurrence of limitations constraining the gravitational number.

The, at least mathematical, existence of multiple solutions for a given gravitational number is also perplexing. As representative example, we show the two solutions relative to $N=2.4$ in Figs. 7 and 8; the curves corresponding to $N=N_{m}$ (dashed lines) are included as reference. The b-solution is more extreme than the a-solution. Central cores and peripheral layers have more or less commensurate extension ($\bar{\mbox{{\tengr h}}}_{a}\simeq 0.68$ and $\bar{\mbox{{\tengr h}}}_{b}\simeq 0.59$) but the density of the b-solution is almost six times greater at the sphere center and about 0.66 times smaller at the shell internal wall with respect to that of the a-solution. The position of the gravitational-field minimum shifts from about 48% to 20% of the radius and the intensity becomes almost 2.5 times stronger.

Inevitable questions arise when looking at the results of Figs. 5–6: (a) What is the physical reason behind the existence of a maximum value of the gravitational number? Regarding this question, we definitely encourage the reader to familiarize with the interesting reflections provided and conclusion drawn by Darwin at page 18 of his communication [50] to the Royal Society of London almost 130 years ago. (b) The physical variables involved in the definition [Eq. (37a)] of the gravitational number are experimentally controllable and can be selected to produce a value greater than $N_{m}$; what happens in that case? (c) Are all multiple solutions corresponding to a given gravitational number physical? If yes, (c1) what circumstance selects the occurrence of one rather than another? If not, (c2) what reason removes the occurrence of the unphysical solutions? One thing is clear at this point: the answer to question (b) cannot be provided by a gravitofluid-static analysis; the latter simply, and only, says that solutions do not exist if $N>N_{m}$. The answer must, therefore, be sought within a gravitofluid-dynamics context. The quest for answers to the other questions brought us to investigate in details the thermodynamics of the physical system (Sec. 3).

2.5 Reflections on the Lane-Emden solution

Before turning to thermodynamics, we wish to conclude Sec. 2 by expressing our point of view, with collegial spirit absolutely free of premeditated criticism, regarding the approach more or less systematically followed in the astrophysical literature that brings to the Lane-Emden solution, particularly for isothermal gas spheres. A rather accurate and extremely informative historical account is given by Chandrasekhar in the bibliographical notes at page 176 of his textbook [58].

The conceptual pathways of the literature approach and of our approach basically separate just after Eq. (28a). The former approach carries on with the idea that the imposition of the central density

[TABLE]

is a legitimate boundary condition; Eq. (68) by itself does not remove the idea of a possible finite extension of the gas sphere but certainly pushes that idea automatically out of the mathematical description’s focus. The replacement of a gravitational boundary condition [Eq. (28b)] with a fluid-dynamic boundary condition [Eq. (68)] does not affect at all the procedure to cast the problem in nondimensional form. We introduce again nondimensional variables [Eqs. (33)] and reach once more Eqs. (34a) and (34b); Eq. (34c), however, must be replaced with the nondimensional version

[TABLE]

of Eq. (68) and, therefore, the characteristic-number set changes from Eq. (35) to

[TABLE]

There are only two independent characteristic numbers this time. We resolve Eq. (70) for the scale factors

[TABLE]

and set the P-type characteristic numbers to the values

[TABLE]

that lead to the scale factors widely used in the literature

[TABLE]

In Eq. (73), we have affixed the subscript “A” (Author) to the scale factors in order to avoid confusion with those [Eq. (40)] we use in our approach. A bit of attention should be put at remembering that there is a factor 3 between the radial-distance scale factors of, respectively, Emden (A$\rightarrow$E) and King (A$\rightarrow$K)

[TABLE]

that reflects also on the nondimensional radial coordinates

[TABLE]

The central density is taken as scale factor

[TABLE]

in both cases so that the nondimensional-density variables coincide

[TABLE]

With these premises, the mathematical problem reads

[TABLE]

An attractive feature of Eqs. (78) is the absence in them of any physical characteristic number; their solution, therefore, should be expected in the form of a universal density profile somehow apt to describe all isothermal gas spheres. On the other hand, the scale factors [Eq. (73)] contain the central density, a neither known a priori nor controllable ​​ variable ​​ in demand of being fixed by additional physical information. Until the latter is not procured somehow, we believe the density profile produced by the numerical integration of Eqs. (78) remains idle; Eq. (68) or Eq. (78c) cannot be fully flagged as boundary conditions but are, at most, scaling conditions [61] and, in turn, the ​​ density profile ​​ should be considered of parametric nature mathematically speaking. As a matter of fact, before setting in motion numerical machineries, we should dedicate a bit of attention to a question hardly asked: from a physical point of view, at which radial distance do we terminate the numerical integration? Well, if we wish to deal with our test case ($r\leq a$) via Eqs. (78) rather than Eqs. (42) then we should stop the numerical integration at the terminal distance

[TABLE]

obviously. Operationally, we need an explicit number for \mbox{{\tengr h}}_{\scalebox{0.5}{A,a}} but we see right away from Eq. (79) that we do not go very far with the expression on its rightmost-hand side because, although radius $a$ and temperature $T$ are under conceptual control, the central density ​​ escapes our reach. This reflection leads inevitably to think that the formulation based on Eqs. (78) is not self-contained. ​​

The central-density hindrance can be circumvented by calculating the mass of gas contained within the terminal distance

[TABLE]

and by making resurface from it the gravitational number

[TABLE]

Equation (81) establishes a precise connection between terminal distance and gravitational number. Moreover, the connection can even be improved because the integral can be transformed with the help of Eq. (78a) and a bit of integration-by-part jugglery as

[TABLE]

so that Eq. (81) goes into the more convenient form

[TABLE]

that spares us the trouble of the numerical evaluation of the integral. In Eq. (83), we should regard the terminal distance as unknown and the gravitational number as prescribed. It comes with no surprise that the latter also fixes nondimensionally the central density

[TABLE]

Equation (84) easily follows from the average-density definition [Eq. (14a)] taking into account Eq. (80), Eq. (82) and Eq. (83). In this way of looking at things, ​​ the gravitational number, an ingredient foreign to the mathematical formulation based on Eqs. (78), sanitizes the lack of self-containedness of that formulation and emerges as the identifier that singles out a specific isothermal gas sphere from the multitude described by the universal solution of Eqs. (78). The path is hence clear. Numerical algorithms can be smoothly launched to solve Eqs. (78) by going mathematically as far as wished from the center; for example, Emden went up to $\mbox{{\tengr h}}_{\scalebox{0.5}{E}}=2\cdot 10^{3}$, we reached $\mbox{{\tengr h}}_{\scalebox{0.5}{E}}=3\mbox{{\tengr h}}_{\scalebox{0.5}{K}}=10^{5}$ (F & FP algorithms only). This action leads to the ​​ density profile illustrated in Fig. 9 showing our numerical results with superposed data collected in the literature and adequately post-processed.

With the scale factors of Eq. (73), the nondimensional gravitational field follows straightforwardly from Eq. (24)

[TABLE]

and its profile, illustrated in Fig. 10, is thus readily obtained.

This is also a universal profile encompassing all isothermal gas spheres. The compliance of the peripheral gravitational field

[TABLE]

with the Gauss theorem’s prescript [Eq. (16)] is readily verified with the help of Eqs. (79) and (83).

With the solution of Eqs. (78) in hand, we can ​​ proceed to build the graph of Eq. (83), shown in Fig. 11, which turns out to be rather revealing because it confirms the existence of gravitational-number upper bound $N_{m}$ and of multiple solutions above $N\simeq 1.84$. The knowledge of the gravitational number allows to enter the graph on the vertical axis and to read off on the horizontal axis the terminal distance that marks the terminus on the profiles in Figs. 9 and 10 corresponding to the isothermal gas sphere identified by the prescribed gravitational number. It is this step that basically fixes the central density [Eq. (84)] and raises the standing of Eqs. (68) and (78c) to the level of legitimate boundary conditions. In order to consolidate these thoughts with an explicit example, we have considered again in Fig. 11 the case with $N=2.4$. The corresponding horizontal line intersects the profile in the points a and b and selects two values of the terminal distance: \mbox{{\tengr h}}_{\scalebox{0.5}{E,a}}=3\mbox{{\tengr h}}_{\scalebox{0.5}{K,a}}\simeq 6.09 and \mbox{{\tengr h}}_{\scalebox{0.5}{E,a}}=3\mbox{{\tengr h}}_{\scalebox{0.5}{K,a}}\simeq 14.7. In turn, these values determine the density-profile portions on the Lane-Emden solution, displayed with thicker lines in Fig. 12, that correspond to the curves a and b of Fig. 7, respectively; for the purpose of verification, we have superposed also the results extracted from Fig. 7 and adequately converted. The situation for the gravitational field is similar and illustrated in Fig. 13. For completeness, we have reproduced in Fig. 14 the diagram of the central density as function of the gravitational number from Eq. (84) for visual comparison with the one shown in Fig. 6: the coincidence is unequivocal.

The graphs in Figs. 12 and 13 may induce in the reader’s mind an impression of equivalence between the Lane-Emden solution founded on the idea of the central-density boundary condition [Eq. (78c)] and our solution founded on the peripheral-gravitational-field boundary condition [Eq. (42c)] because, after all, they lead to same results although with different parameterization; well, maybe yes from a mathematical point but certainly not from a physical point of view. In this regard, we wish to stress that the gravitational number is invisible in the formulation based on Eqs. (78) and the mathematical equivalence is established only after the physical injection of Eq. (81). We are aware that our viewpoint regarding the isothermal Lane-Emden problem and its solution is at variance with the interpretation prevailing in the astrophysical literature for which we quote as representative the statement

A pure isothermal sphere stretches to infinity and has an infinite mass

found at page 323 in Saslaw’s textbook [60], a seemingly straightforward conclusion when looking at density profile of Fig. 9 without the support of Eq. (81). Nevertheless, we do not know how to reconcile this quote, and the interpretation it represents, with the graph in Fig. 11 which clearly indicates that there is only one isothermal gas sphere that conforms with the quote: the asymptotic one identified by the rightmost intersection of the horizontal line at $N=2$ with the curve, whose intercept corresponds to an infinite terminal distance

[TABLE]

As a matter of fact, Eq. (87) implies two possibilities. One is the isothermal gas sphere we are looking for that “stretches to infinity” ($a\rightarrow\infty$) and “has an infinite mass”

[TABLE]

Its average gas density vanishes

[TABLE]

but its central density turns into an indeterminate form

[TABLE]

The other one that springs up from Eq. (87) is the isothermal gas sphere characterized by infinite central density

[TABLE]

with possibly finite but definitely indeterminate size

[TABLE]

mass and average density. We see at this point how the scale-factor choice [Eq. (73)] backfires via the bad repercussions of these peculiarities [Eqs. (90) and (91)]: Eq. (73) becomes idle for the infinitely-stretched isothermal sphere because the central density is an indeterminate form [Eq. (90)] and it goes into the perverse form

[TABLE]

for the other isothermal sphere because the central density becomes infinite [Eq. (91)]. It is rather easy to convince ourselves that Eq. (93) molds the ​​ density profile of Fig. 9 into one single physical point, i.e. the center $r=0$, in which the isothermal sphere is concentrated possibly with finite mass but definitely with infinite density.

We conjecture that the root of the interpretation prevailing in the literature can be traced back to Lane’s paper [45] whose attention concentrated on adiabatic ideal-gas spheres. In the text at page 60, Lane introduced the idea that the differential-equation solution, obviously based on Eq. (68), should also provide a way to determine gas-sphere size and central density. He began by saying that the adiabatic index, i.e. the specific-heat ratio, influences

… first the form of the curve that expresses the value of $\frac{\mbox{{\tengr r}}}{\mbox{{\tengr r}}_{\scalebox{0.6}{0}}}$ for each value of $x$; secondly, the value of the upper limit of $x$ corresponding to $\frac{\mbox{{\tengr r}}}{\mbox{{\tengr r}}_{\scalebox{0.6}{0}}}=0$; and thirdly, the corresponding value of m.

In Lane’s notation, $x$ is the nondimensional radial coordinate [defined in his Eq. (4)] and m is the gas-sphere nondimensional mass [defined in his Eqs. (5) and (6)]. Then Lane continued

These limiting, or terminal, values of $x$ and m cannot be found except by calculating the curve … But when these values have been found … they may be introduced once for all into equations (4) and (5) from which the values of $\mbox{{\tengr r}}_{\scalebox{0.6}{0}}$ and … are at once deduced.

Emden [52] extended the analysis to polytropic gas spheres and, according to the results of his calculations, summarized their behavior at page 155 as

Alle polytropen Gaskugeln von $n=0$ bis $n=\infty$ haben eine Oberfläche, an der Druck, Temperatur und Dichte die Werte Null annehmen. Für $n>5$ liegt dieselbe im Unendlichen.

[All polytropic gas spheres from $n=0$ to $n=\infty$ possess a surface on which pressure, temperature and density vanish. For $n>5$ the surface lies at infinity]

In Emden’s notation, $n$ is the polytropic index, not to be confused with the formal exponent we used in Eq. (72) (top line). There is a very short and straightforward step from the last sentence in the Emden’s statement quoted here to the conclusion for isothermal spheres ($n=\infty$) exemplified in Saslaw’s statement quoted just before Eq. (87).

Lane and Emden had no hesitation to accept the part of solution to their differential equations at the left of the “limiting, or terminal, value of $x$” at which $\mbox{{\tengr r}}/\mbox{{\tengr r}}_{\scalebox{0.6}{0}}=0$ and to dismiss the part of solution yielding negative densities at the right of that value. It seems to us that Lane’s and Emden’s choice, although justified by the necessity of building up a physical model, was simply motivated by convenience without further reflection. Already 116 years ago, Jeans signaled the conceptual precariousness of a vanishing density at page 3 of [51]:

Whether we suppose the thermal equilibrium of the gas to be conductive or adiabatic, we are still met by the difficulty that the gas equations break down over the outermost part of the nebula, through the density not being sufficiently great to warrant the statistical methods of the kinetic theory.

In line with Jeans’ appropriate warning, a fluid dynamicist would hesitate to accept and be suspicious about the physical consistency of solutions that contain vanishing mass density, pressure, temperature, because of the incompatibility of such occurrences with the reason of existence of the governing equations that generated them. In other words, the applicability of the continuum-medium model breaks down in the left vicinity of the “limiting, or terminal, value of $x$” and, therefore, what sense should one make of a result, i.e. $\mbox{{\tengr r}}/\mbox{{\tengr r}}_{\scalebox{0.6}{0}}=0$, yielded by a physical model in a place outside of its applicability domain? Of course, these considerations are absolutely in no way meant to imply rejection of the polytropic-sphere problem as conceived by Lane and Emden but we believe a prudent reflection about whether or not something may be in need of revision in their formulation would not be an idle exercise.

3 Thermodynamics

3.1 Entropy

The entropy of the physical system is the sum of the entropy of the gas and the entropy of the shell

[TABLE]

The determination of the shell contribution $S_{s}$ requires the integration on the shell volume of the relative specific entropy

[TABLE]

The determination of the gas contribution $S_{g}$ requires the integration over the gas-sphere volume of the perfect-gas specific entropy

[TABLE]

The specific entropy is composed by two terms

[TABLE]

The first descends from the translational degrees of freedom of the molecules

[TABLE]

and features a dependence on the specific volume $v(r)=1/\mbox{{\tengr r}}(r)$; the second term, $s_{gi}(T)$, is associated with the internal molecular structure and does not need to be expanded in explicit details because it depends only on temperature. In Eq. (98), we have set for brevity

[TABLE]

with

[TABLE]

The integral on the rightmost-hand side of Eq. (96) is easily performed and leads to the interesting separation

[TABLE]

Equation (100) clearly emphasizes the importance of the gas-density distribution for the purpose of entropy calculations. The reformulation of the gravitational number [Eq. (37a)] in terms of the state variables $T,V,m_{g}$

[TABLE]

gives full evidence of the gravitational correction’s lack of first-degree homogeneity; we see through Eqs. (94) and (100a), therefore, that total entropy and gas entropy are inevitably deprived of the same characteristic. The integral in Eq. (100) can be evaluated numerically after that the solution $\mbox{{\tengr x}}(\mbox{{\tengr h}},N)$ has been obtained, of course; however, with a bit of patience and integration-by-part mastery, it is also possible to obtain the integral in the following analytical form

[TABLE]

that clearly highlights the peripheral density’s role in establishing the gravitational correction. Its nondimensional diagram is shown in Fig. 15;

the pattern is obviously driven by that of the peripheral density (Fig. 5). Indeed, it also attains vertical slope at $N=N_{m}$; thereat, it crosses the level

[TABLE]

marking the boundary between upper and lower branches of fluid-static configurations. Subsequently, at $N\simeq 1.88$, it goes through the absolute minimum $I_{ent}(1.88)=-1.103$ and then spirals around the point $\,N=2,I_{ent}(2)=\ln 3-2\simeq-0.901\,$. From Eq. (94), taking into account Eq. (95), Eqs. (100) and Eq. (3.1), the explicit expression of the system entropy reads

[TABLE]

Incidentally, it may be worthwhile to point out that this entropy expression depends on the state variables $T,V,m_{g}$ and, therefore, it should not be looked at as a fundamental relation [80, 81, 82, 83]. The gravitational correction is always negative; therefore, the presence of the gravitational field always reduces the entropy with respect to the situation when it is absent. In case of multiple solutions, the upper-branch configurations have greater entropy than the lower-branch configurations. This circumstance perhaps gives a somewhat probabilistically privileged status to the former branch but it does not exclude at all the potential realizability of the latter branch.

Some of the lower-branch configurations present physically unusual peculiarities that can be brought to light by the study of the entropy differential. But first, a little simplification. It is clear from Eq. (3.1) that the shell mass $m_{s}$ is a legitimate state parameter and should be treated as such in general. On the other hand, our emphasis is not on the role of the container within the physical problem we are studying; therefore, just to simplify a little bit the forthcoming algebra, we decide at this point to introduce the (irrelevant) constraint

[TABLE]

Another important step required in the procedure to obtain the entropy differential is the logarithmic differentiation

[TABLE]

of the reformulated gravitational number [Eq. (101)]. With that, all the ingredients are in place to expand from Eq. (3.1) the entropy differential

[TABLE]

in terms of the state-variable differentials $dT,dV,dm_{g}$ and to recognize the corresponding partial derivatives by inspection of Eq. (3.1). Therein,

[TABLE]

is the shell specific heat; moreover, the following two new functions of the gravitational number appear

[TABLE]

The separation structure of Eq. (100a) is inherited by all partial derivatives in Eq. (3.1): there is a standard term appropriate when gravitational effects are absent followed by a gravitational correction. The peculiarities mentioned before descend from the capability of the gravitational corrections to flip the sign of the derivatives and this occurrence is of particular relevance for the derivative

[TABLE]

because it is a criterion of thermodynamic stability in the absence of the gravitational field and it may somehow be connected to thermal stability also in the presence of the gravitational field. The scenario is better brought forward visually by the diagram of the function $\mbox{{\tengr f}}_{ent}(N)$ shown in Fig. 16. The profile is monotonic for the upper-branch configurations (Fig. 16a) and goes through a repetitive pattern in between asymptotes for the lower-branch configurations. The asymptotes correspond to the points of vertical slope in Fig. 5 at which the derivative $\mbox{{\tengr x}}^{\prime}(1,N)$ and, by consequence, all the derivatives in Eq. (3.1) become infinite. This happens for the first time at $N=N_{m}$. The function $\mbox{{\tengr f}}_{ent}(N)$ is always negative and the derivative $\left({\partial S}/{\partial T}\right)_{V,m_{g}}$ is always positive for the upper-branch configurations (Fig. 16a); thermodynamically speaking, these are well behaved and all seems in order. The sign of the derivative can change for some of the lower-branch configurations according to where they are positioned with respect to the level established by the sum of the specific-heat terms [Eq. (110), top line] which depends on temperature and, ultimately, on the gravitational number if Eq. (101) is accordingly resolved to read

[TABLE]

For the mere purpose of visualization, we have considered a diatomic gas at sufficiently low temperature ($c_{vi}/R=1$) and set the shell contribution to a qualitative, although meaningful as order of magnitude, value ($c_{s}/R\simeq 1.573$; mild steel confining molecular nitrogen, for example); in the sequel, we will refer to these values as “visualization specific-heat settings”. In this case, the specific-heat level, labeled “shell (qualitative)” in Fig. 16, does not change with temperature or gravitational number and is obviously horizontal. So, all lower-branch configurations above the “shell (qualitative)” level have negative derivative $\left({\partial S}/{\partial T}\right)_{V,m_{g}}$, a peculiarity that we perceive as undeniable omen of thermal instability. The danger, therefore, exists exclusively for some of the lower-branch configurations among which the one at $N=N_{m}$ inaugurates the first series. The intersections of the several, actually infinite, curves of the function $\mbox{{\tengr f}}_{ent}(N)$ with the “shell (qualitative)” level identify a series of points F1, F2, F${}_{3},\ldots$ in correspondence to which the derivative vanishes

[TABLE]

They represent, therefore, the boundary points between series of well-behaved and peculiar configurations; clearly, the derivative $\left({\partial S}/{\partial T}\right)_{V,m_{g}}$ changes sign an infinite number of times.

These thermodynamic occurrences point out how much interesting and, at the same time, necessary thermodynamic-stability analyses of self-gravitating fluids are. The literature offers several [84, 59, 85, 86, 87, 61, 64, 69] but, in the majority of the references we have explored, the thermodynamic stability is systematically studied within a kinetic-theory framework centered around the Boltzmann’s definition of entropy. We believe it would be beneficial to investigate the thermodynamic stability also within an axiomatic-thermodynamics framework [80, 81, 82, 83] in order to find out whether or not such an alternative conceptual path brings to the same results achieved by kinetic-theory based analyses. Unfortunately, given the priorities of our research, we have no other choice for the time being than to label this task as future work and carry on.

3.2 Energy

The physical system is composed by matter, gas and shell, and gravitational field. Therefore, we consider two forms of energy whose belongingness is obvious: matter energy and gravitational-field energy. Their sum gives the total energy. There is no kinetic energy in static circumstances and matter energy is made up exclusively by the thermodynamic part distributed in space with density

[TABLE]

In Eqs. (113), $u_{gi}(T)$ is the contribution to specific thermodynamic energy of the gas due to internal molecular structure and $u_{s}(T)$ is the shell specific thermodynamic energy. The gravitational energy is distributed in space with density [20, 21]

[TABLE]

The expressions tagged (g), (s) and (es) in Eq. (114) are obtained by direct substitution of the gravitational-field solution [Eq. (46); Eq. (45) (shell) and (empty space)] and apply locally in the space inside the gas sphere, inside the shell and outside the shell in empty space, respectively. Our position is conceptually different from the one usually taken in the literature because we do not make use of the concept “potential energy of the gas”. The equivalence between the concepts “matter potential energy” and “field energy” holds only under stationary, a fortiori static, circumstances, as we learn from Feynman’s chapters 8 and 27 in [88]. The definition on the top line of Eq. (114) always holds, within Newton’s theory of gravity, regardless whether a dynamic or a stationary or a static situation is prevailing. The justification for our opting in favor of the “field energy” interpretation even under static circumstances is not as drastically severe as the devastating motivation for rejection expressed by Heaviside [10, 12]

Potential energy, when regarded merely as expressive of the work that can be done by forces depending upon configuration, does not admit of much argument. It is little more than a mathematical idea, for there is scarcely any physics in it. It explains nothing.

It simply resides in the wish to be self-consistent all the way along our study path because, as explained in the introduction (Sec. 1), our research objectives foresee also gravitofluid-dynamic analyses in which the equivalence does not hold.

The total-energy density is given by the sum

[TABLE]

and its integral extended to all space provides the total energy of the physical system

[TABLE]

[TABLE]

[TABLE]

We draw the reader’s attention to the fact that the presence of the shell, always systematically ignored in the literature, affects not only the matter energy via the contribution $u_{s}(T)$ but also the gravitational energy via the function $\Phi$.

There are two crosschecks that we can carry out on Eq. (117c) to verify if we are on the right track. If $N\rightarrow\epsilon\ll 1$ then Eq. (61) is applicable and we can substitute it in Eq. (117c) to obtain

[TABLE]

If the shell contribution also turns out to be negligible ($\Phi\ll 6/5$) then the gravitational energy reduces further to

[TABLE]

a well known textbook result, although the limitations ($N\rightarrow\epsilon\ll 1$ and $\Phi\ll 6/5$) of its applicability are hardly emphasized. The second crosscheck is slightly more elaborate. We make the total energy [Eq. (117a)] nondimensional by dividing it with the gravitational-energy order of magnitude $Gm_{g}^{2}/a$. Moreover, we conform to the “literature settings” by assuming a point-particle gas [$u_{gi}(T)=0$] and by neglecting the shell contributions to both matter energy [$u_{s}(T)\ll RT$] and gravitational energy ($m_{s}\ll m_{g}\Rightarrow\Phi\ll 1$). In this way, Eq. (117a) reduces to the simplified nondimensional form

[TABLE]

that corresponds numerically to the expressions deduced within the Lane-Emden’s approach and invariably used and commented throughout the literature. As visual proof, we show the diagram of Eq. (121) in Fig. 17a styled according to Fig. 2 at page 344 of LABEL:phc2002aa and in Fig. 17b styled according to Fig. 4.4 at page 317 of LABEL:tp1990pr. The graphs in Fig. 17 clearly illustrate that the profile spirals around the point $N=2\;\mbox{or}\;1/N=0.5,aE/Gm_{g}^{2}=-0.25$, and, at $N\simeq 2.03$ or $1/N\simeq 0.492$, it goes through the notorious minimum $aE/Gm_{g}^{2}\simeq-0.335$ discovered by Antonov [84] in 1962 and declared verge of the gravothermal catastrophe by Lynden-Bell and Wood [59]. Equation (121) has an undesirable feature: it diverges if $N\rightarrow 0$. This occurrence would persist even if we repristinate the neglected terms. The reader should rest assured that this divergence is not the hallmark of an additional catastrophe but is, simply, a warning that the order of magnitude $Gm_{g}^{2}/a$ is not an adequate, or perhaps convenient, scale factor of the total energy [Eq. (117a)] for small values of the gravitational number. From our point of view, the matter-energy order of magnitude $m_{g}RT$ is a more appropriate scale factor because the ensuing gravitational energy

[TABLE]

is better behaved under all circumstances. The diagram of Eq. (122) in nondimensional form is shown in Fig. 18. Before anything else, we notice that the gravitational field in the shell affects the profile and its characteristic points via the function $\Phi$. The profile is basically linear

[TABLE]

for $N\leq 1$, then it slopes down until $N=N_{m}$ where it achieves the gravitational-energy level

[TABLE]

The center of the spiral is located at $N=2$ and at the gravitational-energy level

[TABLE]

Equation (122) also presents an absolute minimum whose position depends on the function $\Phi$. The extrema are fixed by the vanishing of the derivative with respect to the gravitational number which leads to a condition

[TABLE]

that unfortunately is not exploitable analytically because it is not resolvable explicitly. In our calculations, $\Phi\simeq 0.441$ [Eq. (118c)] and we find numerically the absolute-minimum location at $N\simeq 2.32$ with minimal gravitational energy

[TABLE]

According to Eq. (117a), the explicit expression of the total energy is obtained by adding matter [Eq. (117b)] and gravitational-field [Eq. (122)] contributions

[TABLE]

The discovery of thermodynamic peculiarities connected with the study of the entropy differential [Eq. (3.1)] discussed in Sec. 3.1 auto-suggests analogous action with respect to the total-energy differential. The latter can be obtained straightforwardly from Eq. (3.2) by taking into account the same ingredients [Eq. (101), Eqs. (105) and (106)] used to derive the entropy differential plus the additional (irrelevant) constraint

[TABLE]

that we introduce by invoking the same justification we adduced for adopting Eq. (105). The final result reads

[TABLE]

and contains two new functions

[TABLE]

The specific heats on the top line of Eq. (3.2), given in Eqs. (108) in terms of the specific entropies, can be expressed also in terms of the corresponding specific energies

[TABLE]

The partial derivatives of the total energy with respect to the state variables are easily identified by inspection of Eq. (3.2) and once again we recognize in them the separation structure already encountered when dealing with the entropy partial derivatives [Eq. (3.1)]. We put on hold discussing them here, however, and just save the result in view of the thermodynamic considerations of the next section because, we believe, the discussion seems more appropriate within that context.

3.3 Axiomatic-thermodynamics perspective

3.3.1 Fundamental relation

The astrophysical literature systematically regards total energy as entitled to play the role of thermodynamic fundamental relation in the energetic representation: $E=E(S,V,m_{g})$; its inverse $S=S(E,V,m_{g})$ provides the fundamental relation in the equivalent entropic representation [80, 81, 82, 83]. If we conform to this idea then both fundamental relations are available in parametric form via Eq. (101), Eq. (3.1) and Eq. (3.2), the parameter being the temperature. Regrettably, the latter variable is not eliminable analytically within the set of the indicated equations in general and, therefore, the fundamental relation may not be obtainable in explicit form save for simplified circumstances such as that corresponding to the visualization specific-heat settings introduced just after Eq. (111). Nevertheless, we believe it is worthwhile to dig up the details for this simplified case because of the insight into the physics of the problem that can be gained. The integration of Eqs. (108) and Eqs. (132) is straightforward; the integration constants can be conveniently brought to the left-hand side of Eqs. (3.1) and (3.2) and incorporated as harmless reference levels of total energy and entropy. After that, the elimination of the temperature follows smoothly. The diagram illustrating the global view of energy versus entropy for given volume and gas mass is shown in Fig. 19. The profile has a typical monotonic trend and all seems in order and as expected for the upper-branch configurations. The peculiar part of the profile lies, of course, in the neighborhood of $N_{m}$ but it is invisible in the scale of Fig. 19. Before zooming in, however, it is appropriate first to investigate the partial derivative $\left({\partial E}/{\partial S}\right)_{V,m_{g}}$.

3.3.2 State equation $\left({\partial E}/{\partial S}\right)_{V,m_{g}}$

The state equations are provided by the fundamental relation’s partial derivatives and the latter can be obtained in a straightforward manner from the entropy and total-energy differentials [Eqs. (3.1) and (3.2)] without specifying any simplification regarding the specific-heat terms. The recipe is rather simple. Suppose we wish to obtain $\left({\partial E}/{\partial S}\right)_{V,m_{g}}$: then first we set $dV=dm_{g}=0$ in Eqs. (3.1) and (3.2) and after eliminate $dT$ from the reduced expressions of the differentials. This algebraic manipulation leads to a result

[TABLE]

that provokes some anxiety because, as it stands, it seems to deny the variables $(T,S)$ the status of couple of conjugated variables in the energetic representation. Actually, Eq. (3.3.2) teaches a good lesson: in the presence of the gravitational field, we should not make the mistake of taking for granted thermodynamic notions we are accustomed to in the absence of the gravitational field. Incidentally, this is the reason why we refrained to put instinctively the stamp of thermal-stability criterion on Eq. (110), as we can do in the absence of the gravitational field, and awaited, instead, to follow through with the standard analysis of the fundamental relation and its first/second derivatives. Of course, the first thought that comes to mind when looking at Eq. (3.3.2) is to compare the functions $\mbox{{\tengr f}}_{ent}(N)$ and $\mbox{{\tengr f}}_{ene}^{\mbox{\tiny+}}(N)$ [Eqs. (109a) and (131a)], although they popped in rather independently. To this aim, we have superposed the function $\mbox{{\tengr f}}_{ene}^{\mbox{\tiny+}}(N)$ on the profiles of Fig. 16. The outcome is shown in Fig. 20 and reveals an unexpected surprise: the two functions coincide

[TABLE]

The immediate consequence is the reduction of Eq. (3.3.2) to the more reassuring form

[TABLE]

that sanctions the status of $(T,S)$ as couple of conjugated variables and, at the same time, dissipates the anxiety provoked by Eq. (3.3.2). With reassurance [Eq. (135)] that the slope of the profile in Fig. 19 is the temperature, we can swiftly determine also the curvature

[TABLE]

and look back now at Eq. (110) as legitimate criterion of thermal stability also in the presence of the gravitational field. On the basis of Eq. (134), with a view to the differentials of entropy [Eq. (3.1)] and total energy [Eq. (3.2)], we can affirm that

[TABLE]

and recognize that

[TABLE]

for the thermally unstable configurations. With regard to Eq. (138), we ​​ are aware of Lynden-Bell and coauthors’ crusade [89, 59, 90] ​​ in support of existence and physical plausibility of a negative heat capacity of the gas and the controversy it spawned in the literature [91, 92, 93, 94, 95]. However, we refrain from discussing virial-theorem applications and statistical-ensemble nonequivalence because it would take us far away from the mainstream of the present discourse and postpone the exposition of our viewpoint to future communications. Having taken that commitment, here we just state that we do not interpret ​​ $\left({\partial E}/{\partial T}\right)_{V,m_{g}}$ as heat capacity, in general, and of the gas, in particular, because its more explicit expression

[TABLE]

adapted from Eq. (3.2) includes not only terms belonging to the gas [Eq. (3.3.2), top line] but also properties belonging to the shell [Eq. (3.3.2), middle line] and a correction belonging to the gravitational field [Eq. (3.3.2), bottom line]; therefore, what sense would it make to regard $\left({\partial E}/{\partial T}\right)_{V,m_{g}}$ as a physical property of the gas? ​​

From our standpoint, we view the mentioned terms on top and middle lines of Eq. (3.3.2) as, respectively, gas and shell heat capacities; they are always positive, by definition. It is the gravitational correction on the bottom line of Eq. (3.3.2) that can flip the sign of $\left({\partial E}/{\partial T}\right)_{V,m_{g}}$ and introduce thermal instability. ​​

Another amazing consequence of Eq. (134) is the possibility to obtain an analytical expression for the derivative $\mbox{{\tengr x}}^{\prime}(1,N)$, as anticipated in Sec. 3.1 near Eq. (109c). Substituting Eqs. (109a) and (131a) into Eq. (134) and solving for the derivative yields

[TABLE]

We see in Eq. (140) the analytical confirmation of a numerical result displayed in Fig. 5: all vertical-slope configurations share the peripheral-density level $\mbox{{\tengr x}}(1,N)=1/3$. As further verification, we show the comparison of numerical [Eq. (109c)] versus analytical [Eq. (140)] values of $\mbox{{\tengr x}}^{\prime}(1,N)$ in Fig. 21; the agreement is satisfactory notwithstanding a slight deterioration in accuracy of the D-algorithm (solid circle) for the lower-branch configurations in the vicinity of the asymptote $N=1.84$.

3.3.3 Fundamental relation (reprise).

We return now to the profile in Fig. 19 and resume the discussion regarding the fundamental relation put on hold in the paragraph ending Sec. 3.3.2. Having acquired reliable knowledge about slope [Eq. (135)] and curvature [Eq. (136)], we proceed to expand the axes’ scales to bring the situation in the vicinity of $N_{m}$ into focus. The zoomed-in view illustrated in Fig. 22a reveals a singular zigzag pattern whose cusps correspond to the characteristic points F1, F2, F${}_{3},\ldots$ we have encountered in Fig. 16. At the right of $N_{m}$, as we already know from Fig. 19, the configurations are thermally stable [$\left({\partial S}/{\partial T}\right)_{V,m_{g}}>0$] and the profile has positive curvature

[TABLE]

At $N_{m}$, the derivative $\left({\partial S}/{\partial T}\right)_{V,m_{g}}$ diverges to $\pm\infty$ [Eq. (110) and Fig. 16a] and the profile goes through an inflection point

[TABLE]

where the curvature transitions from positive to negative. At the left of $N_{m}$, there is the first series of thermally unstable configurations [$\left({\partial S}/{\partial T}\right)_{V,m_{g}}<0$] with negative curvature

[TABLE]

which terminates in F1 where the slope, i.e. the temperature, is continuous but the curvature flips from $-\infty$ to $+\infty$ reintroducing thermal stability. The sequence “thermally stable configurations - inflection point - thermally unstable configurations” repeats in between each couple of cusps endlessly, in consistent agreement with the picture conveyed by Fig. 16. Each thermally unstable configuration, conceded though that there is a feasible way to prepare the physical system in such desired configuration, is clearly the brink of a potential gravothermal catastrophe because any thermal disturbance, however small, will knock the physical system off and make it drift away from the initial configuration: if heated up/cooled down then the physical system will demand/yield more energy. The location of the cusp F1 is important because it signals the existence of bounding values of entropy and total energy: no fluid-static configuration exists to the left of and/or below F1. Therefore, if the physical system turns out to be in a thermodynamic state with entropy and/or total energy beyond the bounds indicated by F1 then it will necessarily be in a gravitofluid-dynamic condition. In order to emphasize more incisively this aspect, we show in Fig. 22b the zoomed-in view relative to the calculation of the fundamental relation with the literature settings because, in it, the reader will feel comfortable to recognize the Antonov’s well-publicized total-energy minimum (-0.335) and more receptive about the existence of a, perhaps less-publicized, entropy minimum (-2.094).

We ought to point out that the thermal-stability conclusions we reached along the axiomatic-thermodynamics guidelines differ to some extent from what reported in the literature. For example, Padmanabhan [61] referred to and extended Antonov’s analysis [84]; with reference to his Fig. 4.4, reproduced here in Fig. 17b, he wrote at page 317

As we shall see later, the branch CD is unstable and is not physically realisable. (Point C corresponds to a density contrast of 709; thus isothermal spheres with a density contrast in the range (32,709) have a negative specific heat.)

We share Padmanabhan’s opinion that thermally unstable configurations may not be physically realizable. On the other hand, our findings indicate that not all the configurations on his CD branch, which corresponds to the zigzagging arcs F1F2, F2F${}_{3},\ldots$ in Fig. 22b, are thermally unstable. Moreover, the quoted statement seems to exclude his BC branch, which represents the “isothermal spheres with a density contrast in the range (32,709)” and corresponds to the arc $N_{m}$F1 in Fig. 22b, from the set of thermally unstable configurations. This impression is corroborated by statement (iii) at page 323 of [61]

Systems with $RE/GM^{2}>-0.335$ and $\mbox{{\tengr r}}_{c}<709\,\mbox{{\tengr r}}_{e}(R)$ can form isothermal sphere which are local maxima of entropy.

In Padmanabhan’s notation, $R$ is the sphere radius, corresponding to our $a$, $\mbox{{\tengr r}}_{c}$ and $\mbox{{\tengr r}}_{e}(R)$ are central and peripheral densities respectively. Padmanabhan’s latter quote is aligned with the more explicit conclusion drawn by Katz at page 768 of [86]

The branch of the curve between $h=1$ (at infinity) and $h=709$ is a branch of stable configurations.

while referring to his Fig. 3 which illustrates the curve $1/T$ against $E$; in Katz’s notation, $h$ is the density contrast. However, the text in parentheses after the semicolon in Padmanabhan’s former quote explicitly acknowledges the applicability of Eq. (138), unambiguous mark of thermal instability according to axiomatic thermodynamics, to his BC branch. We perceive irreconcilable contradictions in the quoted statements.

3.3.4 State equation $\left({\partial E}/{\partial V}\right)_{S,m_{g}}$

We move on now to next partial derivative $\left({\partial E}/{\partial V}\right)_{S,m_{g}}$. The application of the recipe described in the paragraph just before Eq. (3.3.2) and adapted to the present case yields

[TABLE]

Apart the minus sign, the right-hand side of Eq. (144) shows the difference between the peripheral pressure and a term due to the shell’s presence. This state equation, with proper consideration of Eqs. (14) and Eq. (101), expands into the form

[TABLE]

bearing a strong resemblance to the state equation that Bonnor [57] displayed as Eq. (1.2) and claimed to Terletsky (Ref. 7 in Bonnor’s list of references); the gravitational effects reside in the presence of $\mbox{{\tengr x}}(1,N)$, which Terletsky did not have because he considered a gas of uniform density, and of the second term on the right-hand side whose structure indeed conforms with the correction suggested by Terletsky and whose origin is clearly attributable to the shell.

Taking into account the definition of the gravitational energy [Eq. (122)], the right-hand side of Eq. (144) can be easily manipulated in order to rephrase the same equation into the, perhaps more meaningful, form

[TABLE]

Equation (146) informs that the profile of $\mathfrak{p}/\bar{p}$ is simply a rescaling of the gravitational-energy profile (Fig. 18) and, apparently, does not add any novel information. However, a curious thermodynamicist would definitely wonder how the isotherms in the $\mathfrak{p},V$ plane look like and for a good reason: the isotherm’s slope $\left({\partial\mathfrak{p}}/{\partial V}\right)_{T,m_{g}}$ constitutes a criterion of thermodynamic stability [57, 96, 60]. We do know this, by the way, strong of the reassuring Eq. (135). Equation (144) requires a bit of adaptation in order to satisfy the thermodynamicist’s curiosity. The step sequence consists in resolving the volume in terms of the gravitational number from Eq. (101)

[TABLE]

then making the average pressure [Eqs. (14)] explicit in Eq. (144)

[TABLE]

and finally replacing the volume in Eq. (148) with Eq. (147) to obtain

[TABLE]

As visual proof of the equivalence between Bonnor’s parametric equations and our Eqs. (147) and (149b), we show the reproduction of Bonnor’s Fig. 1 styled according to the original’s CGI units in Fig. 23a and in the corresponding nondimensional form in Fig. 23b.

Returning now to the more complete case with inclusion of the shell’s term [Eq. (149a)], we show the isotherms’ nondimensional profile in Fig. 24, in logarithmic scale for clearer visualization of the characteristic points labeled A,B,C in compliance with Bonnor’s notation; for convenience, we have also included the cusp F1. Their locations depend, in general, on the function $\Phi$. The spiral center is positioned at $N=2\;\mbox{or}\;1/N^{3}=1/8$ and at the level $8/3(1-\Phi)$; if $\Phi=0$, the level reduces to $8/3$ in full agreement with Bonnor’s results. The left bound B is located at $N=N_{m}\;\mbox{or}\;1/N_{m}^{3}\simeq 0.063$ and at the level $N_{m}^{3}(1-\Phi N_{m}/2)/3$. The determination of the extrema A, C, etc., requires the vanishing of the derivative

[TABLE]

which follows smoothly from Eqs. (147) and (149a) by logarithmic differentiation; the nondimensional profile of the negative isotherm’s slope is shown in Fig. 25 to facilitate thermodynamic-stability considerations.

Unfortunately, the zeros of Eq. (150) cannot be found analytically. In our calculations, $\Phi\simeq 0.441$ [Eq. (118c)] and the numerically found values for the positions of A,C are indicated in Fig. 24; the values corresponding to Bonnor’s case ($\Phi=0$) are indicated in Fig. 23b. The profiles in Figs. 24 and 25 give unequivocal evidence of the existence of a) configurations that can simultaneously be thermally stable/unstable and mechanically unstable/stable and b) upper-branch configurations that are mechanically unstable. Indeed, the configurations to the right of A ($N_{\mbox{\tiny A}}\simeq 2.25$ or $1/N^{3}_{\mbox{\tiny A}}\simeq 0.088$) and on the arc AB belong to the upper branch and are all thermally stable (Figs. 16a and 22a). Those to the right of A are also mechanically stable because $-\!\left({\partial\mathfrak{p}}/{\partial V}\right)_{T,m_{g}}>0$ for them; instead, those on the arc AB have $\,-\!\left({\partial\mathfrak{p}}/{\partial V}\right)_{T,m_{g}}<0$ and are mechanically unstable. The stability condition swaps for the configurations on the arc BF1: they are thermally unstable but mechanically stable because $-\!\left({\partial\mathfrak{p}}/{\partial V}\right)_{T,m_{g}}>0$ again. The described pattern repeats endlessly as the profile keeps spiraling.

We are aware that Bonnor’s [57] and our axiomatic-thermodynamics flavored conclusions differ. We have done our homework to follow his reasoning, presented at page 356, but have not been able to reconcile the differences. First of all, we have verified analytically that his Eq. (2.16) is contained in our Eq. (150); with a bit of attention to notation conversion, the verification is rather easily accomplished by substituting Eq. (140) in Eq. (150) and then setting $\Phi=0$. After that, we could consider with reassurance the graph of Bonnor’s Eq. (2.16), unfortunately not shown in [57], adequately represented by the profile of Fig. 25; the locations of A and C are slightly displaced if $\Phi=0$ but the aspect of the profile in terms of slope and curvature remains similar. Subsequently, we have proceeded to further derivation of Eq. (150) to determine the second derivative

[TABLE]

because it is a key element in Bonnor’s argumentation. The partial derivative on the bottom line of Eq. (151) is the slope of the profile in Fig. 25; it does not need to be expanded for the purpose of the present discussion. Equation (151) informs that, in A and C, such a slope and the second derivative $\left({\partial^{2}\mathfrak{p}}/{\partial V^{2}}\right)_{T,m_{g}}$ are directly proportional because the first derivative $\left({\partial\mathfrak{p}}/{\partial V}\right)_{T,m_{g}}$ vanishes. With these preliminaries in hand, we concentrated on Bonnor’s text. At page 356, just below Eq. (3.2), he considered the situation at A and recognized the negativity of the second derivative

Let $V=v$ be the smallest volume for which (3.1) (and therefore $\partial p/\partial V$) becomes zero; then one easily finds

[TABLE]

although without specifying how to find it analytically. Geometrically, one observes the curvature at A in the profiles of Fig. 23a or Fig. 24; analytically, Eq. (151) settles the second-derivative negative sign because the slope in the profile of Fig. 25 at A is positive. Then Bonnor considered the consequence of a small fluctuation and drew the correct conclusion

Thus for a small fluctuation $dV$ in the volume of the sphere of mean volume $v$, we have

[TABLE]

Therefore if a fluctuation occurs which result in a slight decrease in the volume of the sphere, this will produce a decrease in the pressure inside its boundary which will lead to a further reduction in its volume.

Now, the reasons that have lead Bonnor to draw from the quoted statements both the specific conclusion

Thus an isothermal gas sphere of volume greater than $v$ will be unstable because small fluctuations will make the volume $v$ collapse towards the centre.

and the general conclusion at page 357

…we can now see that any point on the spiral part of the curve below A refers to a state of unstable equilibrium.

meticulously supported by the remark

It is necessary to point this out because along BC, for example, $\partial p/\partial V$ is positive, so it might seem that this part of the curve refers to stable states.

remain obscure to us. As a matter of fact, the repetition of Bonnor’s reasoning for the configuration C brings to the inescapable conclusion of mechanical stability. The volume at C is greater than the volume at A (Bonnor’s $v$). The slope of the profile in Fig. 25 is negative at C and Eq. (151) yields a positive second derivative; therefore, a small volume-decreasing fluctuation around C will produce an increase in the peripheral pressure so that the gas sphere will not ”collapse towards the center” but will return to C.

3.3.5 State equation $\left({\partial E}/{\partial m_{g}}\right)_{S,V}$

For the sake of completeness, we recall that there is a third partial derivative $\left({\partial E}/{\partial m_{g}}\right)_{S,V}$ hardly considered or even mentioned in the literature. We adapt one more time the recipe described in the paragraph just before Eq. (3.3.2) and obtain

[TABLE]

The state equation defined by Eq. (3.3.5) is important because its derivatives $\left({\partial\mathfrak{g}}/{\partial V}\right)_{T,m_{g}}$, $\left({\partial\mathfrak{g}}/{\partial m_{g}}\right)_{T,V}$ together with $\left({\partial\mathfrak{p}}/{\partial V}\right)_{T,m_{g}}$ should be involved in a third criterion of thermodynamic stability connected to the physical-system response to mass variations. However, pending our task (last paragraph of Sec. 3.1) of carrying out a thermodynamic-stability analysis within an axiomatic-thermodynamics framework, here we simply limit ourselves to mention the existence of the third stability criterion and the necessity to explore it but prefer to refrain from engaging prematurely and intuitively into details, necessarily extrapolated from our experience with thermodynamic systems without gravitational fields, that could turn out to be incomplete and/or inaccurate.

3.3.6 Fundamental-relation inhomogeneity

The state-equation definitions [Eq. (135), Eq. (144) and Eq. (3.3.5)] legitimate the Gibbs equation both in the energetic representation

[TABLE]

and in the entropic representation

[TABLE]

The Euler equation, however, does not hold because entropy and total energy do not possess first-order homogeneity; indeed, it requires a simple and straightforward algebra to obtain

[TABLE]

The term on the bottom line of Eq. (3.3.6) measures the inhomogeneity of entropy and energy; its nondimensional profile is shown in Fig. 26.

4 Conclusions and future work

The main objective of this paper, described in the paragraph beginning after Fig. 1, has been achieved: the gravitofluid-static fields have been determined and they are at our disposal as initial fields to start gravitofluid-dynamic calculations. We have followed a solution pathway different [Eq. (28b) and Eq. (42c)] from the widely adopted one on which the Lane-Emden solution is based but most of our results are aligned with those of the latter approach. We have learned that the gravitofluid-static problem [Eqs. (42)] hinges on one single characteristic number [Eq. (37a)] that we call the gravitational number because it measures the importance of gravity effects. With respect to uniform gas conditions and gravitational-field linear profile prevailing when the gravitational number is vanishingly small, density/pressure and gravitational-field profiles (Figs. 2 and 3) vary smoothly and monotonically with increasing gravitational number until, somewhat surprisingly, the latter reaches an upper bound $(N_{m}\simeq 2.51755148)$ beyond which no gravitofluid-static configuration exists. For a prescribed gravitational number belonging to the interval $[1.84,N_{m}]$, we have also found out (Fig. 5) the existence of multiple solutions. These findings are in accordance with what already known in the literature and generated in our mind questions, formulated in the last paragraph of Sec. 2.4, whose answers’ quest lead us to investigate the thermodynamics of the physical system, specifically entropy (Sec. 3.1) and total energy (Sec. 3.2) with the aim in mind to derive the all-governing fundamental relation according to the guidance of axiomatics thermodynamics [80, 81, 82, 83] (Sec. 3.3). In general, thermodynamic properties comply with a separation structure composed by a standard term when gravitational effects are absent plus a gravitational correction [Eqs. (100) for example] that removes first-order homogeneity because of the long-range nature of gravity forces. This was our first contact with non-extensive thermodynamics and we really experienced in full the excitement noted in a comment of Prof. Landsberg [97, end of first paragraph on page 49]

Here we have a thermodynamic system which is different and therefore of exciting novelty.

With the fundamental relation in hand, we obtained both the first derivatives [Eq. (135), Eq. (144), Eq. (3.3.5)], i.e. the state equations, and the second derivatives [Eq. (136), Eq. (150)] that, by conceptual extrapolation from our past familiarity with thermodynamics without gravitational fields, we believe should have a direct bearing on the thermodynamic stability of the physical system. The gravitational correction to the gas entropy (Fig. 15) permits unambiguously to distinguish upper and lower branches of gravitofluid-static configurations. It turns out that all upper-branch configurations are thermally stable (Fig. 16a, Fig. 19). Regarding the lower-branch configurations (Fig. 16, Fig. 22), some of them are thermally unstable and, therefore, not physically realizable; the first series begins with the configuration relative to the upper-bound value $N_{m}$ of the gravitational number; the latter’s origin, therefore, appears to be of thermodynamic nature. We could take the last statement as an, admittedly vague, answer to question (a) (last paragraph of Sec. 2.4) although we also have to surrender humbly to the same defeat acknowledged by Darwin at page 19 of [50]

… I am unable to find any analytical relationships by which the minimum value of $\frac{1}{3}\mbox{{\tengr b}}^{2}$ can be deduced.

130 years ago; in Darwin’s notation, $\frac{1}{3}\mbox{{\tengr b}}^{2}$ corresponds to our $1/N$. Other lower-branch configurations are as legitimately thermally stable as those of the upper branch and we could not find any physical argument within the gravitofluid-statics context preventing their realizability. We entrust our planned gravitofluid-dynamics study with the hope to shed light on the multiple-solution conundrum. In the thermodynamic plane $E,S$ (Fig. 22), series of thermally stable and unstable configurations alternate and are separated by cusps; the cusp F1 is important because it flags limiting values of entropy and total energy to the left of or below which the physical system must necessarily be in a gravitofluid-dynamic condition. Moving on to mechanical stability (Fig. 24, Fig. 25), the striking feature is the existence of upper-branch configurations (arc AB in Fig. 24) that are mechanically unstable, their thermal stability notwithstanding. For the lower-branch configurations, we find on the isotherm’s profile the same alternating pattern of mechanically stable and unstable series as we have seen to exist on the fundamental-relation profile regarding thermal stability. We have mentioned the existence of a third criterion of thermodynamic stability [paragraph after Eq. (3.3.5)] connected to mass variations that should be considered on equal footing with the other two criteria we have discussed. We did not feel comfortable, though, to plunge into third-criterion details and proceed with intuitive elaborations; we preferred to hold on with this matter in view of our planned task to carry out a thorough thermodynamic-stability study along the axiomatic-thermodynamics guidelines. The latter constitutes one step of our future work. The other steps will be gravitofluid-dynamics computational studies of our test case. The first one will concentrate on a spherical-symmetric motion of the gas; the targets are clearly the answers to questions (b) and (c) and, additionally, the verification of the consequences predicted by the thermodynamic-stability criteria. The second one will focus on an axisymmetric motion of the gas compatible with both Netwon’s and gravitomagnetic theories of gravity in the hope to bring to light, at least qualitatively, detectable differences in the respective flow fields.