Complete solution of the Einstein field equations for a spherical distribution of polytropic matter

We determine the complete solution of the Einstein field equations for the case of a spherically symmetric distribution of gaseous matter, characterized by a polytropic equation of state. We show that the field equations automatically generate two sharp boundaries for the gas, an inner one and an outer one, given by radial positions r1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{1}$$\end{document} and r2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{2}$$\end{document}, and thus define a shell of gaseous matter outside of which the energy density is exactly zero. Hence this shell is surrounded by an outer vacuum region, and surrounds an inner vacuum region. Therefore, the solution is given in three regions, one being the well-known analytical Schwarzschild exterior solution in the outer vacuum region, one being determined analytically in the inner vacuum region, and one being determined partially analytically and partially numerically, within the matter region, between the two boundary values r1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{1}$$\end{document} and r2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{2}$$\end{document} of the Schwarzschild radial coordinate r. This solution is therefore somewhat similar to the one previously found for a spherically symmetric shell of liquid fluid, and is in fact exactly the same in the cases of the inner and outer vacuum regions. The main difference is that here the boundary values r1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{1}$$\end{document} and r2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{2}$$\end{document} are not chosen arbitrarily, but are instead determined by the dynamics of the system. As was shown in the case of the liquid shell, also in this solution there is a singularity at the origin, that just as in that case does not correspond to an infinite concentration of matter, but in fact to zero matter energy density at the center. Also as in the case of the liquid shell, the spacetime within the spherical cavity is not flat, so that there is a non-trivial gravitational field there, in contrast with Newtonian gravitation. This inner gravitational field has the effect of repelling matter and energy away from the origin, thus avoiding a concentration of matter at that point.


