format: text

text:
# Spherically symmetric relativistic stellar structures

## 1 Introduction

The line element of a static spherically symmetric spacetime can be written as

$$ds^2=c^2e^{2\mathrm{\Phi }(r)}dt^2+e^{2\mathrm{\Lambda }(r)}dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),$$
(1)
where $`r`$ is the radial area coordinate and $`c`$ the speed of light, which we have chosen to introduce explicitly. A perfect fluid is described by the energy-momentum tensor $`T_{\mu \nu }=\rho u_\mu u_\nu +p(g_{\mu \nu }+u_\mu u_\nu )`$, where $`\rho `$ is the mass-energy density and $`p`$ the pressure as measured in a local rest frame; compatibility with static symmetry requires the unit vector $`u^\mu `$ to be parallel to the static Killing vector, i.e., $`u^\mu =c^1e^\mathrm{\Phi }\delta _0^\mu `$. Einstein’s field equations can be written so that the “relativistic potential” $`\mathrm{\Phi }`$ is given by

$$\frac{d\mathrm{\Phi }}{dr}=\frac{1}{p+\rho c^2}\frac{dp}{dr},$$
(2)
while the remaining equations take the form

$`{\displaystyle \frac{dm}{dr}}`$ $`=4\pi r^2\rho ,`$ (3a)
$`{\displaystyle \frac{dp}{dr}}`$ $`={\displaystyle \frac{Gm\rho }{r^2}}\left(1+{\displaystyle \frac{p}{\rho c^2}}\right)\left(1+{\displaystyle \frac{4\pi r^3p}{mc^2}}\right)\left(1{\displaystyle \frac{2Gm}{c^2r}}\right)^1,`$ (3b)
where $`G`$ is the gravitational constant and where the mass function $`m(r)`$ is defined by $`\mathrm{exp}(2\mathrm{\Lambda })=(1\frac{2Gm(r)}{c^2r})^1`$. To form a determined system, these equations must be supplemented by equations characterizing the properties of the perfect fluid matter, e.g., by a barotropic equation of state $`\rho (p)`$.

Eqs. (3) are the equations that govern hydrostatic equilibrium of a self-gravitating spherically symmetric perfect fluid body, Eq. (3b) is the famous Tolman-Oppenheimer-Volkoff (TOV) equation . The TOV equation is fundamental to the study of relativistic stellar models, nonetheless it is associated with certain problems. One is that the equation is singular at $`r=0`$, which leads to mathematical complications. It was only recently that such central issues as existence and uniqueness of regular static spherically symmetric perfect fluid models could be established . A less known but perhaps even more problematic feature is that the right hand sides of Eqs. (3) are not differentiable in $`p`$ at $`p=0`$ for many barotropic equations of state, e.g., for asymptotically polytropic equations of state, discussed below. A similar problem occurs when $`p\mathrm{}`$; even though this limit is unphysical, this is still important. The relevance of the limit is due to the fact that there are several physical problems that are directly or indirectly governed by the limit $`p\mathrm{}`$. Consider, for example, the Buchdahl inequality $`2GM/c^2R8/9`$, where $`M`$ is the total mass and $`R`$ the surface radius, which implies that no star can have a larger gravitational redshift than $`2`$ at the surface. All relativistic stellar models satisfy the Buchdahl inequality, however, the equality $`2GM/c^2R=8/9`$, is associated with a solution of constant mass-energy density and infinite central pressure. Among other issues relying on differentiability of $`\rho (p)`$ in the limit $`p\mathrm{}`$ there is, e.g., the mass-radius relationships for large central pressures which we will discuss in the present paper. The derivation of the Buchdahl inequality exemplifies that sometimes certain rather unphysical solutions, associated with somewhat unphysical equations of state, are important for physical problems. Hence it is useful, indeed essential for certain problems, to understand the entire solution space associated with Eqs. (3), including especially the limiting solutions with infinitely small or infinitely large pressures. In particular, in order to study generic behavior of solutions, one has to understand the entire solution spaces of large classes of equations of state. As indicated, the TOV equation is rather unsuitable for this task.

The purpose of this paper is to introduce a formulation that avoids the above problems and allows us to (i) obtain a global picture of the entire solution space of a large set of equations of state and to (ii) probe relationships between the equation of state and global features like mass-radius properties. To this end we will develop a framework based on the theory of dynamical systems; Eqs. (3) will be reformulated as a regular autonomous system of differential equations. A formulation that allows such a global treatment of solutions cannot exist for arbitrary barotropic equations of state. However, the large classes of equation of state that can be covered are of fundamental importance, as outlined in the following.

Let us start by making some dimensional comparisons between General Relativity (GR) and Newtonian gravity. The dimensions relevant for self-gravitating perfect fluid models are mass, length, and time. In Newtonian theory and in GR these are linked through the gravitational constant $`G`$. Thus one can choose units so that one unit, e.g., mass, is expressed in terms of the other two. In relativity yet another fundamental constant appears – the speed of light $`c`$. This constant reduces the number of independent units to one, e.g., length. This restriction on the number of independent scales is of great significance since scaling laws are of great importance in gravitational theory, as in all branches of physics (see, e.g.,  for references).

The Newtonian limit of (2) and (3) is obtained by defining the Newtonian potential111If one uses geometrized units $`c=1`$, $`G=1`$, there is no need for a distinction between $`\mathrm{\Phi }`$ and $`\mathrm{\Phi }_N`$. $`\mathrm{\Phi }_N:=c^2\mathrm{\Phi }`$ and letting $`c^2\mathrm{}`$, thus setting the relativistic ”correction” terms to zero. Note that $`\rho `$ stands for the rest-mass density in the Newtonian case. Requiring that these equations be invariant under scalings of time and space forces a barotropic equation of state to be a polytrope $`\rho =Kp^{n/(n+1)}`$, where $`K>0`$ and $`n0`$ is the so-called polytropic index. The appearance of $`c`$ in relativity reduce the number of scale degrees of freedom from two to one since extra terms arise in the equations as compared to Newtonian theory; these terms originate from spacetime curvature and from pressure, which now also contributes source terms, and are proportional to $`c^2`$. Hence polytropes do not leave the relativistic equations invariant under scale transformations. Only the linear and homogeneous equation of state $`p=K\rho `$ does.

Invariance implies symmetry, and symmetry implies mathematical simplification. The existence of scale symmetries allows one to choose all variables except one to be scale-invariant. Doing so leads to that the equation for the single non-scale-invariant variable decouples and thus one obtains a reduced set of equations for the scale-invariant variables. Hence polytropes and linear homogeneous equations of state play a special role in Newtonian theory and relativity, respectively — they are the mathematically simplest perfect fluid models these theories admit, and the associated equations can be formulated as symmetry-reduced problems.222Other equations of state may lead to greater simplifications for problems with extra symmetries, but it is only the polytropic and the linear equations of state that admit symmetry reductions in the general case in Newtonian gravity and relativity, respectively.

The central role these equations of state occupy is not due to gravity alone, but is also based on microscopic considerations of matter models. Again scaling laws can be used to motivate their special status, from a fundamental as well as phenomenological perspective. Consider, e.g., Chandrasekhar’s equation of state for white dwarfs, where degenerate electrons contribute the pressure while the baryons supply the mass. In this case there are two limits — the non-relativistic electron limit, which leads to a polytropic equation of state with index $`n=3/2`$, and the relativistic electron limit, which yields a polytrope with index $`n=3`$. Another example is given by an ideal Fermi gas. For extreme relativistic fermions $`p=\rho c^2/3`$ describes the asymptotic behavior of the equation of state, where the fermions supply both the pressure and the mass-energy.

This discussion suggests that it is natural to consider equations of state that are at least asymptotically scale-invariant, i.e., asymptotically polytropic or linear, and to classify equations of state according to their asymptotic limits; $`p0`$ and $`p\mathrm{}`$. Their special status indicates that asymptotically polytropic/linear equations of state are the mathematically “simplest” equations of state — apart from their exact counterparts.333Again we are referring to generic models, special cases, e.g., special static spherically symmetric models, might admit ”hidden” symmetries, and may thus be mathematically even simpler or lead to explicit solutions, see ,  and references therein.

In GR, asymptotically linear equations of state possess asymptotic invariance properties, however, surprisingly at first glance, also asymptotic polytropes exhibit similar mathematical simplicity in the low pressure regime.444This is satisfactory since in this regime asymptotically polytropic equations of state are physically much more interesting. The reason for this is the following: imposing spherical symmetry and requiring a static perfect fluid equilibrium configuration results in a close connection between Newtonian gravity and GR. The spacetime symmetries prevent the existence of gravitational waves and gravitomagnetic effects, thereby reducing the number of degrees of freedom in GR to the same number as in the Newtonian case. Asymptotically, the relativistic equations of hydrostatic equilibrium coincide with their Newtonian counterparts provided that $`p/\rho c^20`$ when $`p0`$, and hence it follows that variables that are asymptotically scale-invariant in the Newtonian case are also asymptotically scale-invariant in GR.555This idea was first introduced by Nilsson and Uggla , however, it will be exploited further here. Accordingly, asymptotically polytropic equations of state can be naturally included in the relativistic formalism.

The paper is structured as follows. In Sec. 2 the dynamical systems approach to relativistic stellar models is presented: by introducing asymptotically scale-invariant variables we reformulate the equations as a dynamical system. In Sec. 3 we discuss the assumptions on the equations of state that are needed for the dynamical system to become regular. Thereby the notion of asymptotically polytropic/linear is made precise. The section concludes with a discussion about phase transitions. In Sec. 4 we investigate the state space of the dynamical system and study its global properties. The dynamical systems picture is subsequently translated into a physical picture in Sec. 5. In Sec. 6, we prove theorems concerning the relationship between the equation of state and global mass and radius features. To further illustrate our approach, we give several numerical examples in Sec. 7. Finally, in Sec. 8 we give a concluding discussion, while we briefly compare the relativistic and Newtonian cases in the Appendix.

## 2 The dynamical system

Let us begin with some basic assumptions and definitions. Throughout this paper we will assume that the perfect fluid is characterized by a non-negative mass-energy density $`\rho `$ and a non-negative pressure $`p`$, related by a continuous barotropic equation of state $`\rho (p)`$ that is sufficiently smooth for $`p>0`$.666We will see in Sec. 3 that, e.g., $`𝒞^2`$ is sufficient, although this restriction can be weakened. Below we show how to handle even less restrictive situations like phase transitions. Define

$$\mathrm{\Gamma }_N(p):=\frac{\rho }{p}\frac{dp}{d\rho },\mathrm{{\rm Y}}(p):=\mathrm{\Gamma }_N^1(p)=\frac{p}{\rho }\frac{d\rho }{dp},\sigma (p):=\frac{p}{\rho c^2},$$
(4)
where $`\mathrm{\Gamma }_N`$ is the standard Newtonian adiabatic index777However, note that $`\rho `$ stands for the rest-mass density in the Newtonian case. , however, it turns out to be more convenient to use the inverse index $`\mathrm{{\rm Y}}`$, which in addition naturally incorporates incompressible fluids (for which $`\mathrm{{\rm Y}}=0`$) into the formalism. The frequently used polytropic index-function $`n(p)`$ is defined via $`\mathrm{{\rm Y}}(p)=n(p)/(n(p)+1)`$, or equivalently $`\mathrm{\Gamma }_N(p)=1+1/n(p)`$.

###### Assumptions.

For simplicity, although only necessary asymptotically in our framework, we assume that $`d\rho /dp0`$. Then $`\mathrm{{\rm Y}}0`$ and $`\rho (p)`$ is monotonic. If one impose the causality condition $`c_sc`$ on the speed of sound $`c_s=\sqrt{dp/d\rho }`$, then $`\sigma \mathrm{{\rm Y}}`$. Below, further assumptions will be imposed by the dynamical systems formulation.

