The Parameter Planes of the Spherically Symmetric and Static Relativistic Solutions for Polytropes

The Parameter Planes of the Spherically Symmetric and Static Relativistic Solutions for Polytropes Jorge L. deLyra∗ Universidade de S˜ao Paulo Instituto de F´ısica Rua do Mat˜ao, 1371, 05508-090 S˜ao Paulo, SP, Brazil December 20, 2023 Abstract We explore the parameter space of the family of static and spherically symmetric solutions of the Einstein field equations for polytropes, that were presented in a previous paper. This is a four-parameter family of exact solutions, of which one parameter can be factored out, so that there are only three essential free parameters. The solutions are exact in the sense that no approximations are involved, other than those implied by the numerical precision limitations. The primary objectives of this exploration are to establish directly the existence of large collections of specific solutions, and to determine some of their most important properties. For each value of one of the three essential free parameters of the family of solutions, the polytropic index n, taken here, for the sake of simplicity, to be either an integer or a half-integer, we define and explore the parameter planes spanned by the other two essential free parameters. In this way, besides establishing their existence, we are also able to classify the solutions according to their overall matter energy density, as well as in terms of their proximity to solutions that display an event horizon. For four values of nwe successfully establish the allowed regions of the parameter planes, where the solutions not only exist but also correspond to physically acceptable matter. We find that there are solutions within these regions with overall matter energy densities varying all the way from very low to very high, including some that are as close as one may wish to solutions that display an event horizon, and that therefore represent black holes, or extremely dense objects very similar to them. ∗Email: [email protected] ORCID: 0000-0002-8551-9711 1 Keywords: General Relativity, Einstein Equations, Solutions for Polytropes, Event Horizons, Energy Density. DOI: 10.5281/zenodo.7871604 Disclaimer: A version of this article has been accepted for publication, after peer review, but this is not that Version of Record, neither does it reflect any corrections. It does, however, include a few post-acceptance improvements. The Version of Record is available online at: http://dx.doi.org/10.1007/s10714-023-03187-4. 2 1 Introduction In a previous paper [1] we established the static solution of the Einstein field equations for the case of spherically symmetric shells of gaseous fluid located between radial positions r1 and r2of the Schwarzschild system of coordinates. These solutions are exact in the sense that they involve no approximations, but they are not written entirely in closed analytical form. While all relevant functions describing both the geometry and the matter are written analytically in terms of a single function, which we denote by β(r), and although the main properties of this function can also be established analytically, it is still necessary to determine this function in detail by numerical means, including the existence of the internal and external radii r1and r2, and therefore the existence of the solution itself. In this paper we will present a fairly extensive numerical exploration of the parameter space of these solutions, in order to establish the general picture in regards to their existence and to their main properties. Connections with some special cases which are already well known will be pointed out. Some diagnostic observables will be defined and calculated, which were designed to put in evidence some of the most important properties of the solutions. In this paper we will be limited to the case in which the polytropic index n, which is one of the free parameters, is strictly larger that 1, that is, satisfies n > 1. This paper is organized as follows: in Section 2 we review the new class of static and spherically symmetric exact solutions for gaseous shells; in Section 3 we discuss the parameter planes for a few values of the polytropic index n, and examine in particular the allowed regions in each case; in Section 4 we discuss the definition and the behavior of the diagnostic observables over the allowed regions, and thus establish the general character of the solutions; in Section 5 we discuss in general terms how to apply the solutions to problems in astrophysics; in Section 6 we provide some further analysis of the solutions; and in Section 7 we state our conclusions. 2 Review of the Polytropic Solutions In this section we will review the solutions for gaseous fluids presented in [1], in order to establish the notation, as well as the most basic character of the solutions. Both the notation and the ideas presented in [1] rely on some of the notation and ideas presented in the previous paper [2], in which we presented the exact solution for the case of liquid fluids. In this work we will use the time-like signature (+,−,−,−), following [3]. In terms of the coefficients of the metric, for a static invariant interval given in terms of the Schwarzschild coordinates (t, r, θ, φ) by ds2= e2ν(r)c2dt2−e2λ(r)dr2−r2dθ2+ sin2(θ)dφ2,(1) where exp[ν(r)] and exp[λ(r)] are two positive functions of only r, the Einstein field equations reduce to the set of three first-order differential equations n1−2rλ′(r)oe−2λ(r)= 1 −κr2ρ(r),(2) n1 + 2 rν′(r)oe−2λ(r)= 1 + κr2P(r),(3) [ρ(r) + P(r)] ν′(r) = −P′(r),(4) where ρ(r) is the energy density of the matter, P(r) is its isotropic pressure, and the constant κis given by κ= 8πG/c4, where Gis the universal gravitational constant and cis the speed of light. In these equations the primes indicate differentiation with respect to r. 3 It is convenient for the analysis of the solutions to change variables in the field equations from the function λ(r) to a function β(r), which is defined to be such that e2λ(r)=r r−rMβ(r),(5) where rM= 2GM/c2is the Schwarzschild radius associated to the total asymptotic gravitational mass M, which then implies that we have for the corresponding derivatives 2rλ′(r) = −rM β(r)−rβ′(r) r−rMβ(r).(6) Substituting the expressions in Equations (5) and (6) in the component field equation shown in Equation (2) a very simple relation giving the derivative of β(r) in terms of ρ(r) results, β′(r) = κr2ρ(r) rM .(7) Therefore, in any interval where ρ(r) = 0 we have that β(r) is a constant. Since we must have that ρ(r)≥0, 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 must have that β(r)→1 when r→ ∞, it also follows that β(r) is limited from above by 1. In the previous paper [1] we established the static solution of the Einstein field equations for the case of a spherically symmetric shell of gaseous fluid located between the radial positions r1and r2of the Schwarzschild system of coordinates. These positions are not arbitrary, but rather are obtained as part of the solution of the problem. For this problem we assume the hypothesis that the gas satisfies the polytropic equation of state P(r) = K[ρ(r)]1+1/n ,(8) where K, the polytropic constant, is a positive real constant, and n > 1, the polytropic index, is a real number that, merely for simplicity, may be taken to be an integer or halfinteger. For convenience we define the auxiliary quantity F(r) = K[ρ(r)]1/n ,(9) in terms of which the equation of state becomes simply P(r) = F(r)ρ(r).(10) Given the field Equations (2) through (4) and the equation of state shown in Equation (8), the solution is given by λ(r) =                    −1 2lnr+rµ rfor 0 < r ≤r1, −1 2lnr−rMβ(r) rfor r1≤r≤r2, −1 2lnr−rM rfor r2≤r < ∞, (11) 4 ν(r) =                  1 2ln1−rM/r2 1 + rµ/r1+1 2lnr+rµ rfor 0 < r ≤r1, ν(r2)−(n+ 1) ln[1 + F(r)] for r1≤r≤r2, 1 2lnr−rM rfor r2≤r < ∞. (12) The solution introduces into the system the new physical parameter rµwith dimensions of length, which can be associated to a mass parameter µby rµ= 2Gµ/c2. The determination of the function β(r) in the matter region leads to the determination of all the functions that describe both the matter and the geometry of the system within that region, by means of ρ(r) = rMβ′(r) κr2,(13) P(r) = KrMβ′(r) κr21+1/n ,(14) F(r) = KrMβ′(r) κr21/n ,(15) λ(r) = 1 2lnr r−rMβ(r),(16) ν(r) = ν(r2)−(n+ 1) ln[1 + F(r)].(17) This also determines the solutions within the inner and outer vacuum regions, which depend only on rµand rM, through the interface boundary conditions at r1and r2. The four free parameters of the system are K,nand 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 to be related to the matter, since it determines the general scale of the matter energy density, as can be seen from Equation (7). For all sets of parameters for which there is a solution of the differential problem the function β′(r) has a single point of maximum within the matter region, which is the point where β(r) has its single inflection point. For all existing solutions with r1>0 it holds that rµ>0, which implies that β(r) has a single zero within the matter region. The strictly positive value of rµimplies that the solutions have singularities at the origin. However, although the singularity has the invariant character of a curvature singularity, it has also a repulsive rather than attractive character, and thus is not associated to an infinite concentration of matter, but rather to exactly zero matter energy density around that point. Both for the subsequent analysis and for the numerical approach, it is convenient to further 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 r0>0. This mathematical device allows us to define a dimensionless radial variable and a dimensionless parameter associated to the mass Mby ξ=r r0 ,(18) ξM=rM r0 ,(19) as well as to define the dimensionless function of ξ, meant to assume the role of β(r), γ(ξ) = ξMβ(r).(20) 5 0.0 0.1 0.2 0.3 0.4 0.5 0 1 2 3 4 5 The DEC Curve The Tooper Curve The Limit Curve A Black Hole Limit The (C, πe) Parameter Plane for n= 1.5 πe C The DEC Curve The Tooper Curve The Limit Curve Black-Hole Limit Figure 1: The (C, πe) parameter plane for the case n= 1.5. The Tooper curve approaches the origin as an integer power, namely as C3. The limiting curve is given by πe= 1/C1.5. The functions γ(ξ) and β(r) are both determined by the second-order ordinary differential equation π′(ξ) = π(ξ)2 ξ−n n+ 1 1 ξ−γ(ξ) 1 + F(ξ, π) 2F(ξ, π)γ(ξ) ξ+F(ξ, π)π(ξ),(21) where π(ξ) = γ′(ξ) is the derivative of γ(ξ), the primes now indicate derivatives with respect to ξ, and F(ξ, π) is given by F(ξ, π) = Cπ(ξ) ξ21/n ,(22) where C=K/ κr2 01/n is a dimensionless constant associated to the polytropic constant K. This differential system can be interpreted as a pair of first-order coupled ordinary differential equations determining γ(ξ) and π(ξ), the other equation being simply γ′(ξ) = π(ξ).(23) This pair of first-order ordinary differential equations can be used for the numerical integration of this differential system, in order to obtain γ(ξ) and β(r). In terms of these dimensionless variables the equation of state shown in Equation (8) can be written as ¯ P(ξ) = C[¯ρ(ξ)]1+1/n ,(24) where we define the dimensionless pressure and energy density by 6 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0 1 2 3 4 5 6 7 8 The DEC Curve The Tooper Curve The Limit Curve A Black Hole Limit The (C, πe) Parameter Plane for n= 2.0 πe C The DEC Curve The Tooper Curve The Limit Curve Black Hole Limit Figure 2: The (C, πe) parameter plane for the case n= 2.0. The Tooper curve approaches the origin as an integer power, namely as C2. The limiting curve is given by πe= 1/C2.0. ¯ P(ξ) = κr2 0P(ξ),(25) ¯ρ(ξ) = κr2 0ρ(ξ).(26) Relation in Equation (7), when written in terms of the dimensionless variables, becomes π(ξ) = ξ2¯ρ(ξ).(27) Note that the mass Mand the Schwarzschild radius rMno longer appear explicitly in the equations. The mass Mis a free parameter of the system that has effectively been factored out of it. Later on, in Section 5, we will see that this simply means that there are solutions for all possible values of Mand rM. 3 The Parameter Planes We will begin by simply presenting the results of an extensive numerical exploration of the solutions, and subsequently we will discuss each element that is involved in this exploration. In their dimensionless formulation, the solutions are described by the three dimensionless free parameters n,Cand πe=π(ξe), where ξeis the position of the inflection point of γ(ξ), which is also the position of the maximum πeof π(ξ). For each value of n, the solutions can be classified according to a parameter plane spanned by Cand πe. We used a small discrete set of values of n, namely the values 1.5, 2.0, 2.5 and 3.0. For each one of these values we scanned the parameter planes described by the Cartesian coordinates (C, πe), 7 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 2 4 6 8 10 The DEC Curve The Tooper Curve The Limit Curve A Black Hole Limit The (C, πe) Parameter Plane for n= 2.5 πe C The DEC Curve The Tooper Curve The Limit Curve Black Hole Limit Figure 3: The (C, πe) parameter plane for the case n= 2.5. The Tooper curve approaches the origin as a fractional power, namely as C5/3. The limiting curve is given by πe= 1/C2.5. aiming at locating the allowed regions of these parameter planes, that is, the regions where the solutions, first of all exist, and second correspond to physically acceptable matter. In Figures 1 through 4 we present the allowed regions of the parameters planes, for the aforementioned four values of n. As one can see, the four parameter planes shown are generically quite similar. Each point in one of these graphs corresponds to a run of the numerical program [4] with the given values of n,Cand πe, as well as with the value ξe= 1 for the value of the position of the inflection point of γ(ξ), which is the position of the extremum πeof π(ξ). This choice, which since we have that ξe=re/r0can be written as re=r0, means only that from now on we choose the arbitrary parameter r0to be the radial position reof this inflection point. The runs were done on a grid with variations of Cgiven by ∆C= 0.005, starting from C= 0.005. Here is the part of the description of the parameter planes which is common to all four cases. The strong solid lines represent the Tooper curves. On these curves lie the solutions found by Tooper [5] a long time ago. Over these curves we have that ξ1= 0, meaning that these solution are for filled spheres rather than for shells, and also that ξµ= 0, which is the particular condition imposed by Tooper in order to obtain his solutions, a condition which avoids the singularity at the origin. Below these curves there are no acceptable solutions of the differential equation, since the solutions diverge when one integrates towards ξ= 0. Later on, in Subsection 3.4, we will examine in more detail exactly what happens in this case. Therefore, all existing solutions are necessarily above the Tooper curves. The strong dashed lines represent the Dominant Energy Condition [6] curves (or DEC curves for short). Above these curves the stress-energy tensor of the matter does not satisfy 8 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 2 4 6 8 10 The DEC Curve The Tooper Curve The Limit Curve A Black Hole Limit The (C, πe) Parameter Plane for n= 3.0 πe C The DEC Curve The Tooper Curve The Limit Curve Black Hole Limits Figure 4: The (C, πe) parameter plane for the case n= 3.0. The Tooper curve approaches the origin as a fractional power, namely as C3/2. The limiting curve is given by πe= 1/C3.0. The vertical lines mark the anomalous region. The crosses mark points where a divergence was detected when integrating outward, towards infinity. this condition, so that the behavior of the matter is not physically acceptable. Note that this does not mean that the solution of the differential equation itself does not exist above these curves, since in general it does in fact exist, and can be calculated without too much trouble. It just means that there is no physically possible matter that would result in these solutions. It is verified numerically that the DEC curves are always below the limiting curves shown in the graphs, with the light doted lines, which are given by πe(C) = 1 Cn.(28) We will see later, between the end of Subsection 3.1 and Subsection 3.2, that the DEC curve has this same type of dependency with C. Note that this indicates that the allowed regions extend indefinitely to arbitrarily large values of πe, always below the DEC curves, and becoming ever closer to C= 0. Of course, due to this we can only show a finite portion of these allowed regions in the graphs. In short, all physically acceptable solutions are necessarily below the DEC curves, as well as above the Tooper curves. It should be noted that over the C= 0 axis there are also no solutions, since in this case, according to the definition of Cgiven immediately after Equation (22), we have that K= 0, so that according to the equation of state shown in Equation (8) we then have pressureless dust, and therefore there are no stable static solutions. Therefore all acceptable solutions must have C > 0, that is, strictly positive C. It should also be noted that, since πeis the maximum value of the necessarily positive function π(ξ), the πe= 0 axis corresponds to 9 0.0 1.0 2.0 3.0 4.0 -0.5 0.0 0.5 1.0 gamma pi pi’ nu lambda rho A Run Below the Tooper Curve. ξ γ(ξ) π(ξ) π′(ξ) ν(ξ) λ(ξ) ¯ρ(ξ) Figure 6: All the functions obtained from a trial solution of the differential equation with parameters below the Tooper curve. The integration diverges when one integrates towards the origin. The left vertical line marks the initial point of the integration to either side, and the right vertical line marks the position of ξ2as obtained by the integration to the right side. The parameters used were n= 1.5, C= 1/3, ξe= 1.0 and πe= 0.5. to its previous value, decrease the variation ∆πeby some constant factor, and continue the downward search for the position of the DEC curve. This results in a fast exponential search for the position of the DEC curve, since the variation ∆πe, which represents the distance between the position of the search and the position of the DEC curve, decreases exponentially fast. Note that we can start with a rather large value of ∆πe, since it will be decreased exponentially fast. We chose to start with ∆πe= 0.1. The stop criterion for this downward search is the condition that the addition of ∆πeto πeno longer change πe, meaning that the relative precision represented by ∆πe/πefalls below the precision level at which the search program runs. Since that program runs in double precision mode (64-bit word), this guarantees the localization of the DEC curve within that double precision level. Naturally the precision of the results of the search is also affected by the precision with which the differential system is integrated. The integration code runs in quadruple precision mode (128-bit word), and we chose to use an integration interval of ∆ξ= 10−5. Since the integration program uses the Runge-Kutta fourth-order integration algorithm, this results in an excellent level of precision, certainly much greater than what is needed just for the graphs. Once we have obtained the position of the DEC curve for the current value of C, by means of the inner loop over values of πe, we enter an outer loop in which we now change the values of C. We increment Cby a fixed variation ∆C, which will characterize the grid 16 size used for the graphs. We chose to use ∆C= 0.005. Having the new value of C, we repeat the whole search process downward, for the position of the DEC curve. In this next search process we use as the initial value of πethe last value of this variable obtained in the previous search process. Due to the properties of the DEC curve, which has negative derivative everywhere, this starting value of πeis guaranteed to be above the DEC curve at the new value of C, but it is much closer to it than an initial value based on the limiting curve, thus making the search more efficient. The stop criterion for the loop over values of Cis that, when the current exponential search downward over πeconcludes, the integration program detect, in its last run, the presence of a divergence. This detection will also result in the creation of a flag file recording the event, to be used by external programs. The presence of a divergence means that the downward search over πenot only found the DEC curve, but also crossed below the Tooper curve at that value of C. Once the Tooper curve has been crossed or the DEC curve has become sufficiently close to the Tooper curve, there is no longer any chance of finding any point (C, πe) of the parameter plane that is both below the DEC curve and above the Tooper curve, within the precision limitations with which the graphs are being produced. 3.4 Determination of the Tooper Curve The determination of the position of the Tooper curve relies on the detection of divergences in the integration process, when one crosses below that curve. We should therefore examine in more detail, before anything else, what happens when we try to run the program below the Tooper curve. The program would usually stop cold when a divergence is detected, but there is an option to still run it for a little while, in order to plot the graph shown in Figure 6, where one can see a representation of the behavior of the differential system when we run below the Tooper curve, in a simple case, with the parameters n= 1.5, C= 1/3 and πe= 0.5. The behavior shown is quite typical. The graph shows all the relevant quantities as functions of the radial variable ξ. The left vertical dotted line, located at ξ= 1, is the location at which we start the integration of the differential equation shown in Equation (21), to either side. The integration to the right-hand side proceeds normally and produces a value ξ2≈2.9 for the outer radius, as indicated by the right vertical dotted line. To the right of this point we have the exterior Schwarzschild solution, with ξM≈0.97. However, the integration to the left-hand side cannot be completed, a value for ξ1is never found, and the process has to be interrupted when divergences are detected. As one can see in Figure 6, the function γ(ξ) never crosses the ξaxis, and never develops a constant value around the origin. Its derivative π(ξ) does not assume a second zero, and instead turns around and diverges quickly to positive infinity. The second derivative π′(ξ) diverges even faster to negative infinity. The graphs of the metric functions λ(ξ) and ν(ξ) do not cross one another, and in fact diverge to positive and negative infinity respectively, which is the exact opposite of their usual behavior, when the two graphs do cross, and then diverge in the opposite directions. The dimensionless matter energy density ¯ρ(ξ) displays a strong divergence to infinity, and therefore so does ρ(r), thus characterizing a hard singularity of the matter energy density at the origin. In this simple case it can be determined that π(ξ) diverges to infinity faster than 1/ξ, and therefore that neither π(ξ) nor ¯ρ(ξ) are integrable functions around the origin. However, due to the singular behavior of the measure factor exp[λ(ξ) + ν(ξ)] near the origin, it is unclear whether or not this divergence corresponds to either a finite or an infinite total amount of energy around the origin. The detection of the divergence included in the integration program consists of detecting the turnaround of the function π(ξ), and therefore the change 17 of sign of its derivative π′(ξ). 3.5 Algorithm for the Tooper Curve Let us now describe in detail how we determined the position of the Tooper curve. This was done once the DEC curve had already been determined. The exponential search downward towards the Tooper curve proceeded exactly like that for the DEC curve. For each value of nthe determination of the Tooper curve starts at the largest value of Cfor which the position of the DEC curve was successfully determined. At that point, for that value of C, we start the downward search at the value of πegiven by the DEC curve, which is always above the Tooper curve. Since we are looking for a divergence when integrating inward, we integrate the differential system to the left only, towards the origin. If the integration concludes with no divergence, resulting therefore in a definite value for ξ1, we decrease the value of πeby a certain variation ∆πe, and repeat the integration. This proceeds until a divergence is detected, meaning that we are now below the Tooper curve. Once this happens we abort the integration, return πeto its previous value, decrease the variation ∆πeby some constant factor, and continues the downward search for the position of the Tooper curve. Once again this results in a fast exponential search, this time for the position of the Tooper curve, since the variation ∆πe, which represents the distance between the position of the search and the position of the Tooper curve, decreases exponentially fast. Note once again that we can start with a rather large value of ∆πe, since it will be decreased exponentially fast, and again we chose to start with ∆πe= 0.1. Once more the stop criterion for this downward search is the condition that the addition of ∆πeto πeno longer change πe, meaning that the relative precision represented by ∆πe/πefalls below the precision level at which the search program runs. Since this second search program also runs in double precision mode (64-bit word), this guarantees the localization of the Tooper curve within that double precision level. Naturally, once again the precision of the results of the search is also affected by the precision with which the differential system is integrated. We recall that the integration code runs in quadruple precision mode (128-bit word), and again we chose to use an integration interval of ∆ξ= 10−5. Since the integration program uses the Runge-Kutta fourth-order integration algorithm, once again this results in an excellent level of precision, certainly much greater than what is needed just for the graphs. Once we have obtained the position of the Tooper curve for the current value of C, by means of the inner loop over values of πe, we enter an outer loop in which we now change the values of C, which we will now decrease rather than increase. We decrease Cby a fixed variation ∆C, which will characterize the grid size used for the graphs. Again we chose to use ∆C= 0.005. Having the new value of C, we repeat the whole search process downward, for the position of the Tooper curve. In this next search process we use as the initial value of πethe last value of this variable obtained in the previous search process. Due to the properties of the Tooper curve, which has positive derivative everywhere, this starting value of πeis guaranteed to be above the Tooper curve at the new value of C, but it is much closer to it than an initial value based on the DEC curve, again making the search more efficient. The stop criterion for the loop over values of Cis that Cwould become zero or negative at the next decrement. This can be done by checking whether the current value of Cis such that C≤∆C. We therefore stop the process of determining the Tooper curve at C= ∆C > 0. It should be noted that the stop criterion for the loop over πemay fail in some cases, if we are very close to C= 0. This is so because, specially for the case n= 1.5, the Tooper 18 curve approaches C= 0 with zero derivative. This tends to make πeand ∆πealways commensurate in this case, so that the test that the addition of ∆πeto πefails to change πemay never succeed. Due to this it was necessary to implement a maximum number of iterations of the loop over πe, which we chose to be 52. The average number of iterations in order to locate the Tooper curve within the double precision level is a little over 32. 4 Mapping the Observables On a second phase of the exploration we measured certain observables over the allowed regions, in order to characterize the general behavior of the solutions throughout the regions. These can be considered as diagnostic observables, and the results obtained for them can be seen in Figures 7 through 18. Using them we may, for example, classify the solutions as corresponding to either low-density or high-density distributions of matter, or as being either close to or distant from the formation of an event horizon. Each successful run of the program produces four fundamental numbers as results that characterize that solution, given by the dimensionfull radii r1,rµ,r2and rM, and hence by the corresponding dimensionless quantities ξ1,ξµ,ξ2and ξM. These are all encoded in the function γ(ξ); −ξµis the constant value of γ(ξ) within the inner vacuum region, that is, to the left of ξ1;ξ1is the point where π(ξ) = γ′(ξ) becomes zero when integrating inwards; ξMis the constant value of γ(ξ) within the outer vacuum region, that is, to the right of ξ2; and ξ2is the point where π(ξ) = γ′(ξ) becomes zero when integrating outwards. All these quantities have well-defined physical meanings, although we believe that in the case of rµthe complete physical meaning has not yet been completely established; r1 is the inner radius, or the radial coordinate characterizing the inner vacuum region; rµ, at least for the time being, has only the meaning of being zero for the Tooper solutions; r2 is the outer radius, or the radial coordinate characterizing the overall size of the matter distribution; rMis the Schwarzschild radius, associated to the total mass Mand to the total energy Mc2of the bound gravitational system. Next we will define our diagnostic observables and describe their main properties. 4.1 Definition of the Observables We will define, for the purposes of the calculations with the program, three ratios involving the quantities ξ1,ξµ,ξ2and ξM, that will have useful properties regarding the characterization of the solutions. All these ratios will assume values within the interval [0,1]. The first one we will name the Energy Ratio, which is defined as ER=ξµ ξM+ξµ .(59) This quantity has the property that it is zero when ξµ= 0, and it thus characterizes the Tooper curves. It also has the property that if ξµis very large, with ξµ≫ξM, then it approaches the value 1. This ER→1 limit characterizes, therefore, solutions which are maximally different from the Tooper solutions. Presumably these are solutions in which other forms of energy are more important than the energy associated to the total asymptotic gravitational mass M, such as, for example, the gravitational binding energy. The second ratio we will name the Horizon Ratio, defined as HR=ξM ξ2 .(60) 19 ✵✵ ✵✁ ✵✂ ✵✄ ✵☎ ✵✆ ✵✵ ✁✵ ✂ ✵ ✄✵ ☎✵ ✆✵ ✵✵ ✵✂ ✵☎ ✵✝ ✵✞ ✁✵ The ERratio for n= 1.5. ER Cπe Figure 7: The energy ratio ERfor n= 1.5. The contours are spaced vertically by 0.06. ✵✵ ✵✁ ✵✂ ✵✄ ✵☎ ✵✆ ✵✵ ✁✵ ✂ ✵ ✄✵ ☎✵ ✆✵ ✵✵ ✵✂ ✵☎ ✵✝ ✵✞ ✁✵ The HRratio for n= 1.5. HR Cπe Figure 8: The horizon ratio HRfor n= 1.5. The contours are spaced vertically by 0.06. ✵✵ ✵✁ ✵✂ ✵✄ ✵☎ ✵✆ ✵✵ ✁✵ ✂ ✵ ✄✵ ☎✵ ✆✵ ✵✵ ✵✂ ✵☎ ✵✝ ✵✞ ✁✵ The QRratio for n= 1.5. QR Cπe Figure 9: The quantum ratio QRfor n= 1.5. The contours are spaced vertically by 0.06. 20 Since all our solutions here have the property that ξM< ξ2, this quantity is always smaller than 1. However, there may be limits in which it approaches 1. This quantity has the property that for very low-density objects, for which ξM≪ξ2, is approaches zero. Therefore, this observable indicates low-density solutions when it is small, as is the case near the Tooper curves for small values of C. On the other hand, if ξMand ξ2become very close, then this ratio approaches 1. In this case one approaches a situation in which an event horizon forms at the position ξ2, according to the exterior Schwarzschild solution, which is valid outside ξ2, and which in this limit acquires the well-known coordinate singularity at ξ2. Therefore this ratio discriminates between low-density and high-density solutions, including the onset of black hole solutions, in the HR→1 limit. The third ratio we will name the Quantum Ratio, defined as QR=ξ1 ξ2 .(61) Clearly the name chosen requires some explanation. This has the property that it is zero if ξ1= 0, in which case we have solutions for filled spheres, instead of shells. In this way, it characterizes the Tooper curves just like ER, but since there are many other solutions for which ξ1≪ξ2, even if ξ16= 0, this diagnostic is not particularly helpful. On the other hand, if we recall that all the matter in the system is located between ξ1 and ξ2, we see that when the inner radius approaches the outer radius the proper volume containing the matter shrinks to zero, resulting in infinite localized matter energy densities. If we consider the dimensionless radial coordinate ξ, we have that in this case the coordinate thickness ξ2−ξ1of the shell goes to zero, so that the corresponding proper length will eventually become commensurate with the correlation lengths (or wavelengths) of the fields associated to the particles of matter within the shell. In this case it is to be expected that the quantum behavior of the matter will prevent this thickness from decreasing any further. In other words, it is reasonable to expect that the quantum properties of the matter will prevent the shell from degenerating into a true two-dimensional surface. Therefore, the limit QR→1 indicates situations in which the quantum properties of matter come into play. 4.2 Behavior of the Observables Graphs of the observables ER,HRand QRcan be seen in Figures 7 through 18, for all the values of nused. In these graphs we chose to use a step for Cgiven by ∆C= 0.01, and a step for πegiven by ∆πe= 0.1. As one can see in the graphs, the qualitative behavior of the ratios ER,HRand QRas functions of Cand πeis quite similar for all the values of nexplored, thus indicating that the same is at least qualitatively true for ξ1,ξµ,ξ2and ξMas well. In all the parameter plane graphs shown in Figures 1 through 4 there are essentially the same regions with the same specific properties for the solutions, with the possible exception of the rather small anomalous region for the larger values of C, in the case n= 3.0. This particular small region is characterized by very large values of ξ2, and hence by very small values of HRand QR. In this region the program takes a very long time to find the value of ξ2, and thus sometimes overtaxes the available computer resources. The integration interval used was ∆ξ= 10−5in all cases except the case n= 3.0, in which we found it necessary to use ∆ξ= 10−6. Let us start by describing how each one of the three ratios behaves along the allowed regions. The ratio ERgoes to zero at the Tooper curves, as expected, thus pinpointing the locations of these curves. It seems to increase towards 1 as we enter the asymptotic parts of 21 ✵✵ ✵✁ ✵✂ ✵ ✄ ✵☎ ✵ ✆ ✵ ✝ ✵✵ ✁✵ ✂✵ ✄✵ ☎ ✵ ✆✵ ✝✵ ✼✵ ✽✵ ✵✵ ✵✂ ✵☎ ✵ ✝ ✵✽ ✁✵ The ERratio for n= 2.0. ER Cπe Figure 10: The energy ratio ERfor n= 2.0. The contours are spaced vertically by 0.06. ✵✵ ✵✁ ✵✂ ✵ ✄ ✵☎ ✵ ✆ ✵ ✝ ✵✵ ✁✵ ✂✵ ✄✵ ☎ ✵ ✆✵ ✝✵ ✼✵ ✽✵ ✵✵ ✵✂ ✵☎ ✵ ✝ ✵✽ ✁✵ The HRratio for n= 2.0. HR Cπe Figure 11: The horizon ratio HRfor n= 2.0. The contours are spaced vertically by 0.06. ✵✵ ✵✁ ✵✂ ✵ ✄ ✵☎ ✵ ✆ ✵ ✝ ✵✵ ✁✵ ✂✵ ✄✵ ☎ ✵ ✆✵ ✝✵ ✼✵ ✽✵ ✵✵ ✵✂ ✵☎ ✵ ✝ ✵✽ ✁✵ The QRratio for n= 2.0. QR Cπe Figure 12: The quantum ratio QRfor n= 2.0. The contours are spaced vertically by 0.06. 22 the allowed regions, with C→0 and πe→ ∞, thus indicating that ξµbecomes very large as we travel deeper into these regions. The ratio HRhas values smaller than 1, but sometimes not close to zero, over the whole length of the Tooper curves, and very small values near the point C= 0. This means that there may be solutions with significant overall matter density on the Tooper curve, specially in the case N= 1.5. There are also solutions with very small such density near C= 0. This ratio also seems to increase towards 1 as we enter the asymptotic parts of the allowed regions, with C→0 and πe→ ∞, thus indicating that as we travel deeper into these regions we proceed towards high overall densities and the possible formation of event horizons. The ratio QRis zero over the whole length of the Tooper curves, thus indicating that we have ξ1= 0 there, as expected. It is also either zero or very small over most of the allowed regions, thus indicating that in general we have ξ1≪ξ2. It only becomes significantly larger than zero for small C, near the πeaxis. Also, it can be seen that it increases slowly as we travel into the asymptotic regions. This indicates that on limits for which C→0 and πe→ ∞ within these regions the quantum effects of the matter may eventually come into play. 4.3 The General Picture Examining what happens in each part of the allowed regions, some general conclusions can be formulated. To begin with, over or near the Tooper curves we never approach the formation of event horizons. The same is true for the parts of the DEC curves shown in the graphs. Another way to say this is to point out that for finite and non-zero values of C the solutions never form event horizons, and that for the larger values of Cthey never get anywhere even close to that. For each value of nthe asymptotic region, where C→0 and πe→ ∞, is the one region where we can approach the formation of an event horizon. This is also the region where we get extremely high overall matter densities, since this density increases as we enter more and more deeply into these asymptotic regions. The detailed exploration of the black hole limits in these asymptotic regions would require a different numerical approach, with a whole new set of runs, and due to this will have to be postponed to a subsequent paper. The two basic properties of the Tooper curves, that over them we have ξ1= 0 and ξµ= 0, are confirmed by our exploration. Very low density objects seem to exist mostly near the origin, for very small values of both Cand πe. That seems to be the realm of the main-sequence stars. It may be possible to find somewhat similar solutions in the small mostly regular region located beyond the anomalous region of the case n= 3.0, but at this point it is not really clear what the true character of these particular solutions might be. Below the Tooper curves, that is for very small values of πe, given some values of nand C, as well as for sufficiently large C, there are no solutions at all. The upper limit for Cis somewhere between 0.5 and 0.8 depending on the value of n. Close to the DEC curves we have extremely relativistic matter, with overall matter densities that vary from significant to very high. Since the matter involved in this solution is a simple ideal gas of massive particles, there is only one way in which its behavior can violate the dynamic precepts of Relativity, namely the speed of its particles either reaching or surpassing the speed of light. Therefore, it is to be expected that in this region the speed of sound in the gas is close to the speed of light, and also that their average speed of thermal agitation is close to the speed of light. This implies truly immense thermal energies and temperatures, and it can be imagined how this may also help to prevents the formation 23 ✵✵ ✵✁ ✵✂ ✵ ✄ ✵☎ ✵ ✆ ✵✝ ✵✞ ✵✵ ✂✵ ☎✵ ✝✵ ✽✵ ✁✵✵ ✵✵ ✵✂ ✵☎ ✵✝ ✵✽ ✁✵ The ERratio for n= 2.5. ER Cπe Figure 13: The energy ratio ERfor n= 2.5. The contours are spaced vertically by 0.06. ✵✵ ✵✁ ✵✂ ✵ ✄ ✵☎ ✵ ✆ ✵✝ ✵✞ ✵✵ ✂✵ ☎✵ ✝✵ ✽✵ ✁✵✵ ✵✵ ✵✂ ✵☎ ✵✝ ✵✽ ✁✵ The HRratio for n= 2.5. HR Cπe Figure 14: The horizon ratio HRfor n= 2.5. The contours are spaced vertically by 0.06. ✵✵ ✵✁ ✵✂ ✵ ✄ ✵☎ ✵ ✆ ✵✝ ✵✞ ✵✵ ✂✵ ☎✵ ✝✵ ✽✵ ✁✵✵ ✵✵ ✵✂ ✵☎ ✵✝ ✵✽ ✁✵ The QRratio for n= 2.5. QR Cπe Figure 15: The quantum ratio QRfor n= 2.5. The contours are spaced vertically by 0.06. 24 of event horizons. It is important to emphasize here that these are relativistically high temperatures, orders of magnitude above what would be considered as high temperatures in the realm of stars and other celestial objects, even for the hottest such objects. Therefore, close enough to the DEC curves these cases seem to be far beyond any objects that one could ever hope to find in nature. 5 Application to Astrophysics Let us now discuss how one can use this family of solutions in order to obtain the solutions for definite objects, say with given mass and outer radius. To start with, let us assume that we want the solution for an object of mass M, and therefore with Schwarzschild radius rM. The problem poses itself of what values of the free parameters n,Cand πeshould be used in order to obtain a solution with the required properties. Note that the value of Mis not directly involved in the choice of these dimensionless parameters. We will see that there is no single solution to this first problem, but rather a significant amount of arbitrary choice involved. Regarding the value of n, one can either choose an appropriate value of nto represent the type of thermodynamic behavior of the matter involved, or one can construct, for example, a two-layer solution, with layers that are connected by means of appropriate interface boundary conditions at a certain interface between two different layers of matter, with different thermodynamic behaviors. A two-layer model of a star with a radiative inner layer, say with n= 3.0, and an isentropic outer layer, say with n= 1.5, comes to mind. Putting aside for now this more complex alternatives, let us just assume, for the purposes of this discussion, that we simply chose some value of nfor a single-layer solution. For any values of Cand πethat one chooses to run the program with, assuming only that the point (C, πe) is within the allowed region of the parameter plane for the value of n to be used, there will result some definite values for ξ1,ξµ,ξ2and ξM. If we have a required value of rM, then we can at once determine the parameter r0, for any values of Cand πe whatsoever, since we have that r0=rM ξM .(62) This determines the position r0of the inflection point of β(r), and therefore the position of the point of maximum of β′(r). Having determined r0, we can now determine all other relevant dimensionfull quantities, since we have that r1=ξ1r0,(63) rµ=ξµr0,(64) r2=ξ2r0.(65) We therefore see that there are solutions for all possible positive values of rM, and hence of M. The fact that the parameter rMwas factored out of the system, when we wrote it in terms of dimensionless parameters, betrays the existence of a scaling transformation among the solutions, parametrized by the value of rM. Given n, any point (C, πe) in the allowed region of the corresponding parameter plane corresponds in fact to an infinite set of similar solutions, for all possible values of rM, and hence of M. At this point we may also determine the polytropic constant K, since we have that K=Cκr2 01/n ,(66) 25 we are led to surmise that solutions representing black holes, or solutions that are in some sense close to them, are located. It was possible to identify families of curves within the allowed regions, that extend indefinitely within these asymptotic regions, and that have the property of approaching black-hole limits. Of course these infinite asymptotic regions will have to be explored using numerical approaches other than mapping the complete regions, and thus such an effort will have to be postponed to a subsequent paper, to be dedicated to this issue. Besides these asymptotic regions, there are other small regions which might deserve further scrutiny, as was pointed out in the text. And we should not forget that the cases with n < 1, which promise to be significantly different from the ones explored here, are still to be examined in any significant amount of detail. In fact, there are many possible extensions and generalizations of the exploration presented here, as well as different studies that could be done with these solutions. Of particular interest would be the study of the temperature in these solutions, specially the observable surface temperature. However, in all fairness we would like to finish by pointing out that perhaps the most significant limitation of this family of solutions is a rather fundamental one, namely that they do not include the angular momentum of rotating objects. The extension of these solutions to that case remains as a very significant and important open challenge. Acknowledgments The author would like to thank Prof. Oscar J. P. ´ Eboli for the use of the computer systems of his research group. Data Availability Statement Data sharing not applicable to this article as no experimental or observational datasets were generated or analyzed during the current study. Conflict of Interest Statement The author hereby certifies that there are no actual or potential conflicts of interest in relation to this article. References [1] J. L. deLyra and C. E. I. Carneiro, “Complete solution of the einstein field equations for a spherical distribution of polytropic matter,” General Relativity and Gravitation, vol. 55, 2023. Article ID: 67; GRG DOI: 10.1007/s10714-023-03115-6; Zenodo DOI: 10.5281/zenodo.5087722. [2] J. L. deLyra, R. de A. Orselli, and C. E. I. Carneiro, “Exact solution of the einstein field equations for a spherical shell of fluid matter,” General Relativity and Gravitation, vol. 55, 2023. Article ID: 68; GRG DOI: 10.1007/s10714-023-03116-5; Zenodo DOI: 10.5281/zenodo.5087611. [3] P. A. M. Dirac, General Theory of Relativity. John Wiley & Sons, Inc., 1975. ISBN 0-471-21575-9. 32 [4] J. L. deLyra, “Program for spherically symmetric and static relativistic polytropes.” DOI/Zenodo. Freely downloadable compressed tar file with Fortran source code, available at the URL https://zenodo.org/records/8360689 or at the URL https://doi.org/10.5281/zenodo.8360689. [5] R. F. Tooper, “General relativistic polytropic fluid spheres,” Astrophys. J., vol. 140, pp. 434–459, 1964. [6] R. Wald, General Relativity. University of Chicago Press, 2010. [7] J. Ni, “Solutions without a maximum mass limit of the general relativistic field equations for neutron stars,” Science China, vol. 54, no. 7, pp. 1304–1308, 2011. [8] L. Nesluˇsan, “The ni’s solution for neutron star and outward oriented gravitational attraction in its interior,” Journal of Modern Physics, vol. 6, pp. 2164–2183, 2015. [9] S. Weinberg, Gravitation and Cosmology. New York: John Wiley and Sons, 1972. 33