The Problem and its Solution

We will present, in the case of a spherically symmetric distribution of gaseous fluid satisfying a polytropic equation of state, the complete static solution, over all the three-dimensional space, of the Einstein field equations of General Relativity. In this work we will use the time-like signature $(+,-,-,-)$, following [#!DiracGravity!#]. We will start from the same differential system already described in [#!LiquidShells!#], which we will succinctly review here. Just as in [#!LiquidShells!#], the solution will be given in terms of the coefficients of the metric, for an invariant interval given in terms of the Schwarzschild coordinates $(t,r,\theta,\phi)$ by

  $\displaystyle
ds^{2}
=
\,{\rm e}^{2\nu(r)}
c^{2}dt^{2}
-
\,{\rm e}^{2\lam...
...r^{2}
-
r^{2}
\left[
d\theta^{2}
+
\sin^{2}(\theta)
d\phi^{2}
\right],
$ (1)

where $\exp[\nu(r)]$ and $\exp[\lambda(r)]$ are two positive functions of only $r$. As was shown in detail in [#!LiquidShells!#], under these conditions the independent components of the Einstein field equations and the Bianchi consistency condition are equivalent to the set of three first-order differential equations


$\displaystyle \left\{\rule{0em}{2.5ex}1-2\left[r\lambda'(r)\right]\right\}
\,{\rm e}^{-2\lambda(r)}$ $\textstyle =$ $\displaystyle 1-\kappa r^{2}\rho(r),$ (2)
$\displaystyle \left\{\rule{0em}{2.5ex}1+2\left[r\nu'(r)\right]\right\}
\,{\rm e}^{-2\lambda(r)}$ $\textstyle =$ $\displaystyle 1+\kappa r^{2}P(r),$ (3)
$\displaystyle \left[\rho(r)+P(r)\right]
\left[r\nu'(r)\right]$ $\textstyle =$ $\displaystyle -\left[rP'(r)\right],$ (4)

where the primes indicate differentiation with respect to $r$, $\rho(r)$ is the energy density of the matter, $P(r)$ is its isotropic pressure, and where we have the constant $\kappa=8\pi G/c^{4}$, in which $G$ is the universal gravitational constant and $c$ is the speed of light. Note that all the derivatives are written as what we will call homogeneous derivatives, that is, the product of the derivative by a single power of $r$. In our case here the matter distribution will be characterized by four parameters, the two parameters defining the polytropic equation of state, the total asymptotic gravitational mass $M$, associated to the Schwarzschild radius $r_{M}$, and a parameter associated to the value of the energy density $\rho(r)$ at its point of maximum. We will assume that the gas satisfies the polytropic equation of state

  $\displaystyle
P(r)
=
K\left[\rho(r)\right]^{1+1/n},
$ (5)

over the whole three-dimensional space, involving a positive real constant $K$ and the integer or half-integer $n>1$, which we assume to be strictly larger that one. In principle $n$ could be any real number larger than one, and we assume that it is either an integer of a half integer just for simplicity, since this seems to cover all cases of interest. At this point we will introduce an auxiliary function, also just for simplicity, since it will appear repeatedly in all that follows,

  $\displaystyle
F(r)
=
K\left[\rho(r)\right]^{1/n},
$ (6)

which is a dimensionless function, so that the polytropic equation of state reduces to

  $\displaystyle
P(r)
=
F(r)\rho(r).
$ (7)

Note that from its definition we immediately have for the derivative of $F(r)$,

  $\displaystyle
\left[rF'(r)\right]
=
\frac{1}{n}\,
\frac{F(r)}{\rho(r)}\,
\left[r\rho'(r)\right].
$ (8)

Given this, our system of differential equations shown in Equations (2)-(4) can now be written as


$\displaystyle \left\{\rule{0em}{2.5ex}1-2\left[r\lambda'(r)\right]\right\}
\,{\rm e}^{-2\lambda(r)}$ $\textstyle =$ $\displaystyle 1-\kappa r^{2}\rho(r),$ (9)
$\displaystyle \left\{\rule{0em}{2.5ex}1+2\left[r\nu'(r)\right]\right\}
\,{\rm e}^{-2\lambda(r)}$ $\textstyle =$ $\displaystyle 1+\kappa r^{2}F(r)\rho(r),$ (10)
$\displaystyle \left[1+F(r)\right]\rho(r)
\left[r\nu'(r)\right]$ $\textstyle =$ $\displaystyle -\,
\frac{n+1}{n}\,
F(r)\left[r\rho'(r)\right],$ (11)

still in terms of homogeneous derivatives. Our problem is therefore that of finding three functions, $\rho(r)$, $\lambda(r)$ and $\nu(r)$, that solve these equations and that satisfy the correct boundary conditions at asymptotic radial infinity. We will start our analysis by partially solving some of the equations by analytic means, in order to write all relevant quantities in terms of a single function. In order to do this we first change variables from the dimensionless function $\lambda(r)$ to the equally dimensionless function $\beta(r)$, which is defined to be such that

  $\displaystyle
\,{\rm e}^{2\lambda(r)}
=
\frac{r}{r-r_{M}\beta(r)},
$ (12)

which then implies that we have for the corresponding homogeneous derivatives

  $\displaystyle
2\left[r\lambda'(r)\right]
=
-\,
\frac{r_{M}\beta(r)-r_{M}\left[r\beta'(r)\right]}{r-r_{M}\beta(r)}.
$ (13)

The function $\beta(r)$ is analogous to the function $u(r)$ found in Equation (2.9) of the paper by Tooper [#!tooper!#], but it is used here in a completely different context. Note that the asymptotic boundary condition on $\lambda(r)$, that it must behave as the exterior Schwarzschild solution for sufficiently large $r$, translates here as $\beta(r)=1$ under that same condition. Substituting these expressions in the component field equation shown in Equation (9) a very simple relation giving the derivative of $\beta(r)$ in terms of $\rho(r)$ results,

  $\displaystyle
\beta'(r)
=
\frac{\kappa r^{2}\rho(r)}{r_{M}}.
$ (14)

Therefore, the energy density $\rho(r)$ is given in terms of the derivative of $\beta(r)$, and wherever $\rho(r)=0$, characterizing a region where there is a vacuum, we have that $\beta(r)$ is a constant. Since $F(r)$ and $P(r)$ are both given in terms of $\rho(r)$, and since $\lambda(r)$ is given in terms of $\beta(r)$, it follows at this point that, given a function $\beta(r)$, the functions $\rho(r)$, $P(r)$, $F(r)$ and $\lambda(r)$ are all determined. The only function that has yet to be determined in terms of $\beta(r)$ is $\nu(r)$. We can obtain $\nu(r)$ in terms of $F(r)$, and therefore of $\beta(r)$, using the consistency condition in Equation (11), which with the use of Equation (8) can be written as

  $\displaystyle
\nu'(r)
=
-(n+1)\,
\frac{F'(r)}{1+F(r)}.
$ (15)

This can now be integrated from $r_{2}$ to $r$, and recalling that we have that $F(r_{2})=0$, because for $n>0$ we have the boundary condition $\rho(r_{2})=0$, we get

  $\displaystyle
\nu(r)
-
\nu(r_{2})
=
-(n+1)\ln[1+F(r)],
$ (16)

which determines $\nu(r)$ in terms of $F(r)$, up to the integration constant $\nu(r_{2})$, that is to be obtained from the asymptotic boundary conditions at radial infinity, which in the case of $\nu(r)$ is simply $\nu(\infty)=0$. Therefore, given $\beta(r)$, this determines $\nu(r)$ in terms of it, through $F(r)$. Note, for future use, that the fact that we also have that $F(r_{1})=0$, because for $n>0$ we have the boundary condition $\rho(r_{1})=0$, implies that we always have that $\nu(r_{1})=\nu(r_{2})$. We see therefore that the determination of the function $\beta(r)$ leads with no further difficulty to the determination of all the functions that describe both the matter and the geometry of the system,


$\displaystyle \rho(r)$ $\textstyle =$ $\displaystyle \frac{r_{M}\beta'(r)}{\kappa r^{2}},$ (17)
$\displaystyle P(r)$ $\textstyle =$ $\displaystyle K\left[\frac{r_{M}\beta'(r)}{\kappa r^{2}}\right]^{1+1/n},$ (18)
$\displaystyle F(r)$ $\textstyle =$ $\displaystyle K\left[\frac{r_{M}\beta'(r)}{\kappa r^{2}}\right]^{1/n},$ (19)
$\displaystyle \lambda(r)$ $\textstyle =$ $\displaystyle \frac{1}{2}\,\ln\!\left[\frac{r}{r-r_{M}\beta(r)}\right],$ (20)
$\displaystyle \nu(r)$ $\textstyle =$ $\displaystyle \nu(r_{2})-(n+1)\ln[1+F(r)].$ (21)

The free parameters of the system are $K$, $n$ and $M$, all of which describe the nature and state of the matter, and the value of $\beta'(r)$ at its point of maximum, which can also be seen as related to the matter, since it determines the general scale of the matter energy density, as can be seen from Equation (17).

Both for the subsequent analysis and for the numerical approach, it is convenient to transform variables at this point, in order to write everything in terms of dimensionless variables and functions. In order to do this we must now introduce an arbitrary radial reference position $r_{0}>0$. For now the value of this parameter remains completely arbitrary, other than that it must be strictly positive, and has no particular physical meaning. It is only a mathematical device that allows us to define a dimensionless radial variable and a dimensionless parameter associated to the mass $M$ by


$\displaystyle \xi$ $\textstyle =$ $\displaystyle \frac{r}{r_{0}},$ (22)
$\displaystyle \xi_{M}$ $\textstyle =$ $\displaystyle \frac{r_{M}}{r_{0}},$ (23)

as well as to define the dimensionless function of $\xi$, to assume the role of $\beta(r)$,

  $\displaystyle
\gamma(\xi)
=
\xi_{M}\beta(r).
$ (24)

Note that the asymptotic condition that $\beta(r)=1$ for sufficiently large $r$ translates here as the condition that $\gamma(\xi)=\xi_{M}$ for sufficiently large $\xi$. It is important to observe that under the change of variables from $r$ to $\xi$ the homogeneous derivatives transform in a very simple way, for example in the case of $\beta(r)$,

  $\displaystyle
r\,
\frac{d\beta(r)}{dr}
=
\xi\,
\frac{d\beta(\xi)}{d\xi}.
$ (25)

We will also introduce at this point a notation for the derivative of $\gamma (\xi )$, which will be useful in order to deal with the second-order differential equation for $\gamma (\xi )$ that we will arrive at,

  $\displaystyle
\pi(\xi)
=
\gamma'(\xi),
$ (26)

where in this case the prime denotes differentiation with respect to $\xi$. We adopt the convention that whenever there is a prime indicating a derivative, it is to be taken with respect to the explicitly indicated variable of the function. From now on the problem will be formulated in terms of the two functions $\gamma (\xi )$ and $\pi (\xi )$, that will be treated as independent variables. In terms of these new variables we have for the system of differential equations shown in Equations (9)-(11), where from now on the auxiliary function $F(r)$ will be written as $F(\xi,\pi)$,


$\displaystyle 1-2\left[\xi\lambda'(\xi)\right]$ $\textstyle =$ $\displaystyle \frac{\xi[1-\pi(\xi)]}{\xi-\gamma(\xi)},$ (27)
$\displaystyle 1+2\left[\xi\nu'(\xi)\right]$ $\textstyle =$ $\displaystyle \frac{\xi[1+F(\xi,\pi)\pi(\xi)]}{\xi-\gamma(\xi)},$ (28)
$\displaystyle \left[1+F(\xi,\pi)\right]\pi(\xi)
\left[\xi\nu'(\xi)\right]$ $\textstyle =$ $\displaystyle -\,
\frac{n+1}{n}\,
F(\xi,\pi)
\left[\xi\pi'(\xi)-2\pi(\xi)\right],$ (29)

where the primes now indicate differentiation with respect to $\xi$. We therefore have for all the relevant quantities written in terms of $\xi$, $\gamma (\xi )$ and $\pi (\xi )$,


$\displaystyle \rho(\xi)$ $\textstyle =$ $\displaystyle \frac{1}{\kappa r_{0}^{2}}\,
\frac{\pi(\xi)}{\xi^{2}},$ (30)
$\displaystyle P(\xi)$ $\textstyle =$ $\displaystyle \frac{C}{\kappa r_{0}^{2}}\,
\left[\frac{\pi(\xi)}{\xi^{2}}\right]^{1+1/n},$ (31)
$\displaystyle F(\xi,\pi)$ $\textstyle =$ $\displaystyle C\left[\frac{\pi(\xi)}{\xi^{2}}\right]^{1/n},$ (32)
$\displaystyle \lambda(\xi)$ $\textstyle =$ $\displaystyle -\,
\frac{1}{2}\,\ln\!\left[\frac{\xi-\gamma(\xi)}{\xi}\right],$ (33)
$\displaystyle \nu(\xi)$ $\textstyle =$ $\displaystyle \nu(\xi_{2})-(n+1)\ln[1+F(\xi,\pi)],$ (34)

where $C=K/\left(\kappa r_{0}^{2}\right)^{1/n}$ is a dimensionless constant. We see therefore that the determination of $\gamma (\xi )$ and thus of $\pi (\xi )$ leads to the complete solution of the problem. We have therefore reduced the solution of the problem to the determination of the single function $\gamma (\xi )$. In order to determine $\gamma (\xi )$, we obtain an ordinary differential equation for it by eliminating $\nu'(\xi)$ from Equations (28) and (29). Equation (28) can be written as

  $\displaystyle
\left[\xi\nu'(\xi)\right]
=
\frac
{\gamma(\xi)+\xi F(\xi,\pi)\pi(\xi)}
{2[\xi-\gamma(\xi)]},
$ (35)

and Equation (29) can be written as

  $\displaystyle
\left[\xi\nu'(\xi)\right]
=
-\,
\frac{n+1}{n}\,
\frac{F(\xi,\pi)}{1+F(\xi,\pi)}\,
\frac{\xi\pi'(\xi)-2\pi(\xi)}{\pi(\xi)},
$ (36)

so that equating the two right-hand sides we get

  $\displaystyle
\pi'(\xi)
=
\pi(\xi)
\left\{
\frac{2}{\xi}
-
\frac{n}{n+1}...
...,\pi)}
\left[
\frac{\gamma(\xi)}{\xi}+F(\xi,\pi)\pi(\xi)
\right]
\right\}.
$ (37)

Since $\pi (\xi )$ is the derivative of $\gamma (\xi )$, this can be interpreted either as a second-order ordinary differential equation for $\gamma (\xi )$, or as one of a pair of first-order coupled ordinary differential equations determining $\gamma (\xi )$ and $\pi (\xi )$, the other equation being simply

  $\displaystyle
\gamma'(\xi)
=
\pi(\xi).
$ (38)

This second interpretation is the one we will adopt here. This pair of first-order ordinary differential equations can be used for the numerical integration of this differential system, as we will in fact do later on, in Section 4. However, before we plunge into the numerical approach, let us first show that $\gamma (\xi )$ is in fact a very simple function, as is $\beta(r)$, and that it has some properties which are, perhaps, somewhat unexpected.