A central idea of this paper is to obtain a regular dimensionless autonomous system from Eqs. (3) with a compact state space. To this end we first ”elevate” $`r`$ to a dependent variable and introduce $`\xi =\mathrm{ln}r`$ ($`r>0`$) as a new independent variable. We then make a variable transformation from $`(m,p>0,r>0)`$ to the following three dimensionless variables,

$$u=\frac{4\pi r^3\rho }{m},v=\frac{\rho c^2}{p}\left(\frac{\frac{Gm}{c^2r}}{1\frac{2Gm}{c^2r}}\right),\omega =\omega (p).$$
(5)
The new pressure variable $`\omega (p)`$ is a dimensionless continuous function on $`[0,\mathrm{})`$, strictly monotonically increasing and sufficiently smooth on $`(0,\mathrm{})`$. We require that $`\omega (0)=0`$ and $`\omega \mathrm{}`$ when $`p\mathrm{}`$, and hence $`\omega 0`$. The remaining freedom will be used to adapt $`\omega `$ to features exhibited by the classes of equations of state of interest, cf. Sec. 3. The sign of $`u`$ and $`v`$ depends on whether $`m`$ is positive or negative. In the following we restrict our attention to positive masses, i.e., we investigate a perfect fluid solution only in that range of $`r`$ where it possesses a positive mass-function.888The region where $`m(r)<0`$ can be analyzed with the same dynamical systems methods that are going to be used in the following. The treatment turns even out to be considerably simpler, cf. the case of Newtonian perfect fluids . Moreover, we require that $`12Gm/(rc^2)>0`$, i.e., we consider only solutions to (3) that give rise to a static spacetime metric.999If a solution to (3) satisfies $`12Gm/(rc^2)>0`$ initially at $`r=r_0`$ for initial data $`(m_0,p_0)`$, then this condition holds everywhere. This has been proved (for regular solutions), e.g., in . Within the dynamical systems formulation this result can be established quite easily as we will see in Sec. 4. Solutions violating the condition $`12Gm/(rc^2)>0`$ could be treated with the dynamical systems methods presented in this paper as well. However, we refrain from a discussion of such solutions here. Accordingly, we can assume $`u>0`$ and $`v>0`$.

Starting from (3) the transformation to the new variables $`u,v,\omega `$ yields the following system of equations:

$`{\displaystyle \frac{du}{d\xi }}=u(3u+\mathrm{{\rm Y}}h)`$ (6a)
$`{\displaystyle \frac{dv}{d\xi }}=v[(1\mathrm{{\rm Y}})h+(u1)(1+2\sigma v)]`$ (6b)
$`{\displaystyle \frac{d\omega }{d\xi }}=f\omega h,`$ (6c)
where $`h=h(u,v,\omega )={\displaystyle \frac{d\mathrm{ln}p}{d\mathrm{ln}r}}=v(1+\sigma )(1+\sigma u),`$ (6d)
and $`f=f(\omega )=d\mathrm{ln}\omega /d\mathrm{ln}p|_{p(\omega )}`$. Also, $`\sigma `$ and $`\mathrm{{\rm Y}}`$ are understood as functions of $`\omega `$.

We now proceed by defining new bounded variables $`(U,V,\mathrm{\Omega })(0,1)^3`$,

$$U=\frac{u}{1+u},V=\frac{v}{1+v},\mathrm{\Omega }=\frac{\omega }{1+\omega }.$$
(7)
Introducing a new independent variable $`\lambda `$ according to $`d\lambda /d\xi =(1U)^1(1V)^1`$ yields the dynamical system

$`{\displaystyle \frac{dU}{d\lambda }}=U(1U)[(1V)(34U)\mathrm{{\rm Y}}H]`$ (8a)
$`{\displaystyle \frac{dV}{d\lambda }}=V(1V)[(2U1)(1V+2\sigma V)+(1\mathrm{{\rm Y}})H]`$ (8b)
$`{\displaystyle \frac{d\mathrm{\Omega }}{d\lambda }}=f\mathrm{\Omega }(1\mathrm{\Omega })H,`$ (8c)
where $`H=H(U,V,\mathrm{\Omega })=(1+\sigma )V(1U+\sigma U),`$ (8d)
and where $`\mathrm{{\rm Y}}`$, $`\sigma `$, and $`f`$ are now functions depending on $`\mathrm{\Omega }`$.

It turns out to be useful, indeed essential, to include the boundaries in our analysis so that our compactified state space consists of the unit cube $`[0,1]^3`$. To be able to discuss the dynamical system (8) in a straight forward manner (e.g., to do a fixed point analysis) we require the system to be $`𝒞^1`$-differentiable on $`[0,1]^3`$. This is the case if $`\mathrm{\Omega }(1\mathrm{\Omega })f(\mathrm{\Omega })`$, $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$, and $`\sigma (\mathrm{\Omega })`$ are $`𝒞^1`$ on $`[0,1]`$. In the next section we will show that a broad class of equations of state satisfies these requirements.

## 3 Equations of state

In this section the consequences of the $`𝒞^1`$-differentiability requirement are examined in detail. There are three main building blocks: firstly, the requirement that $`\mathrm{\Omega }(1\mathrm{\Omega })f(\mathrm{\Omega })`$ is $`𝒞^1[0,1]`$ is shown to imply that the choice of $`\omega (p)`$ in (5) is subject to certain restrictions; secondly, the assumptions on $`\mathrm{{\rm Y}}`$ are formulated in a precise way: $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ must be $`𝒞^1`$ in an admissible pressure variable $`\mathrm{\Omega }`$; this assumption defines the classes of equations of state we can treat in the dynamical systems framework; thirdly, $`\sigma `$ is made $`𝒞^1`$ by the freedom in choosing $`\omega `$. The section is concluded by illustrative examples.

Part 1. The condition that $`\mathrm{\Omega }(1\mathrm{\Omega })f(\mathrm{\Omega })`$ is $`𝒞^1[0,1]`$ restricts the choice of the pressure variable $`\omega (p)`$ in (5). This condition is satisfied if and only if $`f(\mathrm{\Omega })`$ is $`𝒞^0[0,1]`$ and $`𝒞^1(0,1)`$, such that $`\mathrm{\Omega }df/d\mathrm{\Omega }`$ and $`(1\mathrm{\Omega })df/d\mathrm{\Omega }`$ vanish in the limit $`\mathrm{\Omega }0`$ and $`\mathrm{\Omega }1`$ respectively. Expressed in terms of $`f(p)`$, a pressure variable $`\omega (p)`$ is admissible, if the conditions $`f(p)𝒞^0[0,\mathrm{})`$, $`lim_p\mathrm{}f(p)<\mathrm{}`$, $`f(p)𝒞^1(0,\mathrm{})`$, and $`lim_{p0}d\mathrm{ln}f/d\mathrm{ln}p=lim_p\mathrm{}d\mathrm{ln}f/d\mathrm{ln}p=0`$ are satisfied. It is not necessary, but convenient, to choose $`\omega `$ in such a way that $`f(\mathrm{\Omega })`$ becomes strictly positive. Thus we require $`f_0=f(0)>0`$ and $`f_1=f(1)>0`$ in the following.101010Note that the Jacobian of the right hand side of (8) contains $`d/d\mathrm{\Omega }[f(\mathrm{\Omega })\mathrm{\Omega }(1\mathrm{\Omega })]`$, which when evaluated at $`\mathrm{\Omega }=0`$ and $`\mathrm{\Omega }=1`$ equals $`f_0`$ and $`f_1`$, respectively. If these numbers are not zero the discussion of the dynamical system becomes easier.

###### Examples.

The simplest example of an admissible pressure variable is provided by $`\omega =(kp)^a`$, where $`a>0`$ is a constant and where $`k`$ is a positive constant carrying the dimension $`[p]^1`$, so that $`\omega `$ becomes dimensionless. In this example $`f(p)=d\mathrm{ln}\omega /d\mathrm{ln}pa`$, which obviously satisfies the imposed conditions. Other simple examples include pressure variables described by $`\omega (p)=(kp)^a(\pm \mathrm{ln}kp)^b`$ for large $`p`$ or small $`p`$, respectively. Here either $`a>0`$ and $`b`$, or $`a=0`$ and $`b<0`$.

Part 2. Let us write the equation of state in the form

$$k\rho c^2=\psi (kp),$$
(9)
where $`k>0`$ is a constant with dimension $`[p]^1`$; hence $`\mathrm{{\rm Y}}(kp)=d\mathrm{ln}\psi (kp)/d\mathrm{ln}(kp)`$. Recall that we assume $`d\rho /dp0`$, i.e., $`\mathrm{{\rm Y}}0`$.

###### Assumptions.

We assume that there exists a variable $`\mathrm{\Omega }`$ (i.e., $`\omega (p)`$) such that $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ is $`𝒞^1[0,1]`$. Moreover, we require $`\mathrm{{\rm Y}}_0=\mathrm{{\rm Y}}(0)<1`$ and $`\mathrm{{\rm Y}}_1=\mathrm{{\rm Y}}(1)=1`$.

###### Remark.

The above assumptions restrict the equations of state to be asymptotically polytropic for low pressures and asymptotically linear at high pressures. Note that if one is only interested in classes of perfect fluid solutions with bounded pressures, $`pp_{\mathrm{max}}`$, then restrictions at $`\mathrm{\Omega }=1`$ are unnecessary, since then $`\mathrm{\Omega }\mathrm{\Omega }_{\mathrm{max}}<1`$. The case $`\mathrm{{\rm Y}}_0=0`$ corresponds to an equation of state that is asymptotically incompressible when $`p0`$. Above we excluded the possibility $`\mathrm{{\rm Y}}_01`$, which includes the asymptotically linear case $`\mathrm{{\rm Y}}_0=1`$. This is not necessary, but models with $`\mathrm{{\rm Y}}_01`$ lead to solutions with infinite masses and radii, and are therefore not particularly interesting from an astrophysical point of view. Note, however, that $`\mathrm{{\rm Y}}(p)`$ might be $`1`$ for some range of $`p`$ (corresponding to a negative polytropic index-function $`n(p)`$), which happens, e.g., for equations of state that cover the phenomenon of neutron drip.

Let us now discuss the consequences of the assumptions for the equation of state. In a neighborhood of $`\mathrm{\Omega }=0`$ let us write $`\mathrm{{\rm Y}}(\mathrm{\Omega })=\mathrm{{\rm Y}}_0+\mathrm{\Omega }\stackrel{~}{\upsilon }(\mathrm{\Omega })`$, or equivalently, $`\mathrm{{\rm Y}}(\omega )=\mathrm{{\rm Y}}_0+\omega \upsilon (\omega )`$. Then the above assumptions can be expressed as follows: $`\upsilon (\omega )`$ is a continuous function, which is $`𝒞^1`$ away from $`\omega =0`$, and satisfies $`lim_{\omega 0}\omega d\upsilon /d\omega =0`$.111111Compare also with the discussion about admissible $`\omega (p)`$. By the definition of $`\mathrm{{\rm Y}}`$ we get $`d\mathrm{ln}(p^{\mathrm{{\rm Y}}_0}\rho )/d\mathrm{ln}p=\omega \upsilon `$, which, using $`f=d\mathrm{ln}\omega /d\mathrm{ln}p`$, gives rise to