Introduction
In a previous paper [1] we established the exact static solution of the Einstein field equations for the case of a spherically symmetric shell of liquid fluid located between two arbitrary radial positions r 1 and r 2 of the Schwarzschild system of coordinates. In this paper we will give the complete solution for a similar problem, that of a spherically symmetric distribution of a gaseous fluid that satisfies the equation of state of a polytrope. We will see that for most sets of values of the physical parameters of the system the Einstein field equations coupled with the polytropic equation of state automatically imply the existence of certain radii r 1 and r 2 where the energy density of the gas becomes exactly zero, thus giving rise to two gas-vacuum interfaces. These two values of the radial variable are not imposed by hand, but are a consequence of the equations describing the dynamics of the system. There are no geometrical free parameters, all the free parameters of the system are those describing the state and properties of the matter.
This puts us in a position, in the current problem, and in a very natural way, which is very similar to the one we had in [1], with a shell of fluid matter surrounding an internal vacuum region and surrounded by an external vacuum region. In these two vacuum regions the solutions of the field equations are known exactly, and were in fact derived and discussed in detail in [1]. Consequently, all that was said in that paper regarding the inner and outer vacuum regions is valid here without any change. However, our current problem in this paper is far less academic in nature, being much closer to the astrophysical applications. The family of solution that we find here can be considered a generalization of the particular family of solutions originally found by Tooper [2].
Results similar to the ones we present here were obtained for the case of neutron stars by Ni [3], including the automatic generation of the inner and outer matter-vacuum interfaces. However, the crucial consideration of the interface boundary conditions was missing from that analysis, thus leading to incomplete results. The discussion of the interface boundary conditions was subsequently introduced by Neslušan [4,5], thus completing the analysis of the case of the neutron stars. Just as in [1] and in the present work, the discussion of the interface boundary conditions led, also in that case, to an inner vacuum region containing a singularity at the origin and a gravitational field leading matter and energy away from that origin.
There is a connection between the purely numerical results presented in [3] and [4] and the ones we present here, which are partly analytical and partly numerical. This is so because the equation of state for a neutron star can be approximated by a polytrope under certain particular conditions. On the other hand, however, the results we present here are not limited to cold neutron stars or some other particular type of object, but can be applied as well to normal stars of any type and size, at any temperature range, such as main sequence stars, red giants and white dwarfs, including configurations with two or more layers, with a different behavior of the matter in each layer.
In regard to the solution within the region containing the polytropic matter, we will present a solution which is partially analytical and partially numerical. We will reduce, by analytic means, all the quantities and functions involved to a single real function that is the solution of a second-order ordinary differential equation. We will show analytically that this function has a very simple behavior, which can be rigorously established without any recourse to numerical means. This function, which we will denote by β(r ), is the single element of the whole system that eventually has to be determined in detail numerically, but its most important properties are analytically established beforehand.
The two radial positions r 1 and r 2 where the energy density becomes zero are soft singular points of the function β(r ), where by "soft" we mean that the function is not analytic, but also does not diverge to infinity at those points. Besides, the positions of these points are not known beforehand, since they are a consequence of the dynamics of the system. Due to all this, it is not really possible to integrate the differential equation numerically from these positions, which constitute the boundary of the region of space containing the matter, into the interior of the matter region. We will therefore have to develop an alternative way to solve this kind of differential problem numerically, since it is significantly different from a typical boundary value problem. This paper is organized as follows. In Sect. 2 we gill give the full statement of the problem and describe the resolution method; in Sect. 3 we will obtain analytically the main properties of the solutions; in Sect. 4 we will present a few examples of complete numerical solutions; in Sect. 5 we will analyze and comment on the numerical results obtained; and in Sect. 6 we will present our conclusions.

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 threedimensional space, of the Einstein field equations of General Relativity. In this work we will use the time-like signature (+, −, −, −), following [6]. We will start from the same differential system already described in [1], which we will succinctly review here. Just as in [1], 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 , θ, φ) by where exp[ν(r )] and exp[λ(r )] are two positive functions of only r . As was shown in detail in [1], 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 where the primes indicate differentiation with respect to r , ρ(r ) is the energy density of the matter, P(r ) is its isotropic pressure, and where we have the constant κ = 8π 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 ρ(r ) at its point of maximum. We will assume that the gas satisfies the polytropic equation of state over the whole three-dimensional space, involving a positive real constant K and the integer or half-integer n ≥ 1, which we assume not to be smaller than 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, which is a dimensionless function, so that the polytropic equation of state reduces to Note that from its definition we immediately have for the derivative of F(r ), Given this, our system of differential equations shown in Eqs. (2)-(4) can now be written as still in terms of homogeneous derivatives. Our problem is therefore that of finding three functions, ρ(r ), λ(r ) and ν(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 λ(r ) to the equally dimensionless function β(r ), which is defined to be such that which then implies that we have for the corresponding homogeneous derivatives The function β(r ) is analogous to the function u(r ) found in Equation (2.9) of the paper by Tooper [2], but it is used here in a completely different context. Note that the asymptotic boundary condition on λ(r ), that it must behave as the exterior Schwarzschild solution for sufficiently large r , translates here as β(r ) → 1 under that same condition. Substituting these expressions in the component field equation shown in Eq. (9) a very simple relation giving the derivative of β(r ) in terms of ρ(r ) results, Therefore, the energy density ρ(r ) is given in terms of the derivative of β(r ), and wherever ρ(r ) = 0, characterizing a region where there is a vacuum, we have that β(r ) is a constant. Since by Eqs. (6) and (7) F(r ) and P(r ) are both given in terms of ρ(r ), and since λ(r ) is given in terms of β(r ), it follows at this point that, given a function β(r ), the functions ρ(r ), P(r ), F(r ) and λ(r ) are all determined. The only function that has yet to be determined in terms of β(r ) is ν(r ). We can obtain ν(r ) in terms of F(r ), and therefore of β(r ), using the consistency condition in Eq. (11), which with the use of Eq. (8) can be written as Note that Eq. (11) cannot be used within the vacuum regions, but only within the matter region. In the vacuum regions one must use Eq. (10) instead, in order to obtain ν(r ). In all cases in which there is an outer matter-vacuum interface at r 2 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 ρ(r 2 ) = 0, we get which determines ν(r ) in terms of F(r ), up to the integration constant ν(r 2 ), that is to be obtained from the asymptotic boundary conditions at radial infinity, which in the case of ν(r ) is simply ν(∞) = 0. Therefore, given β(r ), this determines ν(r ) in terms of it, through F(r ). Note, for future use, that in all cases in which there is an inner matter-vacuum interface at r 1 as well, the fact that we also have that F(r 1 ) = 0, because for n > 0 we have the boundary condition ρ(r 1 ) = 0, implies that we always have that ν(r 1 ) = ν(r 2 ). We see therefore that the determination of the function β(r ) leads with no further difficulty to the determination of all the functions that describe both the matter and the geometry of the system, 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 β (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 Eq. (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 as well as to define the dimensionless function of ξ , to assume the role of β(r ), Note that the asymptotic condition that β(r ) → 1 for sufficiently large r translates here as the condition that γ (ξ) → ξ M for sufficiently large ξ . It is important to observe that under the change of variables from r to ξ the homogeneous derivatives transform in a very simple way, for example in the case of β(r ), We will also introduce at this point a notation for the derivative of γ (ξ), which will be useful in order to deal with the second-order differential equation for γ (ξ) that we will arrive at, where in this case the prime denotes differentiation with respect to ξ . 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 γ (ξ) and π(ξ ), that will be treated as independent variables. In terms of these new variables we have for the system of differential equations shown in Eqs. (9)- (11), where from now on, given that we have the expressions in Eqs. (6) and (14), the auxiliary function F(r ) will be written as F(ξ, π ), where the primes now indicate differentiation with respect to ξ . We therefore have for all the relevant quantities written in terms of ξ , γ (ξ) and π(ξ ), where C = K / κr 2 0 1/n is a dimensionless constant. We see therefore that the determination of γ (ξ) and thus of π(ξ ) leads to the complete solution of the problem. We have therefore reduced the solution of the problem to the determination of the single function γ (ξ). In order to determine γ (ξ), we obtain an ordinary differential equation for it by eliminating ν (ξ ) from Eqs. (28) and (29). Equation (28) can be written as and Eq. (29) can be written as so that equating the two right-hand sides we get Since π(ξ ) is the derivative of γ (ξ), this can be interpreted either as a second-order ordinary differential equation for γ (ξ), or as one of a pair of first-order coupled ordinary differential equations determining γ (ξ) and π(ξ ), the other equation being simply 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 Sect. 4. However, before we plunge into the numerical approach, let us first show that γ (ξ) is in fact a very simple function, as is β(r ), and that it has some properties which are, perhaps, somewhat unexpected.

Main properties of the solutions
Let us start by recalling that, if it turns out that there are indeed inner and outer vacuum regions, then we already know the form of the solutions in those regions. From the derivations in [1], we have that in the inner vacuum region, where ξ < ξ 1 , the solution is given, in terms of our current set of dimensionless variables, by where ξ μ = r μ /r 0 , and where r μ and A are integration constants, and that in the outer vacuum region, where ξ > ξ 2 , the solution is the exterior Schwarzschild solution, which can be written, in terms of our set of dimensionless variables, as where ξ M = r M /r 0 and r M is the Schwarzschild radius. The constant A can be determined with the use of the consistency condition and of the interface boundary conditions for P(ξ ), which as we saw before imply that we always have ν(ξ 1 ) = ν(ξ 2 ) within the matter region. This fact allows us to determine the value of A, using also the fact that ν(ξ ) is a continuous function, thus leading to ν(ξ 1 ) = ν i (ξ 1 ) and ν(ξ 2 ) = ν s (ξ 2 ), which imply that we have This gives us the solution for A in terms of the other quantities, We can therefore write out the complete solution in both vacuum regions, The quantities ξ 1 and ξ 2 used in these sectors of the solutions are obtained during the resolution of the differential equation within the matter region, as are the other two quantities, ξ μ and ξ M , which once γ (ξ) is determined are given by as consequences of the interface boundary conditions related to the continuity of λ(ξ ) across the interfaces. The consequences of the interface boundary conditions related to the continuity of ν(ξ ) are reflected on the condition that the matter energy density and the pressure must be zero at the matter-vacuum interfaces, thus leading to the result shown in Eq. (45). It is important to emphasize that when we use the values of ξ μ and ξ M obtained in the matter region in order to define the geometries of the inner and outer vacua, we are in effect satisfying the interface boundary conditions. We may now derive some crucially important properties of the solutions of the field equations analytically, by obtaining these properties directly from the differential equations.
Property (1). The first important fact about the solution γ (ξ), away from the origin, comes from Eqs. (14) and (24). From the first of these, since all quantities on the right hand side are non-negative, and the only one that can be zero away from the origin is ρ(r ), it follow that β (r ) is a non-negative function, that is zero only within a vacuum region. It therefore follows that β(r ) is a monotonically increasing function, which is a constant if and only if we are within a vacuum region. In addition to this, since according to the asymptotic boundary condition we have that β(∞) = 1, it also follows that β(r ) is limited from above by 1. It then follows from Eq. (24) that similar conclusions can be drawn for γ (ξ), which is therefore a monotonically increasing function which is limited from above, the limit in this case being the parameter ξ M , and which is constant within vacuum regions. There are two more important properties of the solutions γ (ξ) and π(ξ ) of Eqs. (37) and (38), away from the origin, that we can establish analytically. Let us therefore consider Eq. (37), in which F(ξ, π ) is given as in Eq. (32). Let us also recall that we are working under the assumption that we have n ≥ 1.

Property (2).
For the second important property, if we assume that there is a radial position ξ = ξ e , where the subscript e stands for extremum, at which we have that π (ξ e ) = 0 and π(ξ e ) = 0, then there is a definite and unique solution to these conditions that presents itself, namely the expression within braces in Eq. (37) must vanish, leading to Since we already know that γ (ξ) must be a monotonically increasing function that is limited from above by ξ M , the condition π (ξ ) = 0 with π(ξ ) = 0 identifies the inflection point ξ e of γ (ξ), which is also the point of maximum of π(ξ ). The equation above then identifies in a definite way the value of γ (ξ e ) as a certain function of π e = π(ξ e ) at that inflection point, a relation that is given explicitly by This relation is the second important property of the solutions, which identifies the radial position ξ e of the inflection point. This constitutes a radial position ξ e and a pair of values γ (ξ e ) and π(ξ e ) of γ (ξ) and π(ξ ) that can be used as the starting point for a numerical integration of the equation, to either one of the two sides of ξ e . When there is a well-defined matter region, it allows us to perform the numerical integration starting from a regular point in the interior of the matter interval, rather than at one of the two ends of that interval which, as we will see shortly, are soft singular points.

Property (3).
The third important property of the solutions is that, if there is a radial position ξ = ξ b where π(ξ ) = 0, then Eq. (37) implies that π (ξ ) = 0 at that radial position. Since the equation for γ (ξ) is a second-order one, and in this situation we have that both π (ξ ) = 0 and γ (ξ ) = 0 when we approach that radial position integrating from one side, then on the other side of that radial position we will have that both π(ξ ) and γ (ξ) are constant, and in particular that π(ξ ) = 0 is the null constant. Since we have that the energy density ρ(ξ ) is proportional to π(ξ ), this implies that on that other side we have a vacuum, so that the radial position ξ b represents the boundary of the matter. Note that this also implies that, if there is such a radial position ξ b , then we must conclude that the dynamics of the system leads to the formation of a sharp boundary for the matter, even in the case of a gas, so long as it satisfies the polytropic equation of state exactly.
Let us now consider the consequences of Eq. (37) regarding the manner in which such a radial position ξ b is approached by π(ξ ). We must consider separately the cases in which the radial position is approached from the left and from the right. We start by the case in which the radial position ξ b is approached by the integration process from the right, from values ξ > ξ b . We will assume that the function π(ξ ) goes to zero at ξ b as a power m ≥ 1, and verify whether we can find a solution of Eq. (37) in a rightneighborhood of ξ b . Since we know that π(ξ ) must be real and positive, while m is not necessarily an integer, this means that we should assume that in a right-neighborhood where we must have that B ⊕ > 0 and that m ≥ 1, given that π(ξ ) is real and nonnegative and that π (ξ ) must exist. The question is then whether or not we can find values of m and B ⊕ such that the equation is satisfied in that neighborhood. If we denote the difference shown by δξ = ξ − ξ b , we may then write for all the quantities involved in Eq. (37), in a right-neighborhood of ξ b , where D ⊕ = C B m/n ⊕ . Substituting all the quantities in Eq. (37), recalling that we are interested only in the δξ → 0 limit, assuming that γ (ξ b ) = 0, which allows us to use for γ (ξ) just its dominant term γ b = γ (ξ b ), and dropping negligible terms where possible, we get It then follows from the δξ → 0 limit that the only possible solution with m = 0 is m = n, thus confirming that m ≥ 1. Recalling the value of D ⊕ , we then have for B ⊕ Note that, since B ⊕ must be strictly positive, it follows from this that we must have γ b < 0, thus confirming a posteriori that γ b = 0, so long as we are not dealing with just the identically null solution for π(ξ ). Since in the resulting internal vacuum we have that γ b = −ξ μ = −r μ /r 0 , we once more conclude that we must have r μ > 0, just as we did in [1], but in a very different way.
It is important to note that, just as was the case in that previous paper, this result depends on the condition that r 1 > 0, because otherwise there would be no inner vacuum region, and hence no interface boundary conditions at r 1 , so that all the derivations done up to this point would cease to be valid. It is interesting to observe that what we have proved here about r μ can be stated as which then implies that we must also have, just as in the previous paper, that This means that making r μ = 0 takes us from the shell solutions presented here to the Tooper solutions for filled spheres [2]. We may also identify that in this case we have ξ b = ξ 1 , the inner boundary of the matter, and hence we may write that We now consider the case in which the radial position ξ b is approached by the integration process from the left, from values ξ < ξ b . Again we will assume that the function π(ξ ) goes to zero at ξ b as a power m ≥ 1, and verify whether we can find a solution of Eq. (37) in a left-neighborhood of ξ b . Since we know that π(ξ ) must be real and positive, while m is not necessarily an integer, this means that this time we should assume that in a left-neighborhood of ξ b we have where once again we must have that B > 0 and that m ≥ 1. If we now denote the difference shown by δξ = ξ b − ξ , we may then write for all the quantities involved in Eq. (37), in a left-neighborhood of ξ b , where D = C B m/n . As before, substituting all the quantities in Eq. (37), recalling once again that we are interested only in the δξ → 0 limit, assuming that γ (ξ b ) = 0, which allows us to use for γ (ξ) just its dominant term γ b = γ (ξ b ), and dropping negligible terms where possible, we get m = n n + 1 It follows once more from the δξ → 0 limit that the only possible solution with m = 0 is m = n, once again confirming that m ≥ 1. Recalling the value of D , we then have for B Note that, since B must be strictly positive, in this case it follows that we have γ b > 0, once again confirming a posteriori that γ b = 0. Besides, the quantity γ b can be written in terms of the parameters of the resulting external vacuum as γ b = ξ M = r M /r 0 , which also implies that it cannot be zero. We may also identify that in this case we have ξ b = ξ 2 , the outer boundary of the matter, and hence we may write that where we necessarily have that ξ 2 > ξ M since the integration cannot produce the coordinate singularity of a horizon. This gives us the exact asymptotic behavior of π(ξ ), and hence also of γ (ξ), as we approach the radial positions ξ 1 and ξ 2 that delimit the region containing the matter, from within that region. Whether or not the points ξ 1 and ξ 2 exist in each particular case, for various values of the parameters that define the physical system, has to be determined numerically. If they do, then there is more that can be established about γ (ξ). In this case, not only is this a monotonically increasing function that is limited from above, but we now know that it is also limited from below, since it is a constant within the inner vacuum region. In fact, it must go from a constant negative value within the inner vacuum region, to a constant positive value in the outer vacuum region. In addition to this, so long as n ≥ 0 it is a continuous and differentiable function on the whole positive real semi-axis. In particular, it follows from this that γ (ξ) must have a single zero in the open interval (ξ 1 , ξ 2 ), which therefore is always within the matter region. This completes the discussion of the third property.

Property (4).
We can easily determine the circumstances in which the weak and dominant energy conditions [7] are satisfied for the solutions we present here. Since in our case here both g μν and T ν μ are diagonal, with diag[g μν ] = e 2ν(r ) , − e 2λ(r ) , −r 2 , −r 2 sin 2 (θ ) , it follows that, as was explained in the previous paper [1], the weak energy condition is always satisfied, for all values of the parameters of the model. In addition to this, as was also explained in that paper, the dominant energy condition is satisfied so long as one has the condition that ρ(r ) ≥ P(r ), for all r within the matter region. This is equivalent to the condition that In our specific case here we have that the polytropic equation of state between ρ(r ) and P(r ) given in Eq. (5) holds, and that it can be written as which implies that P(r ) ρ(r ) = C κr 2 0 ρ(r ) Recalling from Eq. (30) that we have that κr 2 0 ρ(r ) = π(ξ )/ξ 2 , where ξ = r /r 0 , it follows that the dominant energy condition reduces to This must be true for all values of ξ within the matter region. It suffices therefore to impose this condition only at the point where ρ(r ) is maximum, that is where π(ξ )/ξ 2 is maximum. However, as one can see in Fig. 4, for the larger values of ρ(r ) the points of maximum of ρ(r ) and of π(ξ ) tend to coincide, as do the values of these two quantities at their points of maximum, if we use the value ξ e = 1 for the point of maximum of π(ξ ), which we can do because it corresponds simply to a choice of the arbitrary parameter r 0 as the position r e of the inflection point. Therefore, if we consider the point ξ e = 1, the position of the inflection point of γ (ξ) at which we start the numerical runs, where π(ξ ) = γ (ξ ) has its global maximum π e , then we see that it is a fair estimate for high-density solutions to consider the dominant energy condition only at this point, and thus we arrive at a simple constraint among the free parameters of the model, which therefore establishes a physically valid region within the (π e , C, n) parameter space of the model, for the smaller values of C, and correspondingly larger values of π e . In the general case it may be necessary to find the maximum value of π(ξ )/ξ 2 by numerical means, in order to determine with precision the maximum possible value of π e for each value of C. In any case, we can state with confidence that there is a region in the parameter space of the solutions where the dominant energy condition is satisfied.

Property (5).
Since the solutions presented in this paper have the same form, within the inner vacuum region, as the solutions presented in [1], they therefore have exactly the same singularity at the origin ξ = 0. It follows therefore that the same calculations and comments about the Riemann curvature tensor, the left-dual curvature tensor and the double-dual curvature tensor, as well as about the various curvature scalars that can be derived from them, that were presented in [1], are equally valid here. It is worth the trouble exploring now some of the consequences of the facts just established about γ (ξ) and π(ξ ). Note that the two boundary points ξ 1 and ξ 2 are soft singular points of both functions γ (ξ) and π(ξ ), since these functions cannot be analytic at these points. This follows from the fact that they are constant on one side of the points and behave as a strictly positive power on the other side. Due to this there is no power series centered at ξ b that can converge to these functions on a neighborhood of these points, and therefore they are not analytic. However, both functions are still well-defined at the boundary points, thus characterizing the singularities as soft ones.
Added to the fact that the position of the boundary points is not known beforehand, this makes it impossible to integrate the equation numerically starting at these points. Since the equation for γ (ξ) is a second-order one, and since in this situation we have that both π (ξ ) = 0 and γ (ξ ) = 0 hold at the positions of the boundaries, if we try to integrate from these points into the matter region we will simply get a constant for γ (ξ) and the constant value π(ξ ) = 0, which corresponds to a vacuum rather than to the matter that is there. Therefore, the solution within the matter region is not determined by its values at the boundaries, and hence this problem cannot be characterized as a typical boundary value problem.
Since n ≥ 1, the two functions π(ξ ) and γ (ξ) will be continuous at these boundary points, but there will be some higher-order derivative that does not exist there. For example, if n = 1 the derivative π(ξ ) = γ (ξ ) will be continuous but not differentiable at these points, so that the second derivative π (ξ ) = γ (ξ ) does not exist there. A relevant example is n = 3/2, which is typical of the convective layer of a star, in which case π (ξ ) will not be differentiable at these points, although it is still continuous there, so that the third derivative π (ξ ) of γ (ξ) will not exist, since it diverges at these points. As another relevant example, if n = 3, which is typical of the radiative layer of a star, then the fourth derivative π 3 (ξ ) of γ (ξ) will be discontinuous at these points, so that its fifth derivative π 4 (ξ ) will not exist there.
This same singular character of the interface boundary points has the global consequence that, given definite boundary conditions at radial infinity, the complete solution of the problem is not determined, and is therefore not unique. On the one hand, if one integrates inwards from a point ξ > ξ 2 , then one produces just a continuation of the exterior Schwarzschild solution of the outer vacuum region, until one reaches its horizon, and never any solution associated to the matter. On the other hand, if one integrates outwards from a point ξ < ξ 2 , then there will be many sets of values of the parameters that describe different states and characters of the matter, but which are such that the solution arrives at ξ 2 with the same π(ξ 2 ) = 0 and γ (ξ 2 ) = ξ M . There are therefore many interior solutions that correspond to the same exterior solution. This is simply an instance of the uniqueness theorem of the exterior Schwarzschild solution, shown by Jebsen and Birkhoff [8,9].
The numerical analysis indicates that the boundary points ξ 1 and ξ 2 exist for all values of the parameters of the system, within a wide range of variation of these parameters, within physically reasonable bounds. The boundary point ξ 2 seems to always exist, for all physically allowed values of the parameters. The boundary point ξ 1 seems to exist for most sets of values of the parameters. When it does not exist, the integration for the matter density diverges to infinity towards the origin, so that there is no inner vacuum region, and a truly hard singularity develops at the origin, with an infinite concentration of matter there. In such cases there is no acceptable solution at all. This seems to indicate that, if a static solution exists at all, then it has the property that these two boundary points are present. Next we will describe in detail the numerical approach that leads to this conclusion.

Examples of numerical solutions
Let us now describe in detail our strategy for the numerical solution of Eq. (37). We will consider this to be a system of two coupled first-order differential equations for the functions γ (ξ) and π(ξ ), leading to a pair of numerical development equations, which can be easily obtained, and that are given by where F is the value of the auxiliary function F(ξ, π ), where δξ is a small increment of the variable ξ and where δγ and δπ are the corresponding increments of the functions γ (ξ) and π(ξ ). This system can be integrated by the use of the Runge-Kutta fourthorder algorithm in a straightforward way, so long as one has a sensible starting point for the integration process.
We must now discuss where to start the integration process. There are two alternatives that we have used. The first alternative is to start at the inflection point ξ e of the function γ (ξ). As we saw in Sect. 3, Eq. (37) itself provides us with a fixed relation between γ (ξ e ) and π(ξ e ) at this point. This means that, given a value of π(ξ e ) at this point, which will have the role of a free parameter of the system, we have the corresponding value of γ (ξ e ), which then allows us to start the integration process. The second alternative is to start at the root ξ z of γ (ξ), which we know necessarily exists and is necessarily located within the matter region. Once again the value of π(ξ z ) at this point will have the role of a free parameter of the system, and in this case we just start with γ (ξ z ) = 0. The first alternative seems to work better for the larger initial values of π(ξ e ), which corresponds to the more dense objects, and the second alternative seems to work better for the smaller initial values of π(ξ z ), which corresponds to the less dense objects.
In either case, we must start with a choice for the initial value of ξ . Taking the position ξ e = r e /r 0 of the inflection point as the example of a starting point for this discussion, we note that a choice of value for ξ e is tantamount to a choice of a relation between r e and r 0 , without implying the choice of a definite value for either one. Therefore the choice of the initial value of ξ can be quite arbitrary, without leading to any loss of generality. For example, if we start with ξ e = 1, this only means that we choose the arbitrary parameter r 0 = r e to be the position of the inflection point, without actually choosing a value for r e . Therefore, a solution found in this way is not just a single solution, but a one-parameter family of solutions, parametrized by the values of r 0 . A similar argument can be made in the case in which we start at the position ξ z = r z /r 0 of the zero of γ (ξ).
Note that the physical constants K , κ and r 0 do not appear individually in either the numerical propagation equations or the initial values. They appear only as the combination seen in the definition of the dimensionless constant C, which in turn appears only in the expression for F(ξ, π ), We will therefore adopt C as one of our free parameters for the numerical work, the others being n and the initial value of π(ξ ), which will then result in a certain value of ξ M . Note that the parameters ξ M that contains the mass M also does not appear in either the equation or the initial values. It appears only as the asymptotic boundary condition for γ (ξ), which is obtained only at the end of the integration process in the outward direction.
Departing from the chosen initial point, we then integrate out in both directions, of increasing and decreasing ξ , until π(ξ ) approaches zero in each case, with a high degree of numerical precision. Subsequently we concatenate the results into a single function from ξ 1 to ξ 2 . The values obtained for γ (ξ) at the two boundary points are then used to generate analytically the correct exact solutions in the inner and outer vacuum regions, which are plotted alongside the corresponding numerical solution, from some point to the left of ξ 1 to some point to the right of ξ 2 . In doing this we are in effect satisfying the interface boundary conditions at the matter-vacuum interfaces.
In order to do this, at the outer boundary we just use the fact that as given in Eq. (51), in order to generate the correct exterior Schwarzschild solution as given in Eqs. (41) and (42), thus satisfying the interface boundary conditions. In order to obtain the dimensionfull physical parameters, we must recall that We can then simply choose any positive value that we want for r M , and thus obtain the corresponding value of the parameter r 0 = r M /ξ M . This value of r 0 can then be used to obtain the values of all the other dimensionfull physical parameters of the system. In the case of the inner boundary we use the fact that as given in Eq. (50), as well as the fact that we have for the integration constant A in order to generate the correct interior vacuum solution as given in Eqs. (39) and (40), thus satisfying the interface boundary conditions, with Note that the integration constant for the function ν(ξ ) within the matter region must be chosen so that at the outer interface we have that ν(ξ 2 ) = ν s (ξ 2 ) has the correct value, given by Eq. (42). This is easily accomplished by just correcting the values of ν(ξ ) afterwards, by simply adding a constant to them. From Eq. (32) we can see that, since at the outer interface we have that F(ξ 2 , π 2 ) = 0, where π 2 = π(ξ 2 ), it follows that if we just ignore the integration constant we get ν(ξ 2 ) = 0. Therefore, all that we have to do is to add ν s (ξ 2 ) to ν(ξ ) afterwards, for which we may take advantage of the fact that ν s (ξ 2 ) = −λ s (ξ 2 ) is already known. We have used three sets of values of the parameters in order to generate the data seen on the graphs, the first two starting from the inflection point of the function γ (ξ), and the last one starting from the zero of that function. These sets of parameters are as follows. The first set uses some arbitrary mid-range values of the parameters, that have the property of just displaying in a simple and clear way all the main characteristics of the solutions. The second set of parameters corresponds qualitatively to the configuration of a high-density object such as a white dwarf or neutron star. In these two cases the integration was started from the inflection point ξ e of γ (ξ). The third set of parameters corresponds qualitatively to the configuration of a low-density object, such as a normal main sequence star. In this case the integration was started from the zero ξ z of γ (ξ).

Analysis of the numerical results
The graph seen in Fig. 1 shows most of the functions involved in the solution, for a typical mid-range set of the dimensionless parameters, namely n = 3/2, C = 0.05, and π e = 1.0. The crucial function, from which everything else stems, is the function γ (ξ), which in all cases has the same qualitative behavior shown in that graph, for any values of the parameters for which a solution exists. It is a very simple function, that slopes up monotonically from a constant negative value to a constant positive value. These two constant values determine ξ μ and ξ M respectively. Its derivative π(ξ ) is also a simple function, with a single well-defined point of maximum. It is closely related to the matter energy density ρ(ξ ), which is also shown in the graph, but in its dimensionless version given byρ(ξ ) = κr 2 0 ρ(ξ ), which is the same as π(ξ )/ξ 2 , as one can see in Eq. (30). The second derivative π (ξ ) gives the numerical propagation function. For the value n = 3/2, which is used in this case, it clearly marks the positions of the interface points, where its graph hits the ξ axis at right angles.
The general behavior of the solution is that almost always there are two interface points and thus both an inner vacuum region and an outer vacuum region. However, the parameters can be judiciously adjusted so as to decrease the width of the inner vacuum region to zero, in which case the shell becomes a filled sphere. In this case one gets ξ 1 = 0, so that the inner interface tends to the origin, and then one gets ξ μ = 0 as well, so that we have that γ (ξ) is zero at the origin, that is, the Tooper boundary condition γ (0) = 0 holds, as is the case in most treatments, such as in [7,[10][11][12]. Therefore, this takes us back to the solutions found by Tooper [2].
For each given value of C there is a minimum positive and non-zero value of π e that is allowed, which is the one that gives the Tooper solution. Below that minimum  value of π e a hard singularity of π(ξ ) is generated, corresponding to an infinite matter energy density, and then there is no acceptable solution to the problem. On the other hand, the outer interface point seems to always exist, for all allowed values of the parameters of the system. The only limitation on the values of the parameters due to existence conditions for the solutions are those related to the inner interface.
The graph in Fig. 2 shows the functions ν(ξ ) and λ(ξ ) that describe the geometry of the solution for this same set of input parameters. Outside the outer interface these are the just the functions of the exterior Schwarzschild solution, for which we have that ν(ξ ) = −λ(ξ ), with λ(ξ ) > 0 and ν(ξ ) < 0. Somewhere within the matter region there is a crossing of the graphs of ν(ξ ) and λ(ξ ), so that for sufficiently small ξ these signs are reversed, and we then have that λ(ξ ) < 0 and ν(ξ ) > 0. Inside the inner interface these are the functions of the exact solution for the inner vacuum. Since these two functions have singularities at the origin, the graphs are limited to a region within which the graphs stay below certain maximum absolute values of ν(ξ ) and λ(ξ ). It is worth emphasizing that these are the only two functions involved that The graph in Fig. 3 shows the quantity log(π ) as a function of log(ξ − ξ 1 ), for this first set of values of the parameters. This graph depicts the behavior of the function π(ξ ) in the approach of the inner interface at ξ 1 from within the matter region, which proceeds from the right to the left in the graph. In this limit π(ξ ) behaves indeed as (ξ − ξ 1 ) 3/2 , as one can see from the fact that the slope of the straight line shown in the graphs is exactly 3/2, within the numerical precision level in use. This straight line was obtained by linear regression from the numerical data. This result confirm the analysis made in Sect. 3, and it also shows that π(ξ ) in fact hits zero at ξ 1 , since otherwise this log-log graph could not turn out to be a straight line.
In the graph shown in Fig. 4 we display a high-density case, in which we use the value n = 3, a smaller value of the parameter C, with C = 0.01, and a much larger value of the parameter π e , with π e = 100.0. In this case the inner vacuum is very wide in the radial direction, and all the matter is concentrated within a relatively thin shell located closer to the outer interface than to the inner one.
As one can see in the graph in Fig. 5, in this case the geometry becomes much more extreme. The constant value of γ (ξ) within the inner vacuum regions becomes large and negative. For large but finite values of π e this significantly decreases the physical volume of the inner vacuum region, as compared to its apparent volume. Note that within the inner vacuum region, since the radial lengths shrink, while the angular ones remain invariant, the geometry of a two-dimensional spatial section through the origin, of this four-dimensional geometry, is not embeddable in a flat three-dimensional space, in the way that the exterior Schwarzschild geometry is.
The data for the third set of parameters, shown in the graphs is Figs. 6, 7, and 8, a low-density example with n = 3/2, C = 1/3, and π z = 10 −3 , shows how our solution approaches the behavior that is expected by our classical Newtonian intuition, along most of the matter distribution. This example is well along the limit leading to the Tooper solutions, as indicated by the very small value of ξ μ given in the caption of   Fig. 6. In this case the inner vacuum region and the singularity that it contains become reduced to a very small region near the origin. The smaller the value of π z , the smaller this region becomes, and in the π z → 0 limit it is reduced to a single point. Under appropriate conditions, with small but non-zero π z , this very small singular region may then be washed away by the thermal fluctuations of the system, since these fluctuations certainly violate the spherical symmetry at a sufficiently small length scale, and will cause the small region containing the singular point to fluctuate rapidly around the origin, thus smearing that singular point.
In this case the four-dimensional geometry becomes almost flat almost everywhere, since both ν(ξ ) and λ(ξ ) become very small, as can be seen in the graph shown in Fig. 7. What curvature there is becomes quite smooth, and embeddable in three dimensions, all the way down to the point where λ(ξ ) becomes zero. This embedding only fails to be possible within a very small region near the origin. In Fig. 8 one can see the matter energy density, which has a maximum quite close to the origin, and is monotonically decreasing outward from there, just as is hypothesized in the case of the Tooper solutions. In the π z → 0 limit the radial position of the maximum of the matter energy density tends to the origin.
Note that while in the low-density case the inner vacuum region that contains the singularity becomes ever smaller as π z decreases, in the case of the data for the second set of parameters, with a large value of π e , the inner vacuum region and the effects of the central singularity spread throughout most of the interior of the object, and thus cannot be ignored. Therefore, while for low-density objects our solution does not differ significantly from the Tooper solutions, for dense objects it is rather dramatically different from it.

Conclusions
In this paper we have given the complete static solution of the Einstein field equations for the case of a spherically symmetric distribution of gaseous matter satisfying the equation of state of a polytrope. We have arrived at some of the same important and rather unexpected conclusions that we had already come across in a previous paper [1] on the solution for shells of liquid matter.
One new fact, which is different from that previous case, is that in this case the inner and outer boundaries of the matter, that define a spherical shell, are not imposed by hand, but arise as inevitable consequences of the dynamics of the system, as determined by the Einstein field equations. This makes it impossible to ignore these solutions, for there is no arbitrary choice of a geometrical character involved. All the free parameters of the system are related only to the physical characteristics of the matter, and all geometrical characteristics that arise are inevitable consequences of the equations.
The new solutions converge back to the known Tooper [2] family of solutions in certain limits of the free parameters of the physical system, in which the shell becomes a filled sphere. While for each value of the index n that appears in the polytropic equation of state the Tooper family of solutions is described by a two-dimensional parameter space, the family that we present here is described by a three-dimensional parameter space.
One of the conclusions that is just like those of the previous case involving liquid matter is that all solutions of this type have a singularity at the origin, within the inner vacuum region, that does not, however, lead to any kind of pathological behavior involving the matter. Note that the resulting matter energy density is not imposed by hand, but is instead a consequence of the dynamical equations of the system. Therefore, the fact that the matter energy density is zero at the origin, where the singularity lies, is a consequence of those equations. The other is that, contrary to what is usually thought, a non-trivial gravitational field does exist within a spherically symmetric central cavity, namely the inner vacuum region. This field can be characterized as being repulsive with respect to the origin, which explains why we do not see an infinite concentration of matter at that point.
Unlike the problem examined in the previous paper [1], which can be seen as having a somewhat academic character, the problem we examine here can have direct applications to astrophysical objects. Not only we can use the isentropic case n = 3/2 to represent the external convective layer of a star, and the case n = 3 to represent its internal radiative layer, but we might also consider using both, tied up to one another by means of appropriate interface boundary conditions at an intermediary point, in order to represent a star in a more complete fashion, including the two main layers that are known to exist.
Because these new solutions hold over the entire four-dimensional manifold, in certain limits they have interesting consequences regarding the concept of a black hole, and specially regarding the geometry of its interior region. However, a detailed discussion of that topic would be excessively long to be included here, since it would involve mapping out the whole parameter space of all these possible static solutions. Therefore, that discussion will be presented in a separate paper.