$$k\rho c^2=s(kp)^{\mathrm{{\rm Y}}_0}\mathrm{exp}[f^1(\omega )\upsilon (\omega )d\omega ]=:s(kp)^{\mathrm{{\rm Y}}_0}(1+g(\omega )),$$
(10)
where $`s`$ is a dimensionless constant. The function $`g(\omega )`$ is $`𝒞^1`$, and even $`𝒞^2`$ for $`\omega >0`$, $`g(0)=0`$; moreover, $`\omega d^2g/d\omega ^2`$ vanishes as $`\omega 0`$. To first order we have $`g(\omega )=f^1(\omega )\upsilon (\omega )𝑑\omega `$.

Analogously, in a neighborhood of $`\mathrm{\Omega }=1`$, one can show that the assumptions can be translated to

$$k\rho c^2=(\gamma _11)^1(kp)(1+\stackrel{~}{g}(\omega ^1)),$$
(11)
where $`\gamma _1>1`$ is a dimensionless constant, satisfying $`\gamma _12`$ if we require (asymptotic) causality. Hence we have shown that the above assumptions on $`\mathrm{{\rm Y}}`$ are satisfied if and only if there exists a pressure variable $`\omega `$ and functions $`g`$, $`\stackrel{~}{g}`$ with the above properties, such that the equation of state (9) can be written in the asymptotic form (10) and (11), respectively.

###### Examples.

For a simple example consider the pressure variable $`\omega =(kp)^a`$ and take $`g(\omega )=const\omega `$. Then (10) reads $`k\rho c^2=(kp)^{\mathrm{{\rm Y}}_0}(1+const(kp)^a)`$.

###### Remark.

If the assumptions hold, then it is always possible to choose $`\omega `$ in such a way that $`d\mathrm{{\rm Y}}/d\mathrm{\Omega }`$ vanishes as $`\mathrm{\Omega }0`$ and $`\mathrm{\Omega }1`$. We can even achieve

$$\frac{d\mathrm{{\rm Y}}}{d\mathrm{\Omega }}=O(\mathrm{\Omega }^k)(\mathrm{\Omega }0)\text{and}\frac{d\mathrm{{\rm Y}}}{d\mathrm{\Omega }}=O((1\mathrm{\Omega })^k)(\mathrm{\Omega }1)$$
(12)
for an arbitrary $`k>0`$. This is based on the fact that if $`\omega `$ is an admissible pressure variable then so is $`\stackrel{~}{\omega }=\omega ^l`$ ($`0<l<1`$)121212Its logarithmic derivative $`\stackrel{~}{f}=d\mathrm{ln}\stackrel{~}{\omega }/dp`$ satisfies $`\stackrel{~}{f}=lf`$. , and $`d\mathrm{{\rm Y}}/d\stackrel{~}{\mathrm{\Omega }}=O(\stackrel{~}{\mathrm{\Omega }}^{(1l)/l})`$ for $`\stackrel{~}{\mathrm{\Omega }}0`$, and analogously for $`\stackrel{~}{\mathrm{\Omega }}1`$.

Part 3. The third imposed condition is that $`\sigma (\mathrm{\Omega })`$ is $`𝒞^1[0,1]`$. Firstly, the assumptions on $`\mathrm{{\rm Y}}`$ clearly imply that $`\sigma (\mathrm{\Omega })`$ is continuous with $`\sigma _0=\sigma (0)=0`$ and $`\sigma _1=\sigma (1)=\gamma _11`$. Secondly, $`\sigma (\mathrm{\Omega })`$ is $`𝒞^1[0,1]`$, if $`\mathrm{\Omega }`$ is chosen appropriately: a straight forward computation of $`d\sigma /d\mathrm{\Omega }`$ yields $`d\sigma /d\mathrm{\Omega }=\sigma [\mathrm{\Omega }(1\mathrm{\Omega })]^1f^1(1\mathrm{{\rm Y}})`$. For $`\mathrm{\Omega }1`$ we obtain $`lim_{\mathrm{\Omega }1}d\sigma /d\mathrm{\Omega }=(\gamma _11)f_1^1lim_{\mathrm{\Omega }1}d\mathrm{{\rm Y}}/d(1\mathrm{\Omega })`$, hence $`lim_{\mathrm{\Omega }1}d\sigma /d\mathrm{\Omega }`$ exists by the assumptions on $`\mathrm{{\rm Y}}`$. In the limit $`\mathrm{\Omega }0`$ we obtain $`lim_{\mathrm{\Omega }0}d\mathrm{ln}\sigma /d\mathrm{ln}\mathrm{\Omega }=lim_{\mathrm{\Omega }0}(\mathrm{\Omega }/\sigma )d\sigma /d\mathrm{\Omega }=f_0^1(1\mathrm{{\rm Y}}_0)`$. If $`f_0^1(1\mathrm{{\rm Y}}_0)>1`$, then $`lim_{\mathrm{\Omega }0}d\sigma /d\mathrm{\Omega }`$ exists and $`lim_{\mathrm{\Omega }0}d\sigma /d\mathrm{\Omega }=0`$. If $`f_0^1(1\mathrm{{\rm Y}}_0)<1`$, then $`d\sigma /d\mathrm{\Omega }\mathrm{}`$ for $`\mathrm{\Omega }0`$.131313If $`f_0=(1\mathrm{{\rm Y}}_0)`$, then in many cases $`lim_{\mathrm{\Omega }0}d\sigma /d\mathrm{\Omega }`$ exists. Note, however, that this does not hold in general. Therefore, provided that $`\mathrm{\Omega }`$ is chosen so that $`f_0<(1\mathrm{{\rm Y}}_0)=1/(1+n_0)`$, we conclude that $`\sigma (\mathrm{\Omega })𝒞^1[0,1]`$. In addition, in analogy with the discussion involving $`\mathrm{{\rm Y}}`$, we can always achieve

$$\frac{d\sigma }{d\mathrm{\Omega }}=O(\mathrm{\Omega }^k)(\mathrm{\Omega }0)\text{and}\frac{d\sigma }{d\mathrm{\Omega }}=O((1\mathrm{\Omega })^k)(\mathrm{\Omega }1)$$
(13)
for arbitrarily large $`k`$.

###### Remark.

Consider the equation of state (9). We choose $`\omega `$ as a function of the dimensionless $`kp`$ as already indicated above, i.e., $`\omega =\omega (kp)`$. Then the equation of state can be written in the implicit form

$$kp=\varphi (\omega ),k\rho c^2=\chi (\omega ),$$
(14)
where $`\varphi `$ and $`\chi `$ may contain an arbitrary number of dimensionless parameters. The derived quantities $`f`$, $`\mathrm{{\rm Y}}`$, and $`\sigma `$ then read $`f(\omega )=(d\mathrm{ln}\varphi /d\mathrm{ln}\omega )^1`$, $`\mathrm{{\rm Y}}(\omega )=(d\mathrm{ln}\chi /d\mathrm{ln}\omega )(d\mathrm{ln}\varphi /d\mathrm{ln}\omega )^1`$, and $`\sigma (\omega )=\varphi (\omega )/\chi (\omega )`$. Note that $`f(\omega )`$, $`\mathrm{{\rm Y}}(\omega )`$, $`\sigma (\omega )`$ are functions independent of the dimensional parameter $`k`$; therefore the whole class of equations of state (14), or equivalently (9), parameterized by $`k`$ is described by one dynamical system (8) with specified $`f(\mathrm{\Omega })`$, $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$, and $`\sigma (\mathrm{\Omega })`$.

###### Example.

Relativistic polytropes. The relativistic generalization of the Newtonian adiabatic index $`\mathrm{\Gamma }_N`$ (as defined in (4)) is the adiabatic index

$$\mathrm{\Gamma }(p):=\frac{1}{c^2}\frac{p+\rho c^2}{p}\frac{dp}{d\rho }=(1+\sigma )\mathrm{\Gamma }_N(p).$$
(15)
Clearly, $`\mathrm{\Gamma }`$ reduces to $`\mathrm{\Gamma }_N`$ when $`c\mathrm{}`$.141414Recall, however, that in the Newtonian case $`\rho `$ stands for the rest-mass density. Regarded as a function of $`\mathrm{\Omega }`$, the adiabatic index $`\mathrm{\Gamma }`$ becomes a $`𝒞^1[0,1]`$ function for the considered class of asymptotically polytropic/asymptotically linear equations of state. In particular, $`\mathrm{\Gamma }(0)=\mathrm{\Gamma }_N(0)=\mathrm{{\rm Y}}_0^1=1+1/n_0`$, and $`\mathrm{\Gamma }(1)=1+\sigma _1=\gamma _1`$. Whereas $`\mathrm{\Gamma }_N=const`$ defines the polytropic equations of state, $`\mathrm{\Gamma }=const`$ gives rise to the so-called relativistic polytropes, usually written as

$$\rho c^2=(p/K)^{1/\mathrm{\Gamma }}c^2+p/(\mathrm{\Gamma }1),$$
(16)
where $`K`$ is a positive parameter. Requiring $`d\rho /dp0`$ and causality of the speed of sound $`c_s=\sqrt{dp/d\rho }c`$ implies $`1<\mathrm{\Gamma }2`$, where $`\mathrm{\Gamma }=2`$ corresponds to the case $`c_sc`$ when $`p\mathrm{}`$. The class (16) of equation of state can be expressed nicely in the form (9), if we take $`k=\left(c^2(\mathrm{\Gamma }1)K^{1/\mathrm{\Gamma }}\right)^{\mathrm{\Gamma }/(\mathrm{\Gamma }1)}`$:

$$k\rho c^2=\left(\frac{1}{\mathrm{\Gamma }1}\right)\left((kp)^{1/\mathrm{\Gamma }}+kp\right).$$
(17)
Obviously, the relativistic polytropes are of the form (10) and (11) by the simple choice $`\omega =(kp)^a`$ with $`a(\mathrm{\Gamma }1)/\mathrm{\Gamma }`$, whereby $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ and $`\sigma (\mathrm{\Omega })`$ are $`𝒞^1[0,1]`$. To simplify the expressions for $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ and $`\sigma (\mathrm{\Omega })`$ as much as possible we take $`\omega =(kp)^{(\mathrm{\Gamma }1)/\mathrm{\Gamma }}`$, i.e., $`f=(\mathrm{\Gamma }1)/\mathrm{\Gamma }`$, which results in $`\sigma (\mathrm{\Omega })=(\mathrm{\Gamma }1)\mathrm{\Omega }`$ and $`\mathrm{{\rm Y}}(\mathrm{\Omega })=(1+(\mathrm{\Gamma }1)\mathrm{\Omega })/\mathrm{\Gamma }`$. Hence the right hand side of the dynamical system (8) consists of pure polynomials. Note that only one (dimensionless) parameter appears in the dynamical system, namely $`\mathrm{\Gamma }`$, i.e., the entire one-parameter family parameterized by $`k`$ is represented by a single dynamical system, as expected.

###### Example.

Ideal neutron gas. The equation of state of a degenerate ideal neutron gas is implicitly given by

$`kp`$ $`={\displaystyle \frac{1}{8\pi ^2}}\left(x(1+x^2)^{1/2}(\frac{2}{3}x^21)+\mathrm{ln}[x+(1+x^2)^{1/2}]\right)`$ (18a)
$`k\rho c^2`$ $`={\displaystyle \frac{1}{8\pi ^2}}\left(x(1+x^2)^{1/2}(1+2x^2)\mathrm{ln}[x+(1+x^2)^{1/2}]\right),`$ (18b)
where $`x`$ is essentially the Fermi momentum.151515Namely, $`x=p_F/m_nc`$, where $`m_n`$ is the neutron mass. The constant $`k`$ is given by $`k=\lambda _n^3/m_nc^2`$, where $`\lambda _n=\mathrm{}/m_nc`$, for details see, e.g., , p. 23ff. A straight forward way to treat (18) is to choose $`\omega =x`$, whereby $`\mathrm{\Omega }=x/(1+x)`$: Firstly, $`x`$ is an admissible pressure variable, as $`f(\mathrm{\Omega })`$ is a sufficiently smooth function satisfying $`f=1/5+O(\mathrm{\Omega }^2)`$ ($`\mathrm{\Omega }0`$) and $`f=1/4+O((1\mathrm{\Omega })^2)`$ ($`\mathrm{\Omega }1`$). Secondly, $`\sigma (\mathrm{\Omega })`$ is smooth with $`\sigma =\mathrm{\Omega }^2/5+O(\mathrm{\Omega }^4)`$ and $`\sigma =1/3+O((1\mathrm{\Omega })^2)`$ as $`\mathrm{\Omega }0,1`$ respectively. In particular, $`\sigma _1=1/3`$, i.e., $`\gamma _1=4/3`$. Thirdly, $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ interpolates smoothly between $`\mathrm{{\rm Y}}_0=3/5`$, i.e., $`n_0=3/2`$, and $`\mathrm{{\rm Y}}_1=1`$. Moreover, $`d\mathrm{{\rm Y}}/d\mathrm{\Omega }=O(\mathrm{\Omega })`$ and $`d\mathrm{{\rm Y}}/d\mathrm{\Omega }=O(1\mathrm{\Omega })`$ for $`\mathrm{\Omega }0,1`$. Therefore, the ideal neutron gas is easily described in the dynamical systems formalism, moreover, our formalism also covers more sophisticated equations of state that are used to model neutron stars.

To conclude this section we outline how phase transitions and composite equations of state can be incorporated into the dynamical systems framework. Consider an equation of state $`\rho =\rho (p)`$ that is piecewise continuous on $`[0,\mathrm{})`$ and piecewise sufficiently smooth on $`(0,\mathrm{})`$. This implies that at certain values $`p_j`$ of the pressure, $`\rho (p)`$ makes a jump or a kink, corresponding to a phase transition of first or second kind; $`\mathrm{{\rm Y}}(p)`$ and $`\sigma (p)`$ exhibit similar behavior. The associated “weak solutions” $`p(r)`$ and $`m(r)`$ of equations (3) are required to be continuous, while $`\rho (r)`$ is continuous only if $`\rho (p)`$ is. Also the variables $`u,v`$ and $`U,V`$ are continuous functions if $`\rho `$ is; if $`\rho `$ jumps from some value $`\rho _1`$ to $`\rho _2`$ at $`p_j`$, then $`u,v`$ make a jump according to $`u_1/u_2=v_1/v_2=\rho _1/\rho _2`$, and analogously for $`U,V`$. Choosing a conventional continuous pressure variable $`\omega (p)`$ we obtain piecewise continuous/smooth orbits, monotonic in $`\mathrm{\Omega }`$, that have possible jumps in $`U,V`$ at $`\mathrm{\Omega }_j=\mathrm{\Omega }(p_j)`$.

For practical purposes a slightly different viewpoint turns out to be more convenient. Assume for simplicity that there is only one jump of $`\rho (p)`$, at $`p_j`$. Choose two smooth equations of state $`\rho _I(p)`$, $`\rho _{II}(p)`$, such that $`\rho \rho _I`$ for $`p<p_j`$ and $`\rho \rho _{II}`$ for $`p>p_j`$. Then the equation of state $`\rho `$ is simply obtained by switching from $`\rho _I`$ to $`\rho _{II}`$ at $`p_j`$. Correspondingly, we deal with two different state spaces, in particular we can have two different pressure variables $`\mathrm{\Omega }_I`$ and $`\mathrm{\Omega }_{II}`$. ¿From this point of view, a fluid solutions associated with $`\rho `$ starts in one state space, and ends in the other state space; the jump in between appears as a map from one state space to the other; besides the jump in $`U,V`$ described above it entails a jump in $`\mathrm{\Omega }`$, from $`\mathrm{\Omega }_I(p_j)`$ to $`\mathrm{\Omega }_{II}(p_j)`$. An explicit example will be given in Sec. 7.

## 4 Dynamical systems analysis

In this section we study the dynamical system (8); in particular we investigate the global dynamics. The arena for the analysis is the state space, i.e., the unit cube $`[0,1]^3`$, which is endowed with different vector fields and associated flows that depend on the equation of state.

Central to the dynamical systems analysis are the equilibrium points of the system. The fixed points as well as their associated eigenvalues161616The eigenvalues of the linearization of the system at the fixed point (together with the eigenvectors) characterize the flow of the dynamical system in a neighborhood of the fixed point (Hartman-Grobman theorem). For an introduction to the theory of dynamical systems, see e.g. . are listed in Table 1. Note that the equation of state enters only via its asymptotic properties: both the location of the fixed points and the eigenvalues depend only on $`n_0`$ and $`\gamma _1`$. The state space and the fixed points are depicted in Fig. 1.

We observe that all six faces of the cube are invariant subspaces. On each of the four side faces, the induced dynamical system possesses a simple structure; the orbits on the side faces171717Note that whereas the orbits in the interior of the cube represent perfect fluid solutions, the orbits on the side faces cannot be interpreted in physical terms. This is essentially because the variable transformation (5) cannot be inverted for, e.g., $`u=0`$ or $`v=0`$. In a modified dynamical systems formulation, where $`\omega `$ is replaced by a variable depending on $`\mathrm{\Phi }`$, certain side faces may be interpreted as representing vacuum solutions. are depicted in Fig. 2.

The dynamical system on the $`\mathrm{\Omega }=0`$ subset is in fact Newtonian, cf. Appendix A: since $`\sigma =0`$ on $`\mathrm{\Omega }=0`$, $`U`$ and $`V`$ coincide with Newtonian homology variables and thus the system represents the Newtonian homology invariant equations for an exact polytrope with index $`n_0`$. On the $`\mathrm{\Omega }=0`$ subset the flow of the dynamical system exhibits rather complicated features. In particular, $`n_0`$ serves as a bifurcation parameter: there exist bifurcations for $`n_0=0`$, $`n_0=3`$, and $`n_0=5`$, as indicated by the properties of the fixed points. The case $`n_0=5`$ is exactly solvable: the orbits can be represented implicitly by

$$\mathrm{\Phi }=\sqrt{U}V^{3/2}(U(89V)+7V6)d(1U)^{3/2}(1V)^{5/2}=0,$$
(19)
where the parameter $`d`$ is required to satisfy $`d9/8`$. The fixed point $`B_3`$ is represented by $`d=9/8`$; $`9/8<d<0`$ defines a one-parameter family of closed orbits, which we will henceforth denote by $`C`$; $`d=0`$ characterizes the orbit that connects $`L_2`$ and $`B_1`$; the remaining orbits correspond to positive values of $`d`$, cf. Fig. 2(i). Typical orbits for a selection of values of $`n_0`$ are depicted in Fig. 2; cf. also the analogous non-compact figures in the standard literature, e.g., p.201 in .

The dynamical system on the $`\mathrm{\Omega }=1`$ subset can be regarded as describing the relativistic scale-invariant equations for linear equations of state $`p=(\gamma _11)\rho c^2`$, a fact that reflects itself in the decoupling of the $`(U,V)`$-equations from the system (8) for purely linear equations of state. Typical orbits are depicted in Fig. 2 (cf. also, e.g., , where an analogous non-compact phase portrait is given).

The global dynamics of the dynamical system (8) is determined by the fixed points on the boundaries of the state space, which enables us to describe the asymptotics of the interior orbits. The statement is made precise in the following theorem.

###### Theorem 4.1.

The $`\alpha `$-limit181818For a dynamical system $`\dot{x}=f(x)`$, the $`\omega `$-limit set ($`\alpha `$-limit) of a point $`x`$ is defined as the set of all accumulation points of the future (past) orbit of $`x`$. Limit sets thus characterize the asymptotic behavior of the dynamical system. of an interior orbit is a fixed point on $`L_1`$, $`L_2`$, or the fixed point $`T_3`$. The $`\omega `$-limit of an interior orbit is always located on $`\mathrm{\Omega }=0`$, and it is one of the fixed points $`B_1`$, $`B_2`$, or $`B_4`$, when $`n_05`$. If $`n_0=5`$, then the $`\omega `$-limit can also be an element of the 1-parameter set of closed orbits $`C`$.

###### Proof.

The proof of the theorem is based on the monotonicity principle : if there exists a function $`Z`$ (defined on a closed, bounded, future (past) invariant set $`S`$) that is monotonic along the flow of the dynamical system $`\dot{x}=f(x)`$, then $`\omega (S)`$ ($`\alpha (S)`$) is contained in $`\{xS|\dot{Z}=0\}`$.

The function $`\mathrm{\Omega }`$ is such a monotone function on the state space, i.e., $`d\mathrm{\Omega }/d\lambda <0`$, except on $`\mathrm{\Omega }=0`$, $`\mathrm{\Omega }=1`$, and $`V=0`$. It follows that the $`\omega `$-limit of an interior orbit is located on $`\mathrm{\Omega }=0`$ or $`V=0`$, while the $`\alpha `$-limit of an interior orbit must lie on $`\mathrm{\Omega }=1`$ or $`V=0`$. However, it follows from the known orbit structure on $`V=0`$ and the local fixed point analysis that $`V=0`$ cannot possess an $`\omega `$-limit. Similarly, it follows that only $`L_1`$ and $`L_2`$ constitute $`\alpha `$-limits for interior orbits.

It remains to determine the $`\omega `$-limit ($`\alpha `$-limit) sets for interior orbits on $`\mathrm{\Omega }=1`$ and $`\mathrm{\Omega }=0`$. Again we make use of the monotonicity principle, however, to simplify the discussion we use the uncompactified variables $`(u,v)`$.

On $`\mathrm{\Omega }=1`$, a monotone function is given by191919This monotone function was found through Hamiltonian methods, see  and chapter 10 in .

$$Z_3=[(\gamma _11)h_1(u,v)+3\gamma _12]^2(uv)^{p_1}(1+2(\gamma _11)v)^{p_11},$$
(20)
where the non-negative function $`h_1`$ and the non-negative constant $`p_1`$ are defined by

$$h_1(u,v):=\gamma _1(1+(\gamma _11)u)v=h|_{\mathrm{\Omega }=1},p_1=\frac{4(\gamma _11)^2}{\gamma _1(5\gamma _14)}.$$
(21)
Its derivative along orbits reads

$$\frac{dZ_3}{d\xi }=\frac{p_1}{2}(3\gamma _12)(2h_1)^2[(\gamma _11)h_1+3\gamma _12]^1Z_3.$$
(22)
Application of the monotonicity principle yields that only $`T_3`$ is a possible $`\alpha `$-limit set on $`\mathrm{\Omega }=1`$.

We now turn our attention to the $`\omega `$-limit sets of interior orbits, which we have shown to lie on $`\mathrm{\Omega }=0`$. For $`n_03`$ the function $`Z_1=uv^3`$ is monotone along orbits, $`dZ_1/d\xi =Z_1(2u+v(34\mathrm{{\rm Y}}_0))`$. For $`n_0>3`$ define

$$Z_2=(uv)^{2/(n_01)}\left[\frac{n_05}{(n_01)^2}+\frac{2v}{(n_0+1)(n_01)}\frac{v^2}{2(n_0+1)^2}\frac{uv}{(n_0+1)^2}\right],$$
(23)
and observe that $`dZ_2/d\xi =(n_05)(n_0+1)^2(n_01)^3((uv)^{2/(n_01)}(2(n_0+1)(n_01)v)^2`$, i.e., $`Z_2`$ is monotonic on orbits for $`3<n_05`$. (Note that the flow is non-vanishing on $`v=2(n_0+1)/(n_01)`$ except for at $`B_4`$.) Thus, by the monotonicity principle, taking also the solution structure on the boundaries into account, only the fixed points $`B_1`$, $`B_2`$, $`B_4`$ can serve as $`\omega `$-limit sets on $`\mathrm{\Omega }=0`$. The exceptional case is $`n_0=5`$: then there also exists a 1-parameter family, $`C`$, of closed orbits (periodic solutions) described by (23) (see also (19)), which acts as an $`\omega `$-limit set for interior solutions.

###### Theorem 4.2.

If $`0<n_03`$ $`(n_0=0)`$, then all orbits end at $`B_2`$ $`(L_4)`$.

###### Proof.

We distinguish two cases, $`n_0<3`$ and $`n_0=3`$. For $`n_0<3`$ the proof is trivial; each interior orbit converges to $`B_2`$ for $`0<n_0<3`$ (to $`L_4`$ when $`n_0=0`$), which follows from the proof of Theorem 4.1 and the local dynamical systems analysis (Table 1).

In the case $`n_0=3`$, $`B_1`$ is not hyperbolic. Nevertheless, by applying center manifold theory, we show that no interior orbit converges to $`B_1`$ as $`\lambda \mathrm{}`$:

It is advantageous to investigate the problem in adapted uncompactified variables $`\{x0,y,z0\}`$, defined by

$$u=x(n_02)/(n_0+1)v=(n_0+1)x+y\omega =z,$$
(24)
The fixed point $`B_1`$ is represented by $`(0,0,0)`$ and the dynamical system (6) reads

$$\begin{array}{ccc}\dot{x}=\hfill & & +N_x(x,y,z)\hfill \\ \dot{y}=\hfill & y& +N_y(x,y,z)\hfill \\ \dot{z}=\hfill & f_0(1+n_0)z& +N_z(x,y,z),\hfill \end{array}$$
(25)
where $`N_x`$, $`N_y`$, $`N_z`$ denote the nonlinear terms. $`N_x`$ is given by

$$N_x=x(34\mathrm{{\rm Y}}(z)x/4+(xy)\mathrm{{\rm Y}}(z))+C_1(x,y,z)\sigma (z)+C_2(x,y,z)\sigma ^2(z).$$
(26)
A center manifold of the system (25), i.e., an invariant manifold tangential to the center subspace at the fixed point, is represented by $`\{(x,h_y(x),h_z(x))|0x<\epsilon \}`$ satisfying

$$_xh_y(x)\dot{x}(x,h_y(x),h_z(x))=\dot{y}(x,h_y(x),h_z(x))_xh_z(x)\dot{x}(x,h_y(x),h_z(x))=\dot{z}(x,h_y(x),h_z(x)),$$
(27)
and the tangency conditions $`h(0)=0`$, $`_xh(0)=0`$. We obtain

$$h_y(x)=x^2/2x^3+O(x^4)h_z(x)0.$$
(28)
Note in particular that the center manifold lies in the Newtonian subset of the state space. The center manifold reduction theorem states that the flow of the full nonlinear system is locally equivalent to the flow of the reduced system (see, e.g., )

$$\begin{array}{ccc}\dot{x}\hfill & =& N_x(x,h_y(x),h_z(x))=x^2/2+O(x^3)\hfill \\ \dot{y}\hfill & =& y\hfill \\ \dot{z}\hfill & =& f_0(1+n_0)z.\hfill \end{array}$$
(29)
The flow given by (29) clearly prevents interior solutions $`(x>0,y,z>0)`$ from converging to $`B_1`$, which proves the claim.

Collecting the results of the local and global dynamical systems analysis yields the following. Interior orbits originate from $`L_1`$ (a two-parameter set); $`L_2`$ (a one-parameter set); and $`T_3`$ (a single orbit). They end at $`B_1`$ (when $`n_0>3`$; a one-parameter set); $`B_2`$ (a two-parameter set); $`L_4`$ (when $`n_0=0`$; a two-parameter set); $`B_4`$ (when $`n_0>3`$; a single orbit when $`n_05`$ and a two-parameter set when $`n_0>5`$); and at the periodic curves $`C`$ (when $`n_0=5`$; a two-parameter set).

An orbit in the state space corresponds to a perfect fluid solution. Note that the perfect fluid solution is only represented in that range of $`r`$ where $`p>0`$, i.e., the fluid body, and where $`m>0`$; cf. Sec. 2. If the configuration has a finite radius $`R`$, then the fluid solution is joined to a vacuum solution, a Schwarzschild solution, at $`r=R`$.

The one-parameter set of regular perfect fluid solutions, i.e., solutions of (3) with a regular center of spherical symmetry, appears in the state space as the one-parameter set of orbits that originate from $`L_2`$. When $`d\rho /dp0`$ is assumed, $`\rho (r)\overline{\rho }(r)`$ holds for regular solutions, where $`\overline{\rho }=3m/(4\pi r^3)`$ is the “average density”. Translated to the state space variables we obtain $`U3/4`$. In order to show this inequality within the dynamical systems framework, consider the surface $`U=3/4`$ and the flow through this surface: $`dU/d\lambda |_{U=3/4}=3\mathrm{{\rm Y}}(1+\sigma )(1+3\sigma )V/640`$, since $`\mathrm{{\rm Y}}0`$ when $`d\rho /dp0`$. Hence, $`U=3/4`$ acts as a “semi-permeable membrane” for the flow of the dynamical system. By noting that the unstable subspaces at $`L_4`$ are characterized by $`U3/4`$ the claim is established.202020We observe that in the incompressible case, $`\rho =const`$ (i.e., $`\mathrm{{\rm Y}}=0`$), the regular solutions are located on $`U=3/4`$. The equations are explicitly solvable in this case. Note, however, that we have to cut the state space at some value of $`\mathrm{\Omega }`$, because $`\mathrm{{\rm Y}}=0`$ is not compatible with the requirement $`\mathrm{{\rm Y}}_1=1`$ to obtain a $`𝒞^1`$ dynamical system up to and including $`\mathrm{\Omega }=1`$.

If a perfect fluid solution satisfies $`12Gm/(rc^2)>0`$ initially at $`r=r_0`$ for initial data $`(m_0,p_0)`$, then this condition holds everywhere. For a proof (involving regular solutions) see, e.g., . In the state space picture this result can be established quite easily. By construction, the initial data is a point in the interior of the state space and the orbit passing through this point represents the associated perfect fluid solution. Since $`V<1`$ for the orbit, we obtain that

$$1\frac{2Gm}{rc^2}=\frac{1V}{1V+2\sigma V}>0$$
(30)
holds everywhere within the fluid body. We can even show that $`12Gm/(rc^2)const>0`$ for all $`r`$: the considered orbit can have three possible end points; suppose that the orbit ends in the fixed point $`B_2`$, which is the only non-trivial case since then $`V1`$ when $`\lambda \mathrm{}`$. We will see in the subsequent section that such an orbit gives rise to a perfect fluid solution possessing a finite radius $`R`$ and a total mass $`M`$, which are related by $`GM/(Rc^2)=B/(2B+s)`$, where $`B`$ and $`s`$ are positive numbers. Hence the claim is established.

By introducing the function

$$\mathrm{\Psi }_1:=1\frac{(3+𝒫)^2}{2(2+𝒫)},$$
(31)
where

$$:=\frac{Gm}{c^2r}=\frac{\sigma (\mathrm{\Omega })V}{1V+2\sigma V};𝒫:=\frac{4\pi G}{c^4}r^2p=\left(\frac{U}{1U}\right)\left(\frac{\sigma ^2V}{1V+2\sigma V}\right),$$
(32)
we can write the relativistic Buchdahl inequality as

$$\mathrm{\Psi }_10.$$
(33)
Note that the surface inequality $`4/9`$ is obtained by setting $`𝒫=0`$. The flow on the surface $`\mathrm{\Psi }_1=0`$ is given by

$$\frac{d\mathrm{\Psi }_1}{d\lambda }=\frac{\sigma (34U)(1U+\sigma U)^3V}{(2(1U)+\sigma U)^2(3(1U)+\sigma U)}$$
(34)
and hence, since $`U3/4`$ for regular solutions, it follows that the derivative is positive (zero) when $`U<3/4`$ ($`U=3/4`$). The Buchdahl surface $`\mathrm{\Psi }_1=0`$ and the $`U=3/4`$ surface are depicted in Fig. 3 for the relativistic polytropes with $`\mathrm{\Gamma }=3/4`$ and $`\mathrm{\Gamma }=2`$ (note that the Buchdahl-surface is affected by the equation of state through $`\sigma (\mathrm{\Omega })`$). Since regular solutions start at $`V=0`$ and since the flow on the Buchdahl surface is directed toward $`V=0`$, the regular solutions cannot pass through the Buchdahl-surface. Hence, for regular solutions only a limited part of the state space is accessible (cf. Fig. 3). This constitutes the dynamical systems proof of the Buchdahl inequality.

## 5 Translating the state space picture to a physical picture

In this section we translate the state space picture to the common physical variables. To this end recall, firstly, that an orbit in the interior of the state space stands for a fluid solution $`(m(r),p(r),\rho (r))`$ ($`r>0`$) inside the fluid body via the transformation (5). Secondly, consider the following auxiliary equations:

$`{\displaystyle \frac{dr}{d\lambda }}=(1U)(1V)r,{\displaystyle \frac{dm}{d\lambda }}=U(1V)m`$ (35a)
$`{\displaystyle \frac{dp}{d\lambda }}=(1+\sigma )V(1U+\sigma U)p,{\displaystyle \frac{d\sigma }{d\lambda }}=(1+\sigma )(1\mathrm{{\rm Y}})V(1U+\sigma U)\sigma `$ (35b)
$`{\displaystyle \frac{Gm}{c^2r}}={\displaystyle \frac{\sigma V}{1V+2\sigma V}}`$ (35c)
$`r^2={\displaystyle \frac{c^2}{4\pi G}}\left({\displaystyle \frac{U}{1U}}\right)\left({\displaystyle \frac{\sigma V}{1V+2\sigma V}}\right)\rho (\mathrm{\Omega })^1`$ (35d)
$`m^2={\displaystyle \frac{c^6}{4\pi G^3}}\left({\displaystyle \frac{U}{1U}}\right)\left({\displaystyle \frac{\sigma V}{1V+2\sigma V}}\right)^3\rho (\mathrm{\Omega })^1.`$ (35e)
Based on an equation of state in explicit form (9), i.e., $`k\rho c^2=\psi (kp)`$, and a pressure variable $`\omega =\omega (kp)`$ we can replace $`\rho ^1`$ in the last two equations by

$$\rho ^1=\frac{kc^2}{\psi (\mathrm{\Omega })}\text{or}\rho ^1=kc^2\sigma (\mathrm{\Omega })\frac{1}{(kp)(\mathrm{\Omega })},$$
(36)
where $`(kp)(\mathrm{\Omega })=\mathrm{\Omega }^{1/a}/(1\mathrm{\Omega })^{1/a}`$ in the case $`\omega =(kp)^a`$.

Equations (35a) can be used to obtain an intuitive picture of where solutions gain mass and radius in the state space. Since these equations are independent of $`\mathrm{\Omega }`$, it is possible to visualize $`d\mathrm{ln}m/d\lambda `$ and $`d\mathrm{ln}r/d\lambda `$ as contour plots on the $`U,V`$ plane, see Fig. 4. The picture suggests that no mass is acquired near the side faces $`U=0`$ and $`V=1`$, and hence the fixed points $`B_1`$ and $`B_2`$ (and $`L_5`$, if present) are attractors for solutions with finite masses. In contrast, $`d\mathrm{ln}r/d\lambda `$ does not vanish near $`B_1`$, but only near $`B_2`$ (and $`L_5`$, if present). Therefore, $`B_2`$ is the only attractor for solutions with finite radii. Solutions converging to $`B_3`$ (or to a periodic orbit $`C`$ when $`n_0=5`$) for $`\lambda \mathrm{}`$ acquire both infinite masses and radii.

To obtain a rigorous derivation of the asymptotic mass and radius features of the various solutions we must combine equations (35) with the local dynamical systems analysis in the neighborhood of the fixed points (and the periodic orbits $`C`$ when $`n_0=5`$). We only state the results here (for a detailed discussion of the corresponding and rather similar Newtonian situation, see ):

Solutions that originate from: (i) $`L_1`$ are solutions with negative masses for small $`r`$; when they have acquired sufficient positive mass, so that the mass function $`m`$ has become zero, they appear in the state space through $`L_1`$, cf. the Newtonian case discussed in ); (ii) $`L_2`$ have a regular center and are therefore solutions of particular physical interest; (iii) $`T_3`$ have a non-regular center ($`\rho r^2`$ for $`r0`$) and describe the limit when the central pressure goes to infinity. One can show that $`T_3`$ is associated with a self-similar solution admitting a $`H_4`$ spacetime transitive symmetry group. It is a special case of a solution found by Tolman , but has been rediscovered many times, see  and references therein.

Solutions that end at: (i) $`B_1`$ ($`n_0>3`$) have infinite radii but finite masses; (ii) $`B_2`$ and (iii) $`L_4`$ ($`n_0=0`$) have finite masses and radii and are therefore the most interesting solutions; (iv) $`B_4`$ ($`n_0>3`$) and (v) $`C`$ ($`n_0=5`$) have infinite radii and masses ($`m(r)r^{(n_03)/(n_01)}`$ in the case of $`B_4`$).

Let us now take a closer look at the 2-parameter family of orbits that converges to $`B_2`$ as $`\lambda \mathrm{}`$. All such orbits correspond to perfect fluid solutions with finite radii $`R`$ and masses $`M`$. To describe the behavior of the physical observables as $`rR`$, we will use the local dynamical systems results together with equations (5), (7), or (35). Recall that the solutions of the dynamical system (8) with specified $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ and $`\sigma (\mathrm{\Omega })`$ represent the totality of perfect fluid solutions corresponding to a one-parameter family of equations of state, cf. (14). For $`p0`$ the leading term of a considered equation of state was shown to be given by $`k\rho c^2=s(kp)^{\mathrm{{\rm Y}}_0}`$, where $`s`$ is a dimensionless constant, cf. (10).

The 2-parameter family of orbits that converges to $`B_2`$ can be conveniently characterized by the constants $`A`$ and $`B`$, defined according to

$$A:=\underset{\lambda \mathrm{}}{lim}U(1V)^{n_0},B:=\underset{\lambda \mathrm{}}{lim}\mathrm{\Omega }^{\frac{1}{f_0(n_0+1)}}(1V)^1.$$
(37)
Expressed in these quantities we obtain the following, when $`rR`$,

$`p(r)`$ $`=k^1B^{n_0+1}\delta r^{n_0+1}+{\displaystyle \frac{n_0+1}{n_0+2}}k^1B^{n_0+1}\left(2(n_0+1)+{\displaystyle \frac{4\pi R^3}{Mc^2}}k^1B^{n_0+1}A^1(3+2n_0)\right)\delta r^{n_0+2}+O()`$ (38a)
$`m(r)`$ $`=M(1A\delta r^{n_0+1})+O()`$ (38b)
$`\rho (r)`$ $`={\displaystyle \frac{1}{4\pi R^3}}MA\delta r^{n_0}+O(),`$ (38c)
where $`\delta r=\frac{1}{n_0+1}\frac{Rr}{R}`$.

The radius $`R`$ and mass $`M`$ for a solution are uniquely determined by the values of $`A`$ and $`B`$ through

$$R^2=\frac{kc^4}{4\pi sG}\frac{AB^{1n_0}}{s+2B},M^2=\frac{kc^8}{4\pi sG^3}\frac{AB^{3n_0}}{(s+2B)^3}.$$
(39)
The dimensionless quotient $`GM/(c^2R)`$ is given by $`GM/(c^2R)=B/(s+2B)`$, so that the line element (1) on the surface $`r=R`$ is determined by $`e^{2\mathrm{\Lambda }}=c^2e^{2\mathrm{\Phi }}=1/(1\frac{2GM}{c^2R})=1+2B/s`$. (Hence, naturally, $`GM/(c^2R)`$ can be used instead of $`B`$ together with $`A`$ to parameterize the different solutions.)

Eq. (39) and the related formulas are useful in, e.g., numerical computations, since it is fairly easy to compute $`A`$ and $`B`$. Below, in numerical applications, we are going to compare solutions with a reference solution, i.e., we are going to consider the dimensionless ratios $`R/R_{\mathrm{ref}}`$, $`M/M_{\mathrm{ref}}`$, where $`R_{\mathrm{ref}}`$ and $`M_{\mathrm{ref}}`$ are the radius and mass of a typical reference solution. Of course, the dimensional factors $`(kc^4)/(4\pi sG)`$ and $`(kc^8)(4\pi sG^3)`$ drop out in this case.

For illustrative purposes consider cylindrical coordinates at $`B_2`$, i.e., the coordinate transformation $`(U,V,\mathrm{\Omega })=(\epsilon \mathrm{sin}\varphi ,1\epsilon \mathrm{cos}\varphi ,h^a)`$ with $`0\varphi \pi /2`$, $`0h`$, $`0\epsilon `$, which leads to

$$A=\underset{\lambda \mathrm{}}{lim}ϵ^{1n_0}\mathrm{sin}\varphi (\mathrm{cos}\varphi )^{n_0},B=\underset{\lambda \mathrm{}}{lim}h^{1/(n_0+1)}(ϵ\mathrm{cos}\varphi )^1.$$
(40)
As every orbit in a neighborhood of $`B_2`$ converges to $`B_2`$ as $`\lambda \mathrm{}`$, we obtain approximate expressions with arbitrary accuracy by choosing $`\epsilon `$ sufficiently small. In this picture, every orbit is uniquely characterized by its intersection point $`(h,\varphi )`$ with the small $`\epsilon `$-cylinder, and accordingly these values uniquely determine $`M`$ and $`R`$ by combining (39) with (40). A similar discussion holds for $`L_4`$ when $`n_0=0`$; for the corresponding Newtonian discussion, see .

## 6 Mass-radius theorems

In this section we will prove several theorems concerning mass-radius properties of solutions, where we focus on regular solutions. The underlying equations of state are as always understood to be asymptotically polytropic for $`p0`$ and asymptotically linear for $`p\mathrm{}`$, as discussed in Sec. 3.

###### Theorem 6.1.

All regular solutions have infinite masses and infinite radii if $`\mathrm{\Gamma }_N\frac{6}{5}`$ and $`\sigma 1`$.

###### Proof.

Consider the function $`\mathrm{\Psi }_2(U,V,\mathrm{\Omega })`$, defined by

$$\mathrm{\Psi }_2:=U(89V)+7V6.$$
(41)
The surface $`\mathrm{\Psi }_2=0`$ coincides with the regular orbit from $`L_2`$ to $`B_1`$ for $`n_0=5`$ when projected onto the $`\mathrm{\Omega }=0`$ plane, cf. (19). Taking the derivative of this function and evaluating it on the surface $`\mathrm{\Psi }_2=0`$ yields

$$\frac{d\mathrm{\Psi }_2}{d\lambda }=V(1V)\left((56\mathrm{{\rm Y}})(1U)^2+2\sigma ((4U+U\sigma 1)(34U)3\mathrm{{\rm Y}}(1U)(1+U\sigma ))\right).$$
(42)
When $`\mathrm{\Gamma }_N\frac{6}{5}`$ ($`\mathrm{{\rm Y}}\frac{5}{6}`$) and $`\sigma 1`$, then $`d\mathrm{\Psi }_2/d\lambda <0`$. Since in addition the unstable subspaces of the fixed points $`L_2`$ lie in the region $`\mathrm{\Psi }_2<0`$ of the state space, it follows that a regular orbit can never leave this region. The only attractor in this part of the state space is the fixed point $`B_4`$, cf. Theorem 4.1. Since $`B_4`$ gives rise to perfect fluid solutions with infinite masses and radii the claim of the theorem is established.

###### Remark.

The assumption $`\sigma 1`$ is just the dominant energy condition, and is satisfied, e.g., for causal relativistic polytropes. In the theorem, the condition $`\sigma 1`$ is a sufficient but not necessary; the statement is valid for much larger $`\sigma `$, but not for arbitrarily large values.

Let us define $`\sigma _l:=\left(4+\sqrt{16+2(5\mathrm{\Gamma }_N6)}\right)/6`$.

###### Theorem 6.2.

All regular solutions with $`\mathrm{\Omega }_c(0,\mathrm{\Omega }_{\mathrm{max}})`$ have finite masses and radii if $`\mathrm{\Gamma }_N>\frac{6}{5}`$ and if $`\sigma <\sigma _l`$ on $`(0,\mathrm{\Omega }_{\mathrm{max}})`$.

###### Proof.

Consider again the function $`\mathrm{\Psi }_2`$ as defined in (41). If $`\mathrm{\Gamma }_N>6/5`$ (i.e., $`\mathrm{{\rm Y}}<5/6`$), then $`d\mathrm{\Psi }_2/d\lambda >0`$ is satisfied, provided that $`\sigma <\sigma _l`$. Under the same conditions, the unstable subspace of a fixed point on $`L_2`$ lies in the region $`\mathrm{\Psi }_2>0`$. It follows that a regular orbit, which originates from a fixed point $`(3/4,0,\mathrm{\Omega }_c)`$ on $`L_2`$, is confined to the region $`\mathrm{\Psi }_2>0`$ of the state space, if $`\sigma (\mathrm{\Omega })<\sigma _l`$ for all $`\mathrm{\Omega }\mathrm{\Omega }_c`$. The only attractor in the region $`\mathrm{\Psi }_2>0`$ is the fixed point $`B_2`$, and since $`B_2`$ generates solutions with finite masses and radii, the theorem is established.

###### Remark.

In  it is shown that if $`6/5<\mathrm{\Gamma }_0<2`$ then the solution possesses a finite radius provided that $`\sigma `$ is sufficiently small (Theorem 4, p.994). In  and  theorems related to Theorem 6.2 are proved with completely different methods. In the above theorem, the lower limit for $`\sigma `$ is not a particularly good lower bound for “most” equations of state. Indeed, as follows from the next theorem, when $`n_03`$, then all solutions have finite masses and radii for all $`\mathrm{\Omega }_c[0,1]`$ and thus all values of $`\sigma `$ are allowed.

###### Theorem 6.3.

(Finiteness of perfect fluid solutions). All regular and non-regular perfect fluid solutions have finite radii and masses if $`n_03`$ (i.e., $`\mathrm{\Gamma }_03/4`$).

###### Proof.

This theorem follows from Theorem 4.2, which shows that all solutions end at $`B_2`$ if $`n_03`$ (at $`L_4`$ if $`n_0=0`$), and from that $`B_2`$ ($`L_4`$) is associated with solutions with finite masses and radii.

###### Remark.

(General validity). The last theorem holds irrespective of the asymptotic behavior of the equation of state at high pressures. Also for the first two theorems the asymptotic high pressure regime is irrelevant if one restricts the attention to solutions with finite central pressures. This is because the proofs in such cases do not rely on the inclusion of the boundary $`\mathrm{\Omega }=1`$ in the state space. This is in contrast to the next theorem which makes use of the $`\mathrm{\Omega }=1`$ subset and relies on our assumptions about differentiability of the equation of state when $`\mathrm{\Omega }1`$ ($`p\mathrm{}`$).

###### Theorem 6.4.

(Spiral structure of the $`(M,R)`$-diagram). Let $`0<n_03`$.212121More generally, we can consider equations of state with arbitrary $`n_0`$ as long as the orbit that originates from $`T_3`$ ends at $`B_2`$ (or $`L_4`$ if $`n_0=0`$).  For sufficiently high central pressures the mass-radius diagram exhibits a spiral structure, i.e., $`(R(p_c),M(p_c))`$ is given by

$$\left(\begin{array}{c}R(p_c)\\ M(p_c)\end{array}\right)=\left(\begin{array}{c}R_O\\ M_O\end{array}\right)+\left(\frac{1}{p_c}\right)^{\tau _1}𝒥(\tau _2\mathrm{log}p_c)b+o(\left(\frac{1}{p_c}\right)^{\tau _1}),$$
(43)
where $`R_O`$ and $`M_O`$ are constants, $``$ is a matrix with positive determinant, and $`b`$ a non-zero vector. The matrix $`𝒥(\phi )\mathrm{SO}(2)`$ describes a (positive) rotation by an angle $`\phi `$, and the constants $`\tau _1`$ and $`\tau _2`$ are given by

$$\tau _1=\frac{3\gamma _12}{4\gamma _1},\tau _2=\frac{1}{4\gamma _1}\sqrt{c},$$
(44)
where $`\sqrt{c}=\sqrt{36+44\gamma _1\gamma _1^{\mathrm{\hspace{0.25em}2}}}`$ see Table 1.

Sketch of proof. The basic observation is that the regular orbit on $`\mathrm{\Omega }=1`$ forms a spiral as it converges to $`T_3`$, and that this spiral subsequently reflects itself as a spiral in the $`(R,M)`$-diagram.

Choose $`\epsilon `$ small and set $`\omega _\epsilon =1\epsilon `$. The dynamical system asymptotically decouples as $`\mathrm{\Omega }1`$ (recall from (12) and (13) that we can always choose $`\omega `$ in such a way that $`\frac{d\mathrm{{\rm Y}}}{d\mathrm{\Omega }}=O((1\mathrm{\Omega })^k)`$ and $`\frac{d\sigma }{d\mathrm{\Omega }}=O((1\mathrm{\Omega })^k)`$ as $`\mathrm{\Omega }1`$), which enables us to approximately solve the system on $`[0,1]^2\times [\mathrm{\Omega }_\epsilon ,\mathrm{\Omega }]`$: let $`U_\mathrm{L}(\lambda )`$, $`V_\mathrm{L}(\lambda )`$ denote the regular solution of the system (8a) and (8b) with $`\mathrm{\Gamma }_N=1`$, $`\sigma =\sigma _1=\gamma _11`$. Then $`(U_\mathrm{L}(\lambda ),V_\mathrm{L}(\lambda ),\mathrm{\Omega }(\lambda ))`$ is an approximate solution of (8) when $`\frac{d\mathrm{\Omega }}{d\lambda }=f_1(1\mathrm{\Omega })H_\mathrm{L}`$, where $`H_\mathrm{L}=H(U_\mathrm{L},V_\mathrm{L},1)`$.

For large $`\lambda `$ the orbit $`(U_\mathrm{L}(\lambda ),V_\mathrm{L}(\lambda ))`$ has the form of a spiral, i.e.,

$$\left(\begin{array}{c}U_\mathrm{L}(\lambda )\\ V_\mathrm{L}(\lambda )\end{array}\right)=\left(\begin{array}{c}U_{T_3}\\ V_{T_3}\end{array}\right)+\mathrm{exp}(\delta _1\lambda )^{}𝒥(\delta _2\lambda )b^{},$$
(45)
where $`U_{T_3},V_{T_3}`$ are the coordinates of the fixed point $`T_3`$; $`^{}`$ is a matrix with $`det^{}>0`$, and $`b^{}`$ a non-zero vector; $`\delta _1=\frac{\gamma _1}{4(\gamma _1^2+2)}(3\gamma _12)`$ and $`\delta _2=\frac{\gamma _1}{4(\gamma _1^2+2)}\sqrt{c}`$. The constants $`\delta _1`$ and $`\delta _2`$ are the real and imaginary part of the complex eigenvalue that is associated with $`T_3`$, cf. with Table 1.

Inserting (45) into the equation for $`\mathrm{\Omega }`$ we obtain the approximation

$$1\mathrm{\Omega }(\lambda )=(1\mathrm{\Omega }_c)const\mathrm{exp}\left(f_1H_{T_3}(\lambda \stackrel{~}{\lambda })\right),$$
(46)
where $`H_{T_3}=H(U_{T_3},V_{T_3},1)=\gamma _1^2/(2+\gamma _1^2)`$, and $`\stackrel{~}{\lambda }`$ and $`const`$ are independent of $`\mathrm{\Omega }_c`$ (the value of $`\mathrm{\Omega }(\lambda )`$ as $`\lambda \mathrm{}`$).

Consider the orbit that originates from $`T_3`$ into the state space: it intersects the plane $`\mathrm{\Omega }=\mathrm{\Omega }_\epsilon `$ in a point $`I`$ with approximately the same $`(U,V)`$-coordinates as $`T_3`$. A regular solution with $`\mathrm{\Omega }_c`$ sufficiently close to $`1`$ can be described by (45) and (46). It intersects $`\mathrm{\Omega }=\mathrm{\Omega }_\epsilon `$ at $`\lambda =\lambda _c`$, where

$$\lambda _c=\frac{1}{aH_{T_3}}\mathrm{log}\frac{1}{1\mathrm{\Omega }_c}+const=\frac{1}{H_{T_3}}\mathrm{log}p_c+const.$$
(47)
For the components $`(U(\lambda _c),V(\lambda _c))`$ of this regular solution we have $`(U(\lambda _c),V(\lambda _c))=(U_{T_3},V_{T_3})+(\delta U_0,\delta V_0)`$ with

$$\left(\begin{array}{c}\delta U_0(\eta _c)\\ \delta V_0(\eta _c)\end{array}\right)=\left(\frac{1}{p_c}\right)^{\tau _1}^{\prime \prime }𝒥(\tau _2\mathrm{log}p_c)b^{\prime \prime },$$
(48)
where $`^{\prime \prime }`$ ($`det^{\prime \prime }>0`$) and $`b^{\prime \prime }`$ have absorbed the constants. From (48) we see that the set of regular solutions intersects $`\mathrm{\Omega }=\mathrm{\Omega }_\epsilon `$ in a spiral in a neighborhood $`U_I`$ of $`I`$ in this plane.

The orbit that originates from $`T_3`$ and passes through $`IU_I`$ eventually intersects the mass-radius cylinder (cf. (40)) at some point $`OU_O\mathrm{cylinder}`$, where the cylindrical coordinates of $`O`$ determine the mass $`M_O`$ and the radius $`R_O`$ of the solution (cf. (39)). The flow of the dynamical system induces a diffeomorphism $`U_IU_O`$, so that the spiral (48) in $`U_I`$ is mapped to a distorted (positively oriented) spiral in $`U_O`$. Since the cylindrical coordinates $`(\varphi ,h)`$ in the small neighborhood $`U_O`$ are related to $`(R,M)`$ by a positively oriented linear map, the spiral (48), with different $``$ and $`b`$, appears also in the $`(R,M)`$ diagram.

###### Remark.

Note that the theorem does not rely on the behavior of the equation of state in the intermediate non-asymptotic regimes. Hence the theorem describes a universal phenomenon and connects it with the self-similar solution that corresponds to $`T_3`$. A similar theorem has been proved in  using quite different methods. However, an advantage with the present approach, apart from less restrictive assumptions, brevity of proof, and clarity, is that it visually shows the importance of key solutions like the self-similar one associated with $`T_3`$.

## 7 Examples

In this section the dynamical systems approach to relativistic stellar models will be illustrated by several examples. The presented results also demonstrate the usefulness of the new framework in numerical computations.

### 7.1 Relativistic polytropes

As discussed in Sec. 3, a relativistic polytrope with index $`\mathrm{\Gamma }`$ can be described by the dynamical system (8) with $`\sigma =(\mathrm{\Gamma }1)\mathrm{\Omega }`$ and $`\mathrm{{\rm Y}}=(1+(\mathrm{\Gamma }1)\mathrm{\Omega })/\mathrm{\Gamma }`$. The numerical integration of this system is straight forward, when we use the local analysis near the fixed points. The radius $`R`$ and the total mass $`M`$ of a solution can be obtained by various methods, in particular one can use equations (35d) and (35e), or one can include the auxiliary equations (35a) in the numerical integration of the dynamical system; Eq. (39) provides a different approach. Mass-radius diagrams for the regular solutions associated with different values of $`\mathrm{\Gamma }`$ are shown in Fig. 5.

The theorems of Sec. 6 guarantee that all regular solutions possess finite $`R`$ and $`M`$ for $`n_03`$ ($`\mathrm{\Gamma }4/3`$). For $`5>n_0>3`$ ($`\mathrm{\Gamma }<4/3`$), regular solutions with small central pressures must be finite fluid bodies. The numerical investigation improves these results: all regular solutions associated with $`\mathrm{\Gamma }1.2971`$ ($`n_03.366`$) possess finite $`R`$ and $`M`$. A similar analysis was performed in  for purely polytropic equations of state: for a polytropic index $`n3.339`$ only finite fluid bodies appear. For an extensive discussion of the methods involved in establishing these results we refer the reader to that paper.

Letting $`\mathrm{\Gamma }1`$ for relativistic polytropes implies that $`\sigma 0`$, and hence the dynamical system (8) approaches the Newtonian system, cf. Appendix A. Moreover, the equation for $`\mathrm{\Omega }`$ decouples from the $`U`$\- and $`V`$-equations, and the flows on the $`\mathrm{\Omega }=0`$ and the $`\mathrm{\Omega }=1`$ subset coincide in the limit. It follows that the regular surface is a surface perpendicular to the $`\mathrm{\Omega }=const`$ plane; its projection to $`\mathrm{\Omega }=const`$ is depicted in Fig. 6(c). Note also that the limit $`\mathrm{\Gamma }1`$ reveals a relationship between the Newtonian self-similar solution associated with $`B_4`$ (discussed, e.g., in ) and the relativistic self-similar solution associated with $`T_3`$.

Regular solutions associated with an equation of state characterized by an arbitrary adiabatic index $`\mathrm{\Gamma }(p)`$ close to $`1`$ lie close to the regular ”soft-limit-surface”.222222Consider an arbitrary one-parameter class of equations of state characterized by a parameter $`\mu >0`$, where $`\mathrm{\Gamma }_0`$ and $`\gamma _1`$ converge to $`1`$ as $`\mu 0`$. Then the associated regular surface approaches the soft-limit-surface as $`\mu 0`$. Hence, deviations from the soft-limit-surface measure both deviations of the equation of state from the soft limit and relativistic effects. Stiffer models are located ”further out” compared to the soft-limit-surface and increasing stiffness implies that the $`U=3/4`$ surface is approached, and it is on this “stiff-limit-surface” the incompressible perfect fluid solutions reside, as discussed previously. Note also that the Buchdahl surface, depicted in Fig. 3, for increasingly soft equations of state moves toward the $`V=1`$ surface, thereby becoming increasingly less restrictive. This illustrates that the Buchdahl inequality is a purely relativistic effect, since the surface disappears from the interior state space in the soft, and thereby Newtonian, limit.

### 7.2 The ideal neutron gas

The equation of state of the ideal neutron gas provides a simple model for neutron star matter; it is given in implicit form in equation (18). As discussed in Sec. 3, the ideal neutron gas can be naturally described in the dynamical systems framework, and is easy to handle numerically.

In Fig. 7 regular orbits for various values of $`\mathrm{\Omega }_c`$ and the regular orbits on $`\mathrm{\Omega }=0`$ and $`\mathrm{\Omega }=1`$ are depicted.

### 7.3 Composite equations of state

Consider first a stiff perfect fluid characterized by the equation of state $`p=(\rho \rho _{})c^2`$, which can be rewritten as

$$k\rho c^2=1+kp,$$
(49)
where $`k=(\rho _{}c^2)^1`$. If we choose $`\omega =kp`$, we obtain $`f=1`$, $`\mathrm{{\rm Y}}=\mathrm{\Omega }`$, and $`\sigma =\mathrm{\Omega }`$, i.e., the dynamical system (8) becomes very simple. In this context, it is of interest to note that the solution that corresponds to the orbit from $`T_3`$ to $`L_4`$ (which is a solution with finite radius and mass but infinite central pressure) is explicitly known: it is a special case of a Tolman solution . In the present formalism, the corresponding orbit is given by $`U=3/(5+\mathrm{\Omega })`$ and $`V=(2+\mathrm{\Omega })/(2+3\mathrm{\Omega }+4\mathrm{\Omega }^2)`$. (For interesting features of this solution, see .)

As an example of a composite equation of state, let us continuously join a stiff equation of state $`p=(\rho \rho _{})c^2`$ (when $`p>p_j`$) to a relativistic polytrope (when $`p<p_j`$).232323If an equation of state is only known for low pressure (as is true in reality), then mass estimates for stellar models can be obtained by extending the equation of state as a stiff fluid, see e.g., . Such an equation of state possesses a kink at $`p_j`$.

As outlined in Sec. 3, we work with two state spaces and dynamical systems, one pertaining to the relativistic polytrope and one representing the stiff fluid. The jump at $`p_j`$ reflects itself in a map between the two state spaces. Since $`\rho `$ is continuous, the values of $`U,V`$ are continuous under the map, however, there is a jump in $`\mathrm{\Omega }`$. Since $`\rho _{\mathrm{relpol}}(p_j)=\rho _{\mathrm{stiff}}(p_j)`$ we have

$$\frac{1}{\mathrm{\Gamma }1}\left(\frac{(kp_j)^{1/\mathrm{\Gamma }}}{k}+p_j\right)=\rho _{}c^2+p_j.$$
(50)
Dividing by $`p_j`$ and using that $`\omega _{\mathrm{relpol}}(p)=(kp)^{(\mathrm{\Gamma }1)/\mathrm{\Gamma }}`$ and $`\omega _{\mathrm{stiff}}(p)=(\rho _{}c^2)^1p`$, we obtain that

$$\mathrm{\Omega }_{\mathrm{relpol}}=\frac{\mathrm{\Omega }_{\mathrm{stiff}}}{\mathrm{\Gamma }1}$$
(51)
at the jump. If $`\mathrm{\Gamma }=2`$, i.e., if the relativistic polytrope is asymptotically stiff, then there is no jump at all. A typical solution with a jump is illustrated in Fig. 8.

Note that, in general, if one starts with a regular solution in the first state space, then this solution has to be matched with a non-regular one associated with the other (extended) equation of state in the second state space.

## 8 Concluding remarks

In this paper we have derived a dynamical systems formulation for the study of spherically symmetric relativistic stellar models. The method of ”homology invariants” known from the theory of Newtonian polytropes has been generalized both to a general relativistic context and to the broad class of asymptotically polytropic/linear equations of state. For previous related work see  (Newtonian asymptotic polytropes) and  (polytropes in GR).

The present formulation has turned out to be advantageous in many respects: (i) it provides a clear visual representation of the solution spaces associated with broad classes of equations of state; in particular it is revealed how the global qualitative properties of the solution space are influenced by the equation of state; (ii) the formulation makes the theory of dynamical systems available, which makes it possible to prove a number of theorems and describe the qualitative behavior of solutions; and (iii) the framework is particularly suited for numerical computations, since the numerics is supported by the local dynamical systems analysis.

The idea of using dimensionless variables and exploiting asymptotic symmetries and properties to derive a dynamical systems formulation is a quite versatile one. As another example, in a future paper, we will show how one also can treat a collisionless gas. The ”dynamical systems approach” can even be applied to problems, in general relativity and other areas, when no symmetries exist at all , as will be shown in another set of papers. Thus the ideas in the present work should be seen in a quite broad context, and it is likely that there are many problems in very different areas that could benefit from the type of ideas and the approach presented in this paper.

### Acknowledgements

The authors would like to thank Alan Rendall and Robert Beig (J.M.H.) for helpful discussions. This work was supported by the DOC-program of the Austrian Academy of Sciences (J.M.H.) and the Swedish Research Council (C.U.).

## Appendix A Appendix: Comparison with the Newtonian theory

The Newtonian equations of hydrostatic equilibrium are obtained from (2) and (3) by letting $`c\mathrm{}`$, however, note that in the Newtonian case $`\rho `$ is interpreted as the rest-mass density and $`c^2\mathrm{\Phi }\mathrm{\Phi }_N`$, the Newtonian potential. Consistently, since $`\sigma =p/(\rho c^2)`$, in the dynamical system (8) we must set $`\sigma 0`$. Note also that $`v`$ reduces to $`v=G\rho m/(rp)`$. Since $`\sigma 0`$ for $`p0`$, the relativistic system (8) and its Newtonian counterpart coincide on $`\mathrm{\Omega }=0`$.

In the Newtonian case we can treat equations of state $`\rho (p)`$ with less restrictive asymptotic behavior in the regime $`p\mathrm{}`$. This is because it was the existence of $`\sigma (\mathrm{\Omega })`$ on the right hand side of the dynamical system that forced the assumption of asymptotic linearity upon us in the relativistic case. In particular we can naturally include asymptotically polytropic behavior as $`p\mathrm{}`$ in the Newtonian case (cf. ).

In Newtonian gravity there exist two independent dimensional scales, space and time, while in relativity, space and time are connected by the speed of light $`c`$, so that only one single dimensional scale remains. A Newtonian equation of state can be written in the implicit form

$$k_1p=\varphi (\omega ),k_2\rho =\chi (\omega ),$$
(52)
where $`\omega `$ is a dimensionless variable and $`k_1`$ and $`k_2`$ are constants carrying the dimension $`[p]^1`$ and $`[\rho ]^1`$ respectively; $`\varphi `$ and $`\chi `$ are monotone functions in $`\omega `$ containing any number of dimensionless parameters. An explicit representation of (52) is $`k_2\rho =\psi (k_1p)`$.

Since our dynamical system is expressed in terms of purely dimensionless quantities, the dimensional parameters have to drop out. Consequently, a single dynamical system must be capable of describing the entire two-parameter set of equations of state (52) parameterized by $`k_1`$ and $`k_2`$. Indeed, the relevant function entering the dynamical system in the Newtonian case is $`\mathrm{{\rm Y}}(\omega )`$ given by $`\mathrm{{\rm Y}}(\omega )=(d\mathrm{ln}\chi /d\mathrm{ln}\omega )(d\mathrm{ln}\varphi /d\mathrm{ln}\omega )^1`$. Evidently, this function does not depend on $`k_1`$ and $`k_2`$. The same clearly holds for $`f(\omega )=d\mathrm{ln}\varphi /d\mathrm{ln}\omega `$.

###### Example.

As an illustrative example, consider the equation of state $`\rho =C_1p^{\mathrm{{\rm Y}}_0}+C_2p^{\mathrm{{\rm Y}}_1}`$. By rescaling the equation and the two parameters one can write the equation of state as $`k_2\rho =(k_1p)^{\mathrm{{\rm Y}}_0}+(k_1p)^{\mathrm{{\rm Y}}_1}`$, and thus the equation has been brought to the desired form. The choice $`\omega =(k_1p)^a`$ ensures that $`\mathrm{{\rm Y}}(\omega )`$ and therefore also $`\mathrm{{\rm Y}}(\mathrm{\Omega })`$ are independent of $`k_1`$ and $`k_2`$, and interpolate monotonically between $`\mathrm{{\rm Y}}_0`$ and $`\mathrm{{\rm Y}}_1`$. In analogy to the detailed discussion in Sec. 3, the constant $`a`$ must be chosen sufficiently small so that a $`𝒞^1`$-differentiable dynamical system is obtained.

###### Remark.

(Consistent translatory invariance). As seen previously, it is the freedom to choose $`\omega (p)`$ that allows one to treat a one-parameter class of equations of state simultaneously by one specified dynamical system in the relativistic case. The freedom we were able to exploit in the Newtonian case to cover even two-parameter classes of equations of state reflects itself in the consistent translatory invariance of the Newtonian dynamical system. Namely, if $`𝐱(\xi )`$ ($`𝐱(\xi ):=(u(\xi ),v(\xi ),\omega (\xi ))`$ (recall that $`\xi =\mathrm{ln}r`$) is a solution of the Newtonian dynamical system, then so is $`𝐱_0:=𝐱(\xi \xi _0)`$, and moreover, the translated solution gives rise to a perfect fluid solution associated with an equation of state with the same $`\mathrm{{\rm Y}}(p)`$. In contrast, in the relativistic case this consistency is broken: if $`𝐱(\xi )`$ is a solution of the dynamical system (6), then so is $`𝐱_0(\xi )=𝐱(\xi \xi _0)`$. However, whereas the solution $`𝐱(\xi )`$ can be consistently interpreted as a relativistic perfect fluid solution, associated with an equation of state $`\rho (p)`$ (determined by $`\mathrm{{\rm Y}}`$ and $`\sigma `$), this is not the case with the translated solution. Although the translated solution satisfies differential equations associated with $`\rho (p)`$, the initial data are not consistent with this equation of state. Accordingly, only one particular solution on every orbit of (6), can be interpreted as a relativistic perfect fluid solution. However, in the compactified dynamical system (8) this ‘defect’ is remedied by the appearance of the freedom in the new independent variable $`\lambda `$ (this ‘cure’ could also have been implemented in the $`u,v,\omega `$-formulation by defining $`\xi `$ through $`dr/d\xi =r`$, i.e, by letting $`\xi =\mathrm{ln}r+\xi _0`$, where $`\xi _0`$ is an arbitrary constant, instead of setting $`\xi =\mathrm{ln}r`$, however, the direct relation between $`\xi `$ and $`r`$ is sometimes useful). The reason behind the difference in the Newtonian and relativistic cases is the appearance of $`\sigma `$ in the relativistic dynamical system; $`\sigma `$ uniquely determines the equation of state while $`\mathrm{{\rm Y}}`$ only does so up to a proportionality constant.

To conclude the remarks on the Newtonian case, we consider the Newtonian counterpart of the mass-radius formulas (39). These equations simplify considerably in the Newtonian case. Starting from the two-parameter family (52) we obtain for the surface potential $`\mathrm{\Phi }_N(R)=GM/R`$ the expression $`GM/R=k_1^1k_2B/s`$, and for $`R^2`$ and $`M^2`$,

$$R^2=\frac{k_2^2}{4\pi s^2k_1G}AB^{1n_0},M^2=\frac{k_2^4}{4\pi s^4k_1^3G^3}AB^{3n_0}.$$
(53)