New Solutions of Tolman-Oppenheimer-Volkov-Equation and of Kerr Spacetime with Matter and the Corresponding Star Models

Jan Helm

Technical University of Berlin, Berlin, Germany.

DOI: 10.4236/jhepgc.2022.83052                

Abstract

The Tolman-Oppenheimer-Volkov (TOV) equation is solved with a new ansatz: the external boundary condition with mass M₀ and radius R₁ is dual to the internal boundary condition with density ρ_bc and inner radius r_i, and the two boundary conditions yield the same result. The inner boundary condition is imposed with a density ρ_bc and an inner radius r_i, which is zero for the compact neutron stars, but non-zero for the shell-stars: stellar shell-star and galactic (supermassive) shell-star. Parametric solutions are calculated for neutron stars, stellar shell-stars, and galactic shell-stars. From the results, an M-R-relation and mass limits for these star models can be extracted. A new method is found for solving the Einstein equations for Kerr space-time with matter (extended Kerr space-time), i.e. rotating matter distribution in its own gravitational field. Then numerical solutions are calculated for several astrophysical models: white dwarf, neutron star, stellar shell-star, and galactic shell-star. The results are that shell-star star models closely resemble the behaviour of abstract black holes, including the Bekenstein-Hawking entropy, but have finite redshifts and escape velocity v < c and no singularity.

Keywords

General Relativity, Tolman-Oppenheimer-Volkov Equation, Neutron Stars, Shell Stars

1. Introduction

In General Relativity, one of the most important applications is to calculate the mass distribution and the space-time metric for a given equation-of-state of a stellar model.

Without rotation, one has spherical symmetry and then the Tolman-Oppen-Heimer-Volkov (TOV) equation in radius r, which is derived directly from the Einstein equations (see [1]), is being used. The TOV equation consists originally of 2 coupled non-linear ordinary differential equations (odeq) of degree 1 in r for mass M(r) and density ρ(r), where $4\pi r^2\rho \left(r\right)=M\left(r\right)$, and can be transformed into one ordinary differential equation of degree 2 for M(r) by eliminating ρ(r).

The boundary condition is imposed normally at r = 0 with M(0) = 0 and $\rho \left(0\right)={\rho }_0$, where ρ₀ is the maximal density. Then the TOV-equation for M(r) is solved with this boundary condition at r = 0 for M(r) and M'(r), which gives the total mass M₀(ρ₀) and the total radius R(ρ₀), and a mass-radius relation M₀(R).

The predominant view of the neutron stars and stellar black-holes is, that neutron stars obey an equation-of state (eos) of an interacting-fluid model [2], which solutions of the TOV equation up to about M = 3M_sun. For larger masses, it is assumed that only a black-hole solution remains. This is based on the so-called Oppenheimer limit for the radius of a compact mass

$R_{\lim }=\frac{9}{8}{r}_s=\frac{9MG}{4c^2}$ : for a smaller compact object, the density at the center becomes infinite.

The new ansatz presented here is the extended (inner) boundary condition at $r=r_i$ with the non-zero inner radius r_i, $M\left(r_i\right)=0$ and $\rho \left(r_i\right)={\rho }_0$, i.e. the star becomes a shell-star with an (almost) void interior. With the parameters r_i and ρ₀, this ansatz generates a 2-parametric solution manifold, where, because of energy minimization, the stable physical solution is the one with minimal r_i for a given ρ₀, which determines the total mass $M_0\left(r_i,{\rho }_0\right)$ and the total radius $R\left(r_i,{\rho }_0\right)$. This ansatz circumvents the Oppenheimer limit, because the mass is non-compact, and the density at the center is zero. It yields valid solutions of the TOV-equation with a continuous mass-radius relation $M_0\left(R,r_i\right)$.

The Oppenheimer limit $R_{\lim }=\frac{9}{8}{r}_s=\frac{9MG}{4c^2}$ is deduced from the maximum pressure at r = 0:

$P\left(0\right)={\rho }_0{c}^2\frac{1-\sqrt{1-\frac{r_s}{R}}}{3\sqrt{1-\frac{r_s}{R}-1}}$

The dual (outer) boundary condition is the one at r = R with M = M₀ and $\rho ={\rho }_{bc}$, where ρ_bc depends on the equation-of-state (eos): for neutron stars with interacting nucleon fluid with the equilibrium nucleon density ρ_c, and for nucleon Fermi-gas (stellar shell-stars) ${\rho }_{bc}=0$.

The 2 parameters R and M₀ in the dual outer boundary condition correspond uniquely to the 2 parameters r_i and ρ₀ in the inner boundary condition.

With rotation, one has an axisymmetric model in the variables r and θ (azimuthal angle), and has to solve the Einstein equations in these 2 coordinates. In vacuum, the corresponding solution is the Kerr space-time in r and θ.

With mass, a good starting point is using the extended Kerr space-time in Boyer-Lindquist coordinates with correction-factor functions A0, …, A4 and B0, …, B4 and the mass M(r, θ) as variables and insert this into the Einstein equations.

Setting Bi = 0, 4 of the 10 Einstein equations become trivial and one is left with 6 partial-differential equations (pdeq) in r and θ for the 6 variables A0, …, A4 and M.

The (outer) boundary condition here at the effective star radius R with total mass M₀ is: Ai = 1, M = M₀ and ¶_rAi = 0, ¶_rM = 0, as the density becomes 0 and the space-time becomes the normal Kerr space-time in vacuum.

Now, with rotation, we have a new model parameter, the angular velocity ω, to which corresponds a third parameter in the outer boundary condition: (outer) ellipticity ΔR₁, where $R_{1x}=R_{1y}-\Delta R_1$ and R_1x and R_1y are the equatorial and the polar radius. As in the TOV-case, here to the 3 parameters R_1y, M₀ and ΔR₁ correspond the 3 inner parameters r_iy, ρ₀ and Δr_i.

So here we get a 3-parametric solution manifold, and as in the spherical case, for a given total mass M₀ we have to find the stable physical solution. As before, these will be the ones with minimal riy and among them the one with minimal mean energy density: this defines the inner ellipticity Δr_i. In all considered cases, it can be shown numerically, that such a (non-trivial) minimum exists.

The paper is organized as follows.

In 2 we present the mathematical setup, in 3 the equations for the extended Kerr space-time with rotation, in 4 the solution algorithm for it. In 5 the TOV-equation is introduced, in 6 the equation-of-state for the nucleon fluid and nucleon gas. In 7 the results for the TOV-equation are shown: the parametric solution manifold in 7.1. and the case study for typical stars in 7.2. In 8 the results for the extended Kerr space-time with rotation are presented for three typical star configurations: compact neutron star, stellar shell-star, and galactic shell-star.

2. The Kerr Space-Time, Schwarzschild Space-Time, Einstein Equations

Using the Minkowski metric ${\eta }_{\mu \nu }$ ${\eta }_{\mu \nu }=\text{diag}\left(\text{1},-\text{1},-\text{1},-\text{1}\right)$, the Kerr space-time metric in original Kerr coordinates (u, θ, ϕ) has the line element [1] [3]

$\begin{array}{c}\text{d}{s}^2=\left(1-\frac{r{r}_s}{r^2+a^2{\cos }^2\theta }\right){\left(\text{d}u+a{\sin }^2\theta \text{d}\phi \right)}^2\\ -2\left(\text{d}u+a{\sin }^2\theta \text{d}\phi \right)\left(\text{d}r+a{\sin }^2\theta \text{d}\phi \right)\\ -\left(r^2+a^2{\cos }^2\theta \right)\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)\end{array}$ (1)

where $r_s=\frac{2GM}{c^2}$ is the Schwarzschild radius, and $a=\frac{J}{Mc}$ is the angular momentum radius (amr), a has the dimension of a distance: $\left[a\right]=\left[r\right]$, and J is the angular momentum.

With this line element the Kerr metric tensor $g_{\mu \nu }$ is as follows: [3]

$g_{\mu \nu }=\left(\begin{array}{cccc}1-\frac{r{r}_s}{{\rho }_{12}}& -1& 0& -\frac{r{r}_{s}a{\sin }^2\theta }{{\rho }_{12}}\\ & 0& 0& -a{\sin }^2\theta \\ & & -{\rho }_{12}& 0\\ & & & g_{33}\end{array}\right)$ (2)

with the abbreviations ${\rho }_{12}=r^2+a^2{\cos }^2\theta$ and $g_{33}=-{\sin }^2\theta \left(r^2+a^2+\frac{r{r}_s{a}^2{\sin }^2\theta }{{\rho }_{12}}\right)$.

In the limit $a\to 0$ the Schwarzschild space-time in advanced Eddington-Finkelstein coordinates emerges:

$\text{d}{s}^2=\left(1-\frac{r_s}{r}\right)\text{d}{u}^2-2\text{d}u\text{d}r-r^2\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)$ (3)

This form of the Schwarzschild line element has the advantage in comparison with the original line element

$\text{d}{s}^2=\left(1-\frac{r_s}{r}\right)c^2\text{d}{t}^2-\frac{\text{d}{r}^2}{1-\frac{r_s}{r}}-r^2\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)$ (3a)

that the (apparent) singularity at $r=r_s$ is missing.

The same is valid for the original Kerr space-time: the denominator ${\rho }_{12}$ has no zeros, there is no singularity in $g_{ab}$, which makes it more well-behaved numerically.

Alternatively, in Boyer-Lindquist-coordinates: [3]

$g_{\mu \nu }=\left(\begin{array}{cccc}1-\frac{r{r}_s}{{\rho }_{12}}& 0& 0& \frac{r{r}_{s}a{\sin }^2\theta }{{\rho }_{12}}\\ & -{\rho }_{12}/{\Lambda }_{12}& 0& 0\\ & & -{\rho }_{12}& 0\\ & & & -{\sin }^2\theta \left(r^2+a^2+\frac{r{r}_s{a}^2{\sin }^2\theta }{{\rho }_{12}}\right)\end{array}\right)$ (4)

with the line element

$\begin{array}{c}\text{d}{s}^2=\left(1-\frac{r{r}_s}{r^2+a^2{\cos }^2\theta }\right){\left(\text{d}t\right)}^2+\left(\frac{2r{r}_{s}a{\sin }^2\theta }{r^2+a^2{\cos }^2\theta }\right)\text{d}t\text{d}\phi \\ -\left(\frac{r^2+a^2{\cos }^2\theta }{r^2-r{r}_s+a^2}\right)\text{d}{r}^2-\left(r^2+a^2+\frac{r{r}_s{a}^2{\sin }^2\theta }{r^2+a^2{\cos }^2\theta }\right){\sin }^2\theta \text{d}{\phi }^2\\ -\left(r^2+a^2{\cos }^2\theta \right)\left(\text{d}{\theta }^2\right)\end{array}$ (4a)

with the abbreviation ${\Lambda }_{12}=r^2-r{r}_s+a^2$. Here, Λ₁₂ has zeros at the inner/outer horizon $r=\left(r_s/2\right)\pm \sqrt{{\left(r_s/2\right)}^2-a^2{\cos }^2\theta }$, so for numerical calculations the singularity has to be removed by adding a small ε: ${\Lambda }_{12}=\sqrt{{\left(r^2-r{r}_s+a^2\right)}^2+{\epsilon }^2}$.

In the limit $a\to 0$ the Schwarzschild space-time in the standard form (4) emerges.

The Einstein field equations with the above Minkowski metric are:

$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }{R}_0-\Lambda g_{\mu \nu }=-\kappa T_{\mu \nu }$ (5)

where $R_{\mu \nu }$ is the Ricci tensor, R₀ the Ricci curvature, $\kappa =\frac{8\pi G}{c^4}$, $T_{\mu \nu }$ is the energy-momentum tensor, $\Lambda$ is the cosmological constant (in the following neglected, i.e. set 0),

with the Christoffel symbols (second kind)

${\Gamma }_{\mu \nu }^{\lambda }=\frac{1}{2}{g}^{\lambda \kappa }\left(\frac{\partial g_{\kappa \mu }}{\partial x^{\nu }}+\frac{\partial g_{\kappa \nu }}{\partial x^{\mu }}-\frac{\partial g_{\mu \nu }}{\partial x^{\kappa }}\right)$ (6)

and the Ricci tensor

$R_{\mu \nu }=\frac{\partial {\Gamma }_{\mu \rho }^{\rho }}{\partial x^{\nu }}-\frac{\partial {\Gamma }_{\mu \nu }^{\rho }}{\partial x^{\rho }}+{\Gamma }_{\mu \rho }^{\sigma }{\Gamma }_{\sigma \nu }^{\rho }-{\Gamma }_{\mu \nu }^{\sigma }{\Gamma }_{\sigma \rho }^{\rho }$ (7)

The crucial part of the extended Kerr solution is the expression for the energy-momentum tensor $T_{\mu \nu }$. As usual, one uses the formula for the perfect fluid [1, (45.3)]:

$T_{\mu \nu }=\left(\rho +\frac{P}{c^2}\right)u_{\mu }{u}_{\nu }-P{g}_{\mu \nu }$ (8)

where P and ρ is the pressure and density, $u_{\mu }$ is the covariant velocity 4-vector.

In the Schwarzschild case, when deriving the TOV-equation, one sets the spatial contravariant velocity components to 0: $u^i=0$, in the Kerr case the tangential velocity $u^3=u^{\phi }\ne 0$.

For the velocity one has:

$u^3=r\omega =r\frac{aMc}{I}$, where I is the moment of inertia and ω the angular velocity, $I=f_{I}M{R}_1^2$, so

$u^3=\frac{1}{f_I}\frac{acr}{R_1^2}$

for a homogeneous sphere $I=\frac{2}{5}M{R}^2$ i.e. $f_I=\frac{2}{5}$, for a thin shell $I=\frac{2}{3}M{R}^2$ i.e. $f_I=\frac{2}{3}$, for a disc $I=\frac{1}{2}M{R}^2$.

If we make the obvious assumption that the star rotates as a whole, i.e. with constant angular velocity, then the moment of inertia I becomes r-dependent, like the mass M:

$M\left(r\right)={\displaystyle {\int }_0^r\rho \left(r_1\right)3r_1^2\text{d}{r}_1}$ (9)

$I\left(r\right)={\displaystyle {\int }_0^r\rho \left(r_1\right)3r_1^4\text{d}{r}_1}$

The factor 3 in the integral instead of the usual 4π comes from the dimensionless calculation in “sun units” (see below).

The amr a also becomes r-dependent:

$a\left(r\right)=\frac{J}{Mc}=\frac{I\left(r\right)\omega }{M\left(r\right)c}$

In the relativistic axisymmetric case with rotation with angular velocity ω, u^μ has the form [4]: $u^{\mu }=u^0\left(1,0,0,w\right)$

$u^{\mu }=\left(u^0,0,0,\omega u^0\right)$

Now $u^0$ is calculated from the condition

$c^2=g_{\mu \nu }{u}^{\mu }{u}^{\nu }$

and the covariant velocity from

$u_{\mu }=g_{\mu \nu }{u}^{\nu }$

The resulting expression for $u^0$ is (A_i are the Kerr correction-factors, mass M₁[r₁], moment of inertia I₁[r₁]): (9a)

$\begin{array}{c}{u}^0=(A_0\left(r_1,\theta \right)\sqrt{{\left(\frac{r_1{M}_1{\left(r_1\right)}^3}{{\omega }^2{I}_1{\left(r_1\right)}^2{\cos }^2\left(\theta \right)+r_1^2{M}_1{\left(r_1\right)}^2}-1\right)}^2}\\ -\frac{r_1{\omega }^4{I}_1{\left(r_1\right)}^2{M}_1\left(r_1\right){\cos }^4\left(\theta \right)A_3\left(r_1,\theta \right)}{{\omega }^2{I}_1{\left(r_1\right)}^2{\cos }^2\left(\theta \right)+r_1^2{M}_1{\left(r_1\right)}^2}\\ -\frac{{\omega }^4{I}_1{\left(r_1\right)}^2{\sin }^2\left(\theta \right)A_3\left(r_1,\theta \right)}{M_1{\left(r_1\right)}^2}-r_1^2{\omega }^2{\sin }^2\left(\theta \right)A_3\left(r_1,\theta \right)\\ +{\frac{2{\omega }^2\sqrt{I_1{\left(r_1\right)}^2}{M}_1\left(r_1\right)\sqrt{M_1{\left(r_1\right)}^2}{\sin }^2\left(\theta \right)A_4\left(r_1,\theta \right)}{{\omega }^2{I}_1{\left(r_1\right)}^2{\cos }^2\left(\theta \right)+r_1^2{M}_1{\left(r_1\right)}^2})}^{-1}\end{array}$

The state equation for the pressure P for the nucleon gas has the form

$P=c_1{\rho }^{\gamma }$

or in the dimensionless form with a critical density ${\rho }_c$ and dimensionless pressure $P_1$ and density ${\rho }_1$

$P_1:=\frac{P}{c^2{\rho }_c}=k_1{\left(\frac{\rho }{{\rho }_c}\right)}^{\gamma }=k_1{\left({\rho }_1\right)}^{\gamma }$ (10)

For the horizon, with rotation there is the inner and the outer horizon (M = M₀)

$r_{-}=\frac{M_0}{2}-\sqrt{{\left(\frac{M_0}{2}\right)}^2-a^2}$

$r_{+}=\frac{M_0}{2}+\sqrt{{\left(\frac{M_0}{2}\right)}^2-a^2}$

3. The Equations for the Extended Kerr Space-Time

The solution process starts with the metric tensor $g_{\mu \nu }$ in original Eddington-Finkelstein-coordinates, with 6 non-zero components, corresponding correction-factor functions A0, …, A5, and additive correction functions B0, …, B3 for the zero components.

$g_{\mu \nu }=\left(\begin{array}{cccc}\left(1-\frac{r{r}_s}{{\rho }_{12}}\right)A0& -A3& B1& -\frac{r{r}_{s}a{\sin }^2\theta }{{\rho }_{12}}A4\\ & B0& B2& -a{\sin }^2\theta A5\\ & & -{\rho }_{12}A1& B3\\ & & & -A2{\sin }^2\theta \left(r^2+a^2+\frac{r{r}_s{a}^2{\sin }^2\theta }{{\rho }_{12}}\right)\end{array}\right)$ (11)

and alternatively in Boyer-Lindquist-coordinates, corresponding correction-factor functions A0, …, A4, and additive correction functions B0, …, B4 for the zero components

$g_{\mu \nu }=\left(\begin{array}{cccc}\left(1-\frac{r{r}_s}{{\rho }_{12}}\right)A0& B0& B1& \frac{r{r}_{s}a{\sin }^2\theta }{{\rho }_{12}}A4\\ & -A1{\rho }_{12}/{\Lambda }_{12}& B2& B3\\ & & -{\rho }_{12}A2& B4\\ & & & -A3{\sin }^2\theta \left(r^2+a^2+\frac{r{r}_s{a}^2{\sin }^2\theta }{{\rho }_{12}}\right)\end{array}\right)$ (12)

The equations are the 10 Einstein equations eqR00, eqR11, eqR22, eqR33, eqR12, eqR23, eqR31, eqR01, eqR02, eqR03 in the (dimensionless) variables relative radius $r_1=\frac{r}{r_{ss}}$ and complementary azimuth angle ${\theta }_1=\frac{\pi }{2}-\theta$ with energy tensor $T_{\mu \nu }$ from (8) and the state equation $P_1=k_1{\left({\rho }_1\right)}^{\gamma }$ for the relative pressure $P_1$ and the relative density ${\rho }_1$. We are using the so called “sun units” $r_{ss}=r_s\left(sun\right)$, $M_s=M\left(sun\right)$, ${\rho }_s=\frac{M_s}{4\pi r_{ss}^3}$, $P_s={\rho }_s{c}^2$ for radius r, mass M, density ρ, and pressure P, respectively.

In “sun units” the original angle differential $\text{d}\Omega =4\pi \sin \theta r^2\text{d}r\text{d}\theta$ is transformed into $\text{d}\Omega =3\cos \theta r^2\text{d}r\text{d}\theta$, as for $\theta =0,\cdots ,\pi /2$, $r=0,\cdots ,\text{1}$ : ${\displaystyle \int \text{d}\Omega }=1$.

Also, all equations and variables are symmetric (even) in θ: Ai(−θ) = Ai(θ).

From now on we skip the index of the dimensionless variables and use the original notation, e.g. r instead of r₁.

Furthermore, we adopt the Boyer-Lindquist coordinates and the metric tensor (12).

In sun units, the Boyer-Lindquist metric tensor becomes:

$g_{\mu \nu }=\left(\begin{array}{cccc}\left(1-\frac{r{M}_0}{{\rho }_{12}}\right)A0& B0& B1& \frac{r{M}_{0}a{\cos }^2\theta }{{\rho }_{12}}A4\\ & -A1{\rho }_{12}/{\Lambda }_{12}& B2& B3\\ & & -{\rho }_{12}A2& B4\\ & & & -A3{\cos }^2\theta \left(r^2+a^2+\frac{r{M}_0{a}^2{\cos }^2\theta }{{\rho }_{12}}\right)\end{array}\right)$ (12a)

${\rho }_{12}=r^2+a^2{\sin }^2\theta$

${\Lambda }_{12}=r^2-r{M}_0+a^2$

where M₀ is the mass in sun units.

The 10 Einstein equations have a distinctive structure: there are 6 primary variables A0, A2, A3, A4, B1, B4 with highest derivative ¶_rrand 4 secondary variables A1, B2, B0, B3 with highest derivative ¶_r. Primary variables have boundary conditions for the variable and itsr-derivative, secondary variables only for the variable itself.

This structure is dual in θ: again there are 6 θ-primary and 4 θ-secondary variables.

6 of Einstein equations contain only one 2-derivative ¶_rr of a primary variable: eqR03(¶_rrA4), eqR22(¶_rrA2), eqR00(¶_rrA0), eqR33(¶_rrA3), eqR02(¶_rrB1), eqR23 (¶_rrB4) 3 contain only 2-derivatives of a secondary variable: eqR12, eqR01, eqR31, eqR11 contains all derivatives ¶_rr of a primary variable.

If we make the ansatz Bi = 0, several of the eqRij become identically 0, and we get the 6 equations eqR00, eqR11, eqR22, eqR33, eqR03, eqR12 for the 6 variables Ai and ρ, with the highest derivatives resp. ¶_rrA0, ¶_θθA1, ¶_rrA2, ¶_rrA3, ¶_rrA4, (¶_rrA2, ¶_θθA1).

Thus, we are left with the 6 differential equations degree 2 in r, θ non-linear (quartic) in variables Ai and their 1-derivatives and linear in ρ, ρ^γ.

In total, we have 6 algebro-differential eqs for 6 variables Ai and ρ (ρ enters only algebraically).

We can add 2 dependent equations eqR41 == $D_{\mu }{T}^{\mu 1}=0$ and eqR42 == $D_{\mu }{T}^{\mu 2}=0$, from the covariant continuity equations $D_{\mu }{T}^{\mu \nu }=0$, where $D_{\lambda }{T}^{\mu \nu }={\partial }_{\lambda }{T}^{\mu \nu }+{\Gamma }_{\lambda \kappa }^{\mu }{T}^{\kappa \nu }+{\Gamma }_{\lambda \kappa }^{\nu }{T}^{\mu \kappa }$ is the gravitational covariant derivative.

In eqR41 ρ enters with ¶_rρ, in eqR42 ρ enters with ¶_θρ.

So, alternatively, we have the diff. equations eqR00, eqR11, eqR22, eqR33, eqR03, eqR41, with the highest derivatives resp. ¶_rrA0, ¶_θθA1, ¶_rrA2, ¶_rrA3, ¶_rrA4, ¶_rρ (diff. eq. degree 1 in r for ρ) or the diff. equations eqR00, eqR11, eqR22, eqR33, eqR03, eqR42, with the highest derivatives resp. ¶_rrA0, ¶_θθA1, ¶_rrA2, ¶_rrA3, ¶_rrA4, ¶_θρ (diff. eq. degree 1 in θ for ρ).

In the Schwarzschild spacetime ω = 0 and a = 0, we have spherical symmetry, no dependence on θ, and the TOV-equation can be derived from the non-trivial eqR00, eqR11, eqR22, eqR41.

We impose an r-θ-analytic boundary condition for Ai, ¶_rAi, at r = R₁(R₁ is the star radius):

Ai = 1, ¶_rA0 = 0, ¶_rA2 = 0, ¶_rA3 = 0, ¶_rA4 = 0. For A1, there is no differential boundary condition, as ¶_rA1 is the highest r-derivative, for ρ there is no boundary condition at all, because r is algebraic in the equations, but there is an integral condition:

$M\left(R_1\right)={\displaystyle \underset{0}{\overset{\pi /2}{\int }}{\displaystyle {\int }_0^{R_1}\rho \left(r_1,\theta \right)3r_1^2}\cos \theta \text{d}{r}_1\text{d}\theta }=M_0$ : integral(ρ) = total mass = M₀.

In order to avoid the clumsy integral condition for ρ, we can introduce the mass M as a variable:

$M\left(r,\theta \right)={\displaystyle {\int }_0^r\rho \left(r_1,\theta \right)3r_1^2\text{d}{r}_1}$

$M\left(r\right)={\displaystyle {\int }_0^{\pi /2}M\left(r,\theta \right)\text{d}\theta }={\displaystyle {\int }_0^{\pi /2}{\displaystyle {\int }_0^r\rho \left(r_1,\theta \right)3r_1^2}\text{d}{r}_1\cos \theta \text{d}\theta }$ is the mass of the sphere(r) and

${\partial }_{r}M\left(r,\theta \right)=\rho \left(r,\theta \right)3r^2$

For M(r, θ) we impose the boundary condition at r = R₁:

$M\left(R_{\text{1}},\theta \right)=M_0$, ${\partial }_{r}M\left(R_1,\theta \right)=0$ (i.e. density ρ is zero at boundary, and total mass M₀).

So, if we take the diff. equations eqR00, eqR11, eqR22, eqR33, eqR03, eqR41 and replace $\rho \left(r,\theta \right)$ by ${\partial }_{r}M\left(r,\theta \right)$, we have 6 diff.equations in r, θ of degree 2, for the variables Ai(r, θ) (metric correction factors = mcf) and the mass M(r, θ), with the highest derivative ¶_rrM(r, θ) in M.

According to the Cauchy-Kovalevskaya theorem there exists then a unique solution in a region $R_1>r>r_i$ within the boundary. Inside the region $r_i>r>0$ we can enforce the vacuum Kerr-spacetime with the trivial solution Ai = 1, ρ = 0, i.e. there is no matter there, r_i the inner radius.

The Cauchy-Kovalevskaya theorem guarantees the existence of a mathematical solution outside the horizon, but for a physical solution we must have $r\ge 0$ (meaning ${\partial }_{r}M\left(r,\theta \right)>0$ ) and $M\left(r,\theta \right)\le 0$ for $r\le r_i$ : the mass must become non-positive at the inner radius.

Therefore, for certain {M₀, R₁} values there will be no physical solution, even for the TOV equation.

4. The Solving Process for the Extended Kerr Space-Time

In addition to the fundamental dual parameters {r_i, ρ_i} corresponding to {R₁, M₀} in the rotation-free TOV-case, in the Kerr-case there is the new fundamental parameter Δr_i (inner ellipticity for inner boundary condition), resp. ΔR₁ (outer ellipticity for outer boundary condition), and the angular velocity ω. The outer radii are

$R_{x1}=R_1-\Delta R_1$ and $R_{y1}=R_1$, the latter equality arising from the fact that centrifugal distortion acts only in the x-direction (the y-axis being the rotation axis). The inner radii are correspondingly $r_{xi}=r_i-\Delta r_i$, $r_{yi}=r_i$.

The r-θ-slicing algorithm with an Euler-step obeys the iterative procedure with slice step size h₁ in r, and step size h₂ in θ, starting with the r-boundary at r = R₁ (slice n = 0).

The transition from slice n to n + 1 proceeds as follows.

At slice n all variables and 1-derivatives are known from the previous step, 2-derivatives ¶_rrAi, and ρ are calculated from the 6 equations.

At slice n + 1 the variables and 1-derivatives are calculated by Euler-formula (or Runge-Kutta)

$A{i}_{n+1}=A{i}_n+h_1{\partial }_{r}A{i}_n$

${\partial }_{r}A{i}_{n+1}={\partial }_{r}A{i}_n+h_1{\partial }_{rr}A{i}_n$

The 2-derivatives ¶_rrAi, ¶_rrBi and ρ are again calculated from the 6 significant equations with variables and 1-derivatives inserted from above.

The θ-slicingr-backward algorithm with an Euler-step obeys the iterative procedure with slice step size h₁ in θ as above for r, starting with θ = 0, and solves an ordinary differential equation in r in each θ -step. The boundary condition for the r-odeq is set at r = R₁(θ) (the outer ellipse radius) with Ai = 1, M = M₀My0(θ), ${\partial }_{r}Ai=0$, ${\partial }_{r}M=3{\left(R_1\right)}^2{\rho }_{bc}$, where ρ_bc is the outer boundary value for the density, ρ_bc = 0 for the (non-interacting) neutron-gas in a shell-star and ρ_bc > 0, ρ_bc = ρ_equilibrium for the (interacting) neutron fluid in a neutron star. My0(θ) is the mass-form-factor with the condition ${\displaystyle {\int }_0^{\pi /2}My0\left(\theta \right)\mathrm{Cos}\left(\theta \right)\text{d}\theta }=1$, i.e. the overall mass at the outer boundary is M₀. The inner radius r_i(θ) is reached, when My0(θ) = 0.

The alternative (dual) θ-slicingr-forward algorithm starts with the boundary condition at

$r=r_i(\; \theta \; )$

Ai = 1, M = 0, ${\partial }_{r}Ai=0$, ${\partial }_{r}M=3{\left(r_i\right)}^2{\rho }_{bc}$, where ρ_bc = ρ_i is the inner boundary value for the density, ρ_i is approximately the inner (maximum) density ρ(r_i) from the corresponding TOV-equation, the value must be adapted, so that the resulting total mass is M₀. For the compact neutron star the inner radius r_i(θ) is zero.

In the θ-slicingr-backward algorithm one starts with the outer boundary being the ellipsoid r = R(θ, ΔR₁), where ΔR₁ is the outer ellipticity of the star. In the θ-slicingr-forward algorithm one starts with the inner boundary being the ellipsoid r = r_i(θ, Δr_i), where Δr_i is the inner ellipticity of the star.

At the inner boundary the tangential pressure is uniform, so the density is also uniform and equal to the maximum density, ρ(θ) = ρ_i.

In the actual calculation we were using the θ-slicingr-backward algorithm, because here the boundary condition M = M₀ is achieved automatically, when one starts with My₀(θ) = M₀.

The odeqs in rconsist of the 6 significant Einstein equations eqR00, eqR11, eqR22, eqR33, eqR03, eqR41 for the six variables A0(r, θ), A1(r, θ), A2(r, θ), A3(r, θ), A4(r, θ), M(r, θ) with θ = θ_i and θ-derivatives calculated by Euler-step from the preceding q-slice. For i = 0 i.e. θ = 0 the θ-derivatives are taken from start values for all variables, which normally represent the corresponding TOV-solution (here only A0(r), A1(r), M(r) are non-trivial and do not depend on θ). The odeqs are highly non-linear algebraic differential equations and hard to solve numerically with classical methods for linear odeqs extended by an algebraic equation solver. In the case of a nonlinear odeq-system one uses an Euler or Runge-Kutta method and calculates in each step the highest derivatives with a numerical algebraic equation solver. As an alternative one can use minimization of the least-squares-error in the highest derivatives instead of a numerical algebraic equation solver. Minimization has also the advantage that one can minimize the complete set of Einstein equations plus the 2 additional continuity equations eqR41, eqR42 in the error goal function instead of the 6 significant equations, which improves the stability of the solution (e.g. in case of degeneracy).

The numerical error of the algorithm is calculated from ${\displaystyle {\sum }_i\left\{{\left(e{q}_i\right)}^2,i=1,\cdots ,n\right\}}$ i.e. the Euclidean norm of the equation values (the right side of the Einstein equations being 0). The error is calculated over the lattice {r_i,θ_j} as median, mean or maximum. In the internal loop of the algorithm over r_i at fixed θ_j, the solution of the algebraic discretised Einstein-equation is achieved by square-root error minimization, so it is essential to avoid singularities, e.g. at the horizon and the pseudo-singularity at θ = 0. This is achieved by selecting appropriate analytic convergence factors for the (left side of) the Einstein equations. As the equations are to be zeroed for the solution, the convergence factors do not change the solution of course, but they cancel the numerical singularities, which could otherwise jeopardize the numerical convergence of the algorithm.

The actual calculation was carried out in Mathematica using its symbolic and numerical procedures. In the first stage, the Einstein equations were derived from the ansatz for g_μν from section 2 and simplified automatically. The arising complexity of the equations is such, that it is practically impossible to handle them manually: the Mathematica function LeafCount, which returns the number of terms in the equation, gives the complexity of LeafCount [eqR00] = 17,408, LeafCount [eqR11] = 27,528, LeafCount [eqR22] = 134,929 for the first 3 equations. To verify the equations, the TOV equation was derived by symbolic manipulation for ω = 0a = 0 from eqR00,eqR11,eqR22,eqR41.

The power of Mathematica is sufficient to solve the TOV equation with the single procedure NDSolve. For the full Einstein equations it fails even for the ordinary differential equations (odeq) in rarising for fixed θ. It took us a long time to find an algorithm, which could handle the complexity of the equations and solve them in an acceptable time (3 - 5 minutes on a 16 × 8 lattice) on a PC-desktop and converge in the required region with an acceptable error of around 0.01.

For the second numerical stage we tried several slicing algorithms, and the best alternative proved to be the θ-slicingr-forward algorithm implemented by hand in Mathematica. The solution of the resulting odeq in each r-step was calculated using NDSolve.

Also, for every star model and parameter set, the TOV solution with ω = 0a = 0was calculated first with the algorithm and compared with the exact TOV solution.

5. The TOV Equation as the Limit ω → 0 for the Extended Kerr Space-Time

In the Schwarzschild spacetime ω = 0 and a = 0, we have spherical symmetry, no dependence on θ, then the TOV-equation can be derived from the remaining non-trivial Einstein equations eqR00, eqR11, eqR22, eqR41.

The TOV-equation is in the standard form:

$P\text{'}\left(r\right)=-\left(\frac{GM\left(r\right)\rho \left(r\right)}{r^2}\right)\left(1+\frac{P\left(r\right)}{\rho \left(r\right)c^2}\right)\left(1+\frac{4\pi r^{3}P\left(r\right)}{M\left(r\right)c^2}\right){\left(1-\frac{2GM\left(r\right)}{r{c}^2}\right)}^{-1}$ (13)

and using r_s

$P\text{'}\left(r\right)=-\left(\frac{c^2{r}_s\rho \left(r\right)}{2r^2}\right)\left(1+\frac{P\left(r\right)}{\rho \left(r\right)c^2}\right)\left(1+\frac{4\pi r^{3}P\left(r\right)}{M\left(r\right)c^2}\right){\left(1-\frac{r_{s}M\left(r\right)}{r{M}_t}\right)}^{-1}$,

where

M_t is the total mass, furthermore

$4\pi r^2\rho \left(r\right)=M\left(r\right)$, $P\left(r\right)=k_1\rho {\left(r\right)}^{\gamma }$

In order to make the variables dimensionless, one introduces “sun units”

$r_{ss}=r_s\left(sun\right)=\frac{2G{M}_{sun}}{c^2}=3\text{km},{\rho }_s=\frac{M_{sun}}{4\pi r_{ss}^3/3}=1.76\times {10}^{16}\frac{\text{g}}{{\text{cm}}^{\text{3}}},P_s={\rho }_s{c}^2$

where r_ss Schwarzschild-radius of the sun, ρ_s the corresponding Schwarzschild-density and P_s the corresponding Schwarzschild-pressure.

In “sun units” TOV-equation transforms into

$\begin{array}{l}{P}_1\text{'}\left(r_1\right)r_1^3\left(r_1-M_1\left(r_1\right)M_0\right)\\ =-\frac{1}{2}\left(\frac{M_1\text{'}\left(r_1\right)M_0}{3}+P_1\left(r_1\right)r_1^2\right)\left(M_1\left(r_1\right)M_0+3P_1\left(r_1\right)r_1^3\right)\end{array}$ (14)

with the normalized mass M₁(r₁), and $M_1\left(R_1\right)=1$, or

$P_1\text{'}\left(r_1\right)r_1^3\left(r_1-M\left(r_1\right)\right)=-\frac{1}{2}\left(\frac{M\text{'}\left(r_1\right)}{3}+P_1\left(r_1\right)r_1^3\right)\left(M\left(r_1\right)+3P_1\left(r_1\right)r_1^3\right)$

where $M_0=\frac{M_t}{M_{sun}}$, $M\left(r_1\right)$ is the mass within the radius r, M(r₁) = M₀M₁(r₁) in dimensionless variables r₁, $\rho \left(r_1\right)=\frac{M\text{'}\left(r_1\right)}{3r_1^2}$, M, $P_1=k_1\rho {\left(r_1\right)}^{\gamma }$ and R₁ is the dimensionless radius of the star.

With the replacement $P=k_1{\rho }^{\gamma }$ for the pressure from the equation of state and $\rho =\frac{M_{0}M\text{'}}{3r^2}$ we obtain a diff. equation for M degree 2 in r and we impose the boundary condition inr = R₁:

$M\left(R_1\right)=M_0$, $M\text{'}\left(R_1\right)=0$ for non-interacting Fermi-gas and for an interacting Fermi-gas: $M\left(R_1\right)=M_0$, $M\text{'}\left(R_1\right)=\rho \left(R_1\right)3R_1^2$, $\rho \left(R_1\right)={\rho }_e$, where ρ_e is the equilibrium density in the minimum of V_nn and $P_1\text{'}\left({\rho }_e\right)=0$ (here an equivalent boundary condition is $\rho \text{'}\left(R_1\right)=\infty$ ).

6. The Equation of State and Rotation Parameters

6.1. The Equation of State for an (Non-Interacting) Nucleon Gas

Here, $P=k_1{\rho }^{\gamma }$ is the equation of state of the star, derived from the thermodynamic Fermi gas equation at T = 0 ( [1], chap. 48).

$P=-\frac{\partial E}{\partial V}=8\pi P_0\left(\frac{x_F^3}{3}\sqrt{1+x_F^2}-f\left(x_F\right)\right)$ (15)

$P_0=\frac{m{c}^2}{{\lambda }_c^3}=\frac{m^4{c}^5}{h^3}$, where λ_c is the de-Broglie wavelength of the Fermi gas with particle mass m, ${\lambda }_c=\frac{h}{mc}$.

$x_F=\frac{p_F}{mc}=\frac{{\lambda }_c}{2}{\left(3\pi \right)}^{1/3}{n}^{1/3}$, where x_F is the Fermi-angular-momentum, n the particle density

$f\left(x_F\right)={\displaystyle {\int }_0^{x_F}\text{d}x{x}^2\sqrt{1+x^2}}$

The resulting approximate equations of state for P are

$P=8\pi P_0\left(\begin{array}{l}\frac{x_F^5}{15}\\ \frac{x_F^4}{12}\end{array}\right)=\left(\begin{array}{l}{K}_1{\rho }^{5/3}\rho \ll {\rho }_c\\ K_2{\rho }^{4/3}\rho \gg {\rho }_c\end{array}\right)$ (16)

valid for the density ρ and the critical density ρ_c

${\rho }_c=\frac{m}{{\lambda }_c^3}\frac{8\pi }{3}$

The full expression for P, including temperature T, is as follows ( [5], chap.15).

Here, we use dimensionless variables (r₁ distance unit de-Broglie-wavelength λ_c, V₁ volume unit ${\lambda }_c^3$, n₁ particle density unit $1/{\lambda }_c^3$, E₁ energy unit $E_0=\frac{\hslash c}{{\lambda }_c}=\frac{m{c}^2}{2\pi }$, inverse thermal energy ${\beta }_1=\frac{E_0}{kT}$, chem. potential μ₁ in E₀)., for the gas model we use the Debye model with the state density $D_1\left(E_1\right)=\frac{1}{4{\pi }^{7/2}}\sqrt{E_1}$, maximum energy ${\epsilon }_{F1}=\frac{3^{2/3}{\pi }^{1/3}}{4}{n}_1^{2/3}$, the resulting particle density is

$\begin{array}{c}{n}_1=\frac{N_{op}}{V_1}=\frac{2}{V_1}{\displaystyle {\int }_0^{\infty }\text{d}{\omega }_1}\frac{D_1\left({\omega }_1\right)}{1+\exp \left({\beta }_1\left({\omega }_1-{\mu }_1\right)\right)}\\ =\frac{1}{2{\pi }^{7/2}}{\displaystyle {\int }_0^{\infty }\text{d}{\omega }_1}\frac{\sqrt{{\omega }_1}}{1+\exp \left({\beta }_1\left({\omega }_1-{\mu }_1\right)\right)}\end{array}$

From this relation the chem. potential μ₁ can be calculated, an approximation formula is

${\mu }_1={\epsilon }_{F1}-\frac{{\pi }^2}{12{\beta }_1^2{\epsilon }_{F1}}={\mu }_1(\; n\; 1\; )$

Finally, the resulting pressure (=energy density) $p_1\left({\beta }_1,n_1\right)$ :

$p_1\left({\beta }_1,n_1\right)=\frac{4{\pi }^{3/2}}{3{\pi }^2}{\displaystyle {\int }_0^{\infty }\text{d}{\omega }_1}\frac{{\omega }_1^{3/2}}{1+\exp \left({\beta }_1\left({\omega }_1-{\mu }_1\left(n_1\right)\right)\right)}$ (17)

Below a 3D-diagram of $p_1\left({\beta }_1,n_1\right)$ in dimensionless variables for a nucleon gas (m = m_n, density $\rho =\frac{E_{0}n}{c^2}$ in sun units, E₀ = 149.4 MeV) is depicted in Figure 1 [6].

Here kT is in E₀ units, and one sees the dependence $P=k_1{\rho }^{\gamma }$ except on the left side, when kT reaches the magnitude of 1 Gev (T = 10¹⁰K).

6.2. The Equation of State for an (Interacting) Nucleon Fluid

For the interacting nucleon gas we take into account the nucleon-nucleon-potential in the form of a Saxon-Woods-potential modeled on the experimental data: [7] - [13]

$V_{sw}\left(r,V_0,r_0,d_{r0}\right)=\frac{V_0}{1+\exp \left(\frac{r-r_0}{d_{r0}}\right)}$

$V_{nn}\left(r\right)=-V_{sw}\left(r_a,V_a,r_a,d_{ra}\right)+V_{sw}\left(r_c,V_c,r_c,d_{rc}\right)$

where $V_{nn}$ is thenucleon-nucleon-potential with an attractive part $-V_{sw}\left(r_a,V_a,r_a,d_{ra}\right)$ and a repulsive core $V_{sw}\left(r_c,V_c,r_c,d_{rc}\right)$, the distance r between the nucleons is $r={\left(E_n/\rho \right)}^{1/3}$, where $E_n={\left(m_n{c}^2/\left(2\pi \right)\right)}^{1/3}=149.4\text{MeV}\approx m_{\pi }{c}^2$ is the nuclear energy scale m_π = pion mass = 140 MeV, m_n = neutron mass = 140 MeV.

The Saxon-Wood potential is shown in Figure 2 below.

The pressure of the interacting nucleon fluid becomes then

$P_1\left(r_1\right)=c_1{\rho }_1\left(r_1\right)V_{nn}\left({\left(E_n/{\rho }_1\left(r_1\right)\right)}^{1/3}\right)$ (18)

The experimental data used here are those from [7], and are shown in Figure 3.

And the hard-core potential from the lattice calculation Reid93 [10] is shown in Figure 4.

Both potentials are fitted with a double Saxon-Woods-potential V_nn in Figure 5:

Figure 1. Fermi pressure.

Figure 2. Saxon-Woodpotential V_nn(r) with energy (MeV), r(fm) [6].

Figure 3. Nucleon-nucleon potential experimental data.

Figure 4. Nucleon-nucleon potential lattice calculation.

From the nucleon-nucleon-potential the pressure is calculated taking into account the low-density Fermi-pressure of the nucleons $P_{nf}\left(\rho \right)=K_1{\rho }^{5/3}$ and nucleon-nucleon pressure $P_{nn}\left(\rho \right)=V_{nn}\left({\rho }^{-1/3}\right)\rho$ is shown in Figure 6 [6], and the total pressure $P_{fg}\left(\rho \right)=K_1{\rho }^{5/3}+P_{nn}\left(\rho \right)$ is shown in Figure 7.

This equation-of-state has a minimum at $\rho ={\rho }_c=0.0417$ and $P\text{'}\left(r\right)=1$ at $\rho ={\rho }_m=0.0544$.

Figure 5. Nucleon-nucleon potential data fitted with Saxon-Wood-formula, with r(fm), V(MeV).

Figure 6. Nucleon-nucleon Fermi-pressure.

Figure 7. Total pressure P_fg(r), pressure P and density rshown in sun-units.

As the sound velocity $v=\frac{\text{d}P\left(\rho \right)}{\text{d}\rho }$, v > 0 and v < 1 (i.e. subluminal), the admissible density range in the neutron-fluid model is ${\rho }_c\le \rho \le {\rho }_m$.

6.3. Maximum Omega-Values in Kerr-Space-Time

We consider here a rotation model with constant angular velocity ω. With this model the resulting 4-velocity u^μ has the form [4] [14] [15]:

$u^{\mu }=\left(u^0,0,0,\omega u^0\right)$

The maximum values for ω are calculated from the minimal zeros in omega of the denominator in u⁰ from (9a), minimized over r₁ and th in their respective regions

$r_i\le r_1\le R_1$ and $0\le \theta \le \pi /2$.

The resulting value is $\omega \le \frac{1}{2R_1\sqrt{{\alpha }_f}}$, where α_f is the form-factor in the moment of inertia I₁.

$I_1={\alpha }_{f}M{R}_1^2$, α_f = 2/3 for a shell, α_f = 2/5 for a sphere.

With non-vanishing density the actual ω_max depends on {Ai, ρ}, and has to be calculated from the above expression for u⁰.

A less stringent limit for omega can be deduced from the limit for the parameter a in Kerr-spacetime (in sun-units and c = 1): $a\le \frac{r_s\left(M_0\right)}{2}=\frac{M_0}{2}$, and $a=\frac{{\alpha }_f{M}_0{R}_1^2\omega }{M_0}$ therefore $\omega \le \frac{M_0}{2R_1^2{\alpha }_f}$

7. The TOV-Equation: A New Ansatz

Generally speaking, the parameters of the solution are (in brackets denomination in program code): angular momentum radius a (=alpha1, =0 for TOV), the factor in the state equation k₁, the power in the state equationγ (=gam), radius R, mass M₀, the relative radius uncertainty $dr{02}_{rel}$ (=dr02rel), the moment of inertia factor $f_I$ (= infac, 0 for TOV), the singularity smoothing parameter ε (=epsi, see below), and the boundary factor n_rmax (=nrmax). Here the boundary factor enters the upper boundary of the TOV differential equation as $r_{max}=R\left(1+n_{rmax}dr{02}_{rel}\right)$.

The dimensionless TOV-equation is a differential equation in the mass M(r) of degree 2, and is highly non-linear, the dimensionless mass-density relation is $\rho =\frac{M\text{'}}{3r^2}$.

The customary way of solving the TOV equation is to impose the boundary condition at r = 0 with M(0) = 0, $M\text{'}\left(0\right)=3r^2{\rho }_0$ where ρ₀ the maximum central density.

In the new ansatz for the mass M(r) we impose the outer boundary condition at r = R₁:

for a pure Fermi-gas without interaction: $M\left(R_{\text{1}}\right)=M_0$, $M\text{'}\left(R_1\right)=\rho \left(R_1\right)3R_1^2$, $\rho \left(R_1\right)=0$ ;

for an interacting Fermi-gas: $M\left(R_1\right)=M_0$, $M\left(R_1\right)=\rho \left(R_1\right)3R_1^2$, $\rho \left(R_1\right)={\rho }_e$, where ρ_e is the equilibrium density in the minimum of V_nn and $P_1\text{'}\left({\rho }_e\right)=0$ (here an equivalent boundary condition is $\rho \text{'}\left(R_1\right)=\infty$ ).

The star parameters mass M₀ and radius R₁, which enter the outer boundary condition determine completely the solution. In general, there will be an inner radius r_i > 0 with the maximum density ${\rho }_0=3r^{2}M\text{'}\left(r_i\right)$ and M(r_i) = 0. The corresponding “dual” parameters are the inner radius r_i and the maximum density ρ₀. One can show that for ${\rho }_0\gg {\rho }_c$ (where ρ_c is the critical density of the equation of state) there is no solution with a compact star r_i = 0, i.e. there is a maximum mass M_c for the TOV equation, in case of compact neutron stars M_c = 3.04M_sun (see below). As we will see, there is in general a solution, if we allow r_i > 0 and impose an outer boundary condition at r = R₁, as long as R₁ is not too close to the Schwarzschild radius r_s = M₀ of the star. In the limit $R_1\to r_s$ there will be no positive zero of M(r), i.e. r_i < 0 and the resulting (mathematical) TOV-solution will be no physical solution. But in general, speaking naively, the gravitational collapse of the star is avoided for large masses (M₀ > M_c), if it has a shell structure with the inner radius r_i and the outer radius R₁ > M₀.

As we will see, this outer boundary condition together with allowing r_i > 0 changes dramatically the resulting manifold of physical solutions.

7.1. The TOV-Equation: The Parametric Solution and Resulting Star Types

By setting-up a parametric solution of the TOV-equation one gets a map of possible physical solutions, i.e. possible star structures. As parameters one can use either (M₀, R₁) in the outer boundary condition at r₁ = R₁ or the dual parameter pair (r_i, ρ_bc) in the inner boundary condition r₁ = r_i.

The pure neutron Fermi-gas model yields for compact neutron stars a maximum mass of M_maxc = 0.93M_sun, which is in disagreement with observations. Therefore, at least for compact neutron stars, a model of interacting neutron fluid must be used. In 6.2 above we have described a Saxon-Wood-potential model for the nucleon-nucleon interaction, which seems to fit the experiment and the theory in the best way. There will be a critical density (dependent on temperature of course), where a transition from interacting fluid to Fermi-gas takes place, it is plausible to set this density equal to the Saxon-Wood critical density

${\rho }_c=0.0417$.

We made calculations with the TOV-equation using these two models for neutron-based stars and we came to the conclusion that compact neutron stars with mass M₀ ≤ 3.04M_sun consist of interacting neutron fluid and neutron shell-stars for M₀ ≥ 5M_sun obey the Fermi-gas model. The underlying calculation is the Mathematica-notebook [6].

This approach yield results, which are described below.

Neutron stars consist of interacting neutron fluid and are compact stars with (M₀, R₁) = (0.14, 1.49), …, (3.04, 3.95) and the maximum density 0.048 ≤ ρ_bc ≤ 0.0544 = ρ_bcmax in sun-units, or shell-stars with ρ_bc ≥ ρ_bcmax and (M₀, R₁) = (3.04, 3.95), …, (4.91, 4.92), neutron star R-M-relation follows approximately a cubic-root-law: R~M^1/3.

Stellar shell-stars consist of (almost) non-interacting Fermi-gasof neutrons and are thin shell-stars with R₁ > r_s, R₁ ≈ r_s, r_i ≈ r_s, i.e. the shell is close to the Schwarzschild-radius and its outer edge outside the Schwarzschild-horizon with max. density 0.0025 ≤ ρ_bc ≤ 0.042, and obey an almost linear R-M-relation (M₀, R₁) = (5.5, 9.1), …, (81.3, 91.2), independent of ρ_bc for ρ_bc ≥ 0.028, with redshift factor around 50 for M = 80M_sun.

Galactic (supermassive) shell-stars are very thin shell-stars, which obey the equation-of-state of a white-dwarf (i.e. gravitation counterbalanced by Fermi-pressure of electron gas) and have an almost linear R-M-relation with redshift factor 20…100.

Neutron stars

The parametric solution of the TOV-equation has been carried out for the parameters (ρ_bc, r_i) at the boundary r = r_i, in the range: density 0.02 ≤ ρ_bc ≤ 0.15 and inner radius 0.01 ≤ r_i ≤ 15, yielding physical solutions for density 0.048 ≤ ρ_bc ≤ 0.0544 = ρ_bcmax andinner radius 0.01 ≤ r_i ≤ 3. The TOV-equation is solved for M(r) and ρ(r), and a physical solution is a mathematical solution with M ≥ 0 and ρ ≥ 0, ρ' ≤ 0 and subluminal equation-of-state within a certain interval r = {r_i,r₀₂}, which reaches a point, where M'(r) = 0 and ρ(r) = 0. The radius R₁ and the total mass M₀ is reached at M'(R₁) = 0, the physical solution ends there.

The validity interval for ρ is explained by the fact, that the sound velocity $v_s\left(\rho \right)=\frac{\partial P\left(\rho \right)}{\partial \rho }$ must be positive and below 1 (subluminal in c-units).

The parametric mapping of the solutions results in the following dependence for M₀(r_i, ρ_bc), R₁(r_i, ρ_bc) (r_i, ρ_bc, M₀,R₁ insun-units) (Figure 8 [6]):

For r_i = 0 the mapping describes the compact neutron stars, resulting in R₁(M₀)function (Figure 9(a), Figure 9(b), Figure 10 [6]):

The R-M-relation follows approximately a cubic-root-law: R~M^1/3, with a range of (M₀, R₁) = (0.14, 1.49), …, (3.04, 3.95), i.e. the resulting maximum compact mass is M_maxc = 3.04 M_sun.

For M₀ ≥ M_maxc the function R₁(r_i = const,ρ_bc) is flat or slightly decreasing with ρ_bc, so one expects the stable configuration to be the one with maximum ρ_bc = ρ_bcmax, with a range of (M₀, R₁) = (3.04, 3.95), …, (4.91, 4.92).

The admissible mass range ends, where the thickness of the shell above the Schwarzschild-radius becomes very small (minimum 0.01).

So in total the R-M-relation for neutron stars becomes

The maximum mass for a repulsive-hardcore-model for the equation-of-state DD2 [16] is 2.42M_sun, from our mapping we have the maximum compact neutron star mass of M_maxc = 3.04M_sun.

Figure 8. Neutron star: mass M and radiusR in dependence of inner boundary density ρ_bc and inner radius r_i.

(a)(b)

Figure 9. (a) Compact neutron star: mass-radius R-M relation; (b) Shell neutron star: mass-radius R-M relation.

Figure 10. Neutron star: mass-radius R-M relation.

The actual theoretical limit for neutron star core density is ρ_max = 3.5 × 10¹⁵ g/cm³ = 0.199 in sun-units [8] [9].

The limit for ρ_bc reached in our mapping is only ¼ of this ρ_bc = ρ_bcmax = 0.0544, due to the subluminal-sound-condition and the use of an (attractive) nucleon-nucleon-potential for the nucleon-fluid instead of a pure repulsive-hardcore-model.

The classical argument for the collapse of a neutron star to a black-hole for ρ_bc > ρ_max, dating back to Oppenheimer [1], is invalidated here by the simple introduction of shell-star models, where r_i > 0, and therefore there is no mass at the center, which means physically, there is only a very diluted nucleon gas there.

Stellar shell-stars (stellar black-holes)

We assume that the underlying equation-of-state state for stellar shell-stars is the Fermi-gas of nucleons with the low-density limit of

$P\left(\rho \right)=K_1{\rho }^{5/3}$.

We make a further plausible assumption that the “edge” of the solution mapping are the physically stable solutions, i.e. the R-M-relation for stellar shell-stars. The edge in this case consists for fixed ρ_bc < 0.0417 = ρ_oc of solutions with maximum r_i (because then the average density in the shell is lowest) and for ρ_bc = 0.0417 = ρ_oc it consists of the solutions (M₀, r_i, R₁) at the right boundary (7 ≤ r_i ≤ 22,M₀ ≥ 4.91 = M_maxn), where M_maxn is the maximum mass for the neutron star interactive neutron fluid. We assume that at ρ_bc = ρ_octhe state transition from the interactive neutron fluid to neutron Fermi gas takes place. This is almost certainly an oversimplification, still we believe that the model describes the reality correctly in principle.

The chosen parameter range of the solution mapping is: density 0.025 ≤ ρ_bc ≤ 0.0417 = ρ_oc and inner radius 0.01 ≤ r_i ≤ 90, where ρ_oc is equilibrium value of the nucleon-nucleon-potentials with Pnf'(ρ_oc) = 0, the transition point from the nucleon-fluid to the nucleon-gas phase.

The “edge” of the mapping yields the (M₀-,r_i, R₁-)-range of (M₀, r_i, R₁) = (5.35, 7, 8.49), …, (81.3, 89.6, 91.2), where the upper limit is in fact mathematically open, but the “thinning-out” of the solutions for small ρ_bc and large r_i makes it physically plausible (see Figure 11 [6] below).

The resulting R-M-relation is practically linear and has a maximum mass value of M_max = 81.3M_sun. (Figure 12(a)).

And the corresponding relative shell thickness dR_rel = dR/M is (Figure 12(b)) and the relative Schwarzschild-distance dR_srel = (R−M)/M is (Figure 13)

The inverse of dR_srelgives roughly the light attenuation factor of {1.7, …, 20.}. Taken the attenuation factor and the small relative shell thickness of around 0.02, these stellar shell-stars have approximately the properties expected of a genuine black-hole, when measured from a distance $r\gg R_1$.

Figure 11. Stellar shell star: mass M and radius R in dependence of inner boundary density ρ_bc and inner radius r_i.

(a)(b)

Figure 12. (a) Stellar shell star: mass-radius R-M relation; (b) Stellar shell star: relative shell thickness.

Figure 13. Stellar shell star: relative Schwarzschild-distance.

Entropy of a thin shell star

The celebrated Bekenstein-Hawking formula for the entropy of a black hole reads [17]:

$S=\frac{k_{B}A}{4L_P^2}$, where A is the surface area, k_B the Boltzmann-constant, and L_P the Planck-length.

The entropy of a (cold) shell-star with radius R and thickness dR in the limit R = r_s, $dR\ll r_s$, with all particles in the lowest possible energy state, can be easily calculated from the Boltzmann-formula

$S=k_B\ln W$, where W is the number of possible micro-states.

With the elementary area $\pi L_{min}^2$, where $L_{min}=L_P$, W becomes $W=2^{A/\left(\pi L_{min}^2\right)}$ (each of the $N=A/\left(\pi L_{min}^2\right)$ area elements can be occupied or empty), so $S=\frac{k_{B}A}{\pi L_P^2}\ln 2$, which is identical to the Bekenstein-Hawking entropy with the factor (ln2)4/π = 0.882.

Galactic (supermassive) shell-stars

The mean density of a black-hole scales with its radius R like $\rho \left(R\right)=\frac{M}{V}=\frac{R}{\left(4/3\right)\pi R^3}=\frac{3}{4\pi R^2}$

i.e. for supermassive black-hole with M = 10⁶M_sun we have $\rho \approx {10}^{-12}$ in sun units (su).

In the following we use the abbreviation MM_sun = 10⁶M_sun.

The density scale of a white-dwarf star is 10⁶ g/cm³ = 5.7 × 10⁻¹¹ su [1]. Therefore it is plausible to try a parametric mapping with the white-dwarf equation-of-state, where the underlying Fermi-pressure is that of an electron gas instead of a nucleon gas, i.e. equation-of-state $P_1\left(r_1\right)=k_1{\rho }_1{\left(r_1\right)}^{\gamma }$ for a pure Fermi gas,γ = 5/3 if the density is below the critical density ρ_c.

The results for M₀(r_i, ρ_bc), R₁(r_i, ρ_bc) are shown below (Figure 14 [6]):

Figure 14. Galactic shell star: mass M and radius R in dependence of inner boundary density ρ_bc and inner radius r_i.

From this result one can draw several consequences: first, the actual density is around 10⁻¹², that is well below the critical density for a white-dwarf of ρ_c = 0.91 × 10⁶ g/cm³ = 5.17 × 10⁻¹¹ su: γ = 5/3 in the equation-of-state is justified. Second, the viable solutions lie to the left of a “ridge” reaching up to masses around 30 MM_sun. Third, a stable solution for a fixed mass will have the highest possible maximum density ρ_bc and that will lie on the “ridge”. So one can calculate the R-M-relation following the “ridge”.

The resulting R-M-relation is as follows (Figure 15).

And the inner radius is (Figure 16).

The R-M-relation is almost linear, as expected, and goes up to 50MM_sun. $d{R}_{rel}=\left(R_{\text{1}}-r_i\right)/M_0$ is the relative thickness (Figure 17), and shows, that the shells are very thin indeed, with a minimum of 0.001. The fourth diagram shows the relative Schwarzschild-distance (Figure 18) $d{R}_{srel}=\left(R_{\text{1}}-M_0\right)/M_0$, which has a minimum at {M₀, dR_srel} = {7, 0.00142857}, so that its reciprocal value (approximate light attenuation factor) is around 700. So the overall result is, that the supermassive shell-stars become ever thinner shells, while the distance from the Schwarzschild-horizon is increasing.

7.2. The TOV-Equation: A Case Study for Typical Star Types

In the nearly-rotation-free case the solution of the TOV-equation was calculated for 4 models in sun units with r_s = Schwarzschild radius

Figure 15. Galactic shell star: R-M-relation.

Figure 16. Galactic shell star: inner radius.

Figure 17. Galactic shell star: relative shell thickness.

Figure 18. Galactic shell star: relative Schwarzschild distance.

$r_{ss}=r_s\left(sun\right)=\frac{2G{M}_{sun}}{c^2}$, ${\rho }_s=\frac{M_{sun}}{4\pi r^3/3}$, $P_s={\rho }_s{c}^2$,

$r_{ss}=\text{3}\text{km}$, ${\rho }_s=1.76\times {10}^{16}\text{g}/{\text{cm}}^{\text{3}}$, $M_{sun}=3\times {10}^{30}\text{kg}$,

• average compact neutron star with mass M₀ = 0.932M_sun, radius R₁ = 2.767r_ss.

• maximum mass neutron shell-star M₀ = 4.91, R = 4.926.

• white dwarf with M₀ = 0.6M_sun, radius R₁ = 3000r_ss.

• stellar black hole with M₀ = 15.69M_sun, radius R₁ = 17.89r_ss, inner radius r_i = 17r_ss.

• galactic black hole with M₀ = 4.367 × 10⁶M_sun, R₁ = 4.380 × 10⁶r_ss, r_i = 4.356 × 10⁶r_ss.

Compact neutron star

parameters = {k1 = 0.40, gam = 5/3, M0 = 0.932, R1 = 2.76, rhobc = 0.0456, ri = 0.01};

The mean density is here ${\rho }_{mean}=\frac{M_0}{R_1^3}=0.04447$.

The critical density of the neutron Fermi gas with neutron mass m_n is ${\rho }_{cn}=\frac{m_n^4{c}^3}{3{\pi }^2{\hslash }^3}=0.35$ (see [17]), so the low-density approximation with γ = 5/3 can be used.

Results TOV:

density rho (Figure 19), mass M (Figure 20) are

As can be seen in the ρ-diagram, the derivative $\rho \text{'}\left(R_1\right)=\infty$, because there the equilibrium density ρ_c with $P\text{'}\left(r_c\right)=0$ in the pressure is reached.

Maximum mass neutron shell-star

parameters = {k1 = 0.40, gam = 5/3, M0 = 4.91, R1 = 4.926, rhobc = 0.0544, ri = 3};

The mean density is here ρ_mean = 0.0530.

Results TOV:

density rho (Figure 21), mass M (Figure 22) is

Stellar shell-star (stellar black hole)

parameters = {k1 = 0.40, gam = 5/3, M0 = 15.69, R1 = 17.89, rhobc = 0.0359, ri = 17};

The mean density is here ρ_mean = 0.0194.

The resulting rho (Figure 23) and M (Figure 24) are:

Figure 19. Typical compact neutron star: density in dependence of radius.

Figure 20. Typical compact neutron star: mass in dependence of radius.

Figure 21. Shell neutron star: density in dependence of radius.

Figure 22. Shell neutron star: mass in dependence of radius.

Figure 23. Typical stellar shell star: density in dependence of radius.

Figure 24. Typical stellar shell star: mass in dependence of radius.

Here the radius R₁ is reached, when M'(r₁ = R₁) = 0, i.e. ρ(R₁) = 0.

White-dwarf star

parameters = {k1 = 1.43 × 10⁶, gam = 5/3, M0 = 0.6, R1 = 3000, rhobc = 2.02 × 10⁻¹¹, ri = 0};

The underlying state equation is that of a small-momentum electron Fermi-gas with the critical density [1] ${\rho }_{cw}=\frac{m_e{m}_n^3{c}^3}{3{\pi }^2{\hslash }^3}=0.517\times {10}^{-10}\text{su}$.

The mean density is here ρ_mean = 2.22 × 10⁻¹¹, the maximum deviation of ρ is Δ_maxρ = 0.21 × 10⁻¹¹, so the density is practically constant, as expected.

The solution of the TOV-equation yields density rho (Figure 25), mass M (Figure 26).

Galactic shell-star (supermassive black-hole)

parameters = {k11(γ = 5/3) = 0.0243 × 10⁶, k12(γ = 4/3) = 0.067 × 10⁴, M0 = 4.367 × 10⁶, R1 = 4.380 × 10⁶, rhobc = 4.934 × 10⁻¹², ri = 4.356 × 10⁶};

TOV equation was solved with an exterior boundary condition r₀₂ = R₁(M(r₀₂) = M₀,M'(r₀₂) = 0), which is equivalent to the interior boundary condition r₀₁ = r_i(M(r₀₁) = 0,ρ(r₀₁) = ρ_bc), and with the full Fermi-gas equation-of-state instead of the simple power law $P\left(\rho \right)=K_1{\rho }^{\gamma }$.

The mean density is here ρ_mean = 3.16 × 10⁻¹².

The “naive” mean density is here ${\rho }_{mean}=\frac{M_0}{R_1^3}=3.16\times {10}^{-12}$, i.e. by a factor 10 lower than the mean density of white dwarf. Therefore, despite its huge mass, the galactic black hole can be described by the state equation of a small-momentum (undercritical) Fermi electron gas with the relative density $x_F=\frac{\rho }{{\rho }_{cn}}=0.0612$ much smaller than that for the white dwarf.

Figure 25. Typical white dwarf star: density in dependence of radius.

Figure 26. Typical white dwarf star: mass in dependence of radius.

TOV-solution for rho (in 10⁻¹² units, Figure 27), M (in 10⁶ units, Figure 28) in r (in 10⁶ units), is:

Here there is an internal “hole” with a radius r_i = 4.356 × 10⁶, maximum ρ = 4.934 × 10⁻¹² at r_i. The inner radius r_i lies a little below the Schwarzschild-radius r_s = M₀. The relative shell thickness

Figure 27. Typical galactic shell star: density in dependence of radius.

Figure 28. Typical galactic shell star: mass in dependence of radius.

$d{R}_{rel}=\left(R_1-r_i\right)/M_0=0.00551$, the relative Schwarzschild-distance $d{R}_{srel}=\left(R_1-M_0\right)/M_0=0.00290$, the light attenuation factor is roughly 1/dR_srel = 344.

Furthermore, r_i is little sensitive to the temperature up to T = 10⁷ K.

As for a stellar black hole, when R converges to r_s = M₀, so does the inner radius r_i, and there is no physical solution (with positive ρ and M) for a boundary within the horizon.

8. The Three Star Models for Kerr-Space-Time with Mass and Rotation

The calculation of Kerr-space-time with mass and rotation was carried out for 3 star models:

- a typical compact neutron star with mass around 1 solar mass;

- a presumably typical stellar shell-star with a mass around 15 solar masses;

- a comparatively small galactic shell-star modelled on the central black-hole in the Milky Way with a mass of around 4 million solar masses.

The angular velocity ω was chosen at 0.36ω_max, i.e. about 1/3 of the maximum value.

We are using the so called “sun units” sun Schwarzschild-radius $r_{ss}=r_s\left(sun\right)=\text{3}\text{km}$, sun mass $M_s=M\left(sun\right)=3\times {10}^{30}\text{kg}$, sun Schwarzschild-density ${\rho }_s=\frac{M_s}{4\pi r_{ss}^3}=1.76\times {10}^{16}\text{g}/{\text{cm}}^{\text{3}}$, sun Schwarzschild-pressure $P_s={\rho }_s{c}^2=1.58\times {10}^{35}\text{J}/{\text{m}}^{\text{3}}$ for radius r, mass M, density ρ, and pressure P, respectively.

The mass element here is M₁(r, θ)drdθ and the ring mass M₁(θ) is the differential mass of the θ-beam $\text{d}\theta M_1\left(\theta \right)=\text{d}\theta M_1\left(R_1\left(\theta \right),\theta \right)$, the density ρ is $4\pi \mathrm{Cos}\left(\theta \right)\rho \left(r,\theta \right)\text{d}r\text{d}\theta ={\partial }_r{\partial }_{\theta }{M}_1\left(r,\theta \right)$.

As for the result values, dthrel is the maximum relative angular deviation (in θ), and error (relative to the test function error): wavefront error is the median (on lattice) algorithm error, the interpolation and the Fourier fit error is the error of the respective fit of the discrete solution on the lattice.

We are using here the θ-slicing r-backward algorithm for the shell-stars. For each of the star models a verification step is run first with the angular velocity ω = 0, the result must be the same as in the corresponding TOV-equation. Then a parameter case-study is made for different outer boundary ellipticities ΔR₁ at the outer boundary condition in order to find ΔR₁with a minimal mean energy density: this is the stable solution of the Kerr-Einstein equations.

The algorithm yields as the result a pointwise array of values and first and second derivatives of the variables. The variables for the 6 Einstein equations are the 5 Kerr correction-factor functions A0(r, θ), …,A4(r, θ), and ρ(r, θ).

We are using the θ-slicing r-forward algorithm for the compact neutron star. The variables for the 6 Einstein equations are the 5 Kerr correction-factor functions A0(r, θ), …, A4(r, θ), and M₁(r, θ).

The value denomination is:

a = alpha1 with $a=\frac{J}{Mc}$ the angular momentum radius (amr) of the Kerr model,

ω = omega1 is the angular velocity, R₁ = R1 = r02 is the outer radius, M₀ = M0 is the total mass,

r₁ radius variable, th angle variable,

M(r, θ) = M1(r1, th) is the mass function, A0(r, θ), …, A4(r, θ) Kerr correction-factor functions,

ρ(r, θ) = rho(r1, th) is the density function,

k₁ is the parameter in the approximate Fermi-gas equation-of-state $P_1\left(r_1\right)=k_1{\rho }_1{\left(r_1\right)}^{\gamma }$,

γ = gam, gam1, gam2 is the exponent,

infac is the moment of inertia factor α_f,

epsi is the singularity cancellation parameter with limit(epsi) = 0 introduced to improve the numerical stability in singularities

r_i = riact is the polar inner radius R_y

Δr_i the inner ellipticity is the difference between the polar r_iy and the equatorial inner radius r_ix, Δr_i = r_iy − r_ix

ΔR₁ Is the outer ellipticity, with outer radii R_x1 = R₁− ΔR₁ and R_y1 = R₁

rilow is the minimal radius r₁ reached in the solution

ρ_bc = rhobcx is the boundary condition density

dthrel is the maximum relative difference of a value dependent on θ, e.g.

dthrel(R₁) = (max(R₁(θ)) − min(R₁(θ)))/mean(R₁(θ))

shell thickness dR₁ = R₁(θ)−Δr_i(θ)

Typical rotating compact neutron star

The underlying star model here is a compact (r_i = 0) neutron star of neutron liquid (i.e. strongly interacting neutrons), mass M₀ = 0.932 sun-masses, radius R₁ = 2.76 sun-Schwarzschild-radii r_ss (r_ss = 3.0 km), ω = 0.108688. The underlying calculation is the Mathematica-notebook [18], the results in [19].

The full parameters are:

{alpha1 = 0.331224, omega1 = 0.108688, k1 = 0.4, gam = 5/3., gam1 = 5/3., gam2 = 4/3., M0 = 0.932, dr02rel = 0.33, infac = 2/5., epsi = 0.1, rilow = 0.001, rhobcx = 0.0456}

The r-forward solution is first calculated with the lattice {n_x = 33, n_y = 17} for the rotation-free TOV-case with a TOV-solution as the initial function.

The solution for the Kerr-case starts with this corrected TOV-solution and yields the values:

outer radius $R_1\left(\theta \right)=\left\{2.83912,\cdots ,2.83722\right\}$, mean = 2.83897, dthrel = 0.00118.

total mass M02eff = 0.932.

error: med(err) = 0.0639 wavefront, = 0.0491 interpolation, = 0.0492 Fourier fit.

mean energy density = 0.0791.

The resulting density is depicted in Figure 29, mass M in Figure 30, and the effective radius over azimuthal angle th in Figure 31.

The density distribution is similar to the TOV-case but with a decrease in θ-direction.

Figure 29. Typical rotating compact neutron star: density over x = radius r₁, y = angle th.

Figure 30. Typical rotating compact neutron star: (ring) mass profile for th = 0.1 (equatorial) and th = 1.4708 (polar), the outer edge is reached when M₁(r₁) = M₀.

Figure 31. Typical rotating compact neutron star: effective radius r_02e over angle th.

The rotation results show very small flattening in the polar direction ofdthrel = 0.00118. The neutron star behaves like a fluid because of its “viscosity”, that is, its nuclear interaction and becomes “pumpkin-like”.

Typical rotating stellar shell-star

The star model here is a shell-star (r_i > Schwarzschild-radius) with mass M₀ = 15.69 sun-masses, radius R₁ = 17.89 sun-Schwarzschild-radii r_ss(r_ss = 3.0 km), and angular frequency ω = 0.0126.

The outer Kerr-horizon is r₊ = 15.21. The underlying calculation is the Mathematica-notebook [20], the results in [19].

The outer ellipticity ΔR₁ is at first a free parameter and calculated from a case-study of minimal mean energy density to ΔR₁ = 0.3 = 0.0168R₁.

The full parameters are

{apha1 = 2.68844, omega1 = 0.0126, k1 = 0.4, R1 = 17.89, gam = 5-/3., gam1 = 5/-3., gam2 = 4-/3., M0 = 15.69, infac = 2/-3., epsi = 0.1, rilow = 15.9, rhobcx = 0.036}

The r-backward solution is first calculated with the lattice {n_x = 17, n_y = 9} for the rotation-free TOV-case as the initial function.

Then a case study with the parameter ellipticity ΔR₁ is carried out in order to find the minimal mean energy density, on the set of values ΔR₁ = {0, −0.3, 0.3, 1., 1.6}

The case study yields a minimum at ΔR₁ = 0.3 (cigar-like outer boundary), with amean energy density = 0.017004.

This energy-minimalsolutionwith ΔR₁ = 0.3 yields the values:

outer radius $R_1\left(\theta \right)=\left\{17.59,\cdots ,17.89\right\}$, mean = 17.739, dthrel = 0.0164.

inner radius $r_i\left(\theta \right)=\left\{16.704,\cdots ,16.982\right\}$, mean = 16.823, dthrel = 0.0162.

total mass M02eff = 15.69, inner boundary max(ρ_bc) = 0.035955.

shell thickness dR₁: mean = 0.916, dthrel = 0.0448.

mean energy density = 0.0170.

error: med(err) = 0.00238 wavefront, = 0.00573 interpolation, = 0.00786 Fourier fit.

It is interesting to make a comparison with the spherical-outer-boundary solution.

with ΔR₁ = 0:

inner radius $r_i\left(\theta \right)=\left\{17.0033,\cdots ,17.0243\right\}$, mean = 17.020, dthrel = 0.00123.

total mass M02eff = 15.69, inner boundary max(ρ_bc) = 0.0375.

shell thickness dR₁: mean = 0.870, dthrel = 0.0241.

mean energy density = 0.017698.

error: med(err) = 0.00224 wavefront, = 0.00545 interpolation, = 0.0098 Fourier fit.

The solution with the next higher ellipticity ΔR₁ = 1.0 has the values:

outer radius $R_1\left(\theta \right)=\left\{16.89,\cdots ,17.88\right\}$, mean = 17.38, dthrel = 0.055.

inner radius $r_i\left(\theta \right)=\left\{16.01,\cdots ,16.955\right\}$, mean = 16.438, dthrel = 0.0556.

total mass M02eff = 15.69, inner boundary max(ρ_bc) = 0.0359.

shell thickness dR₁: mean = 0.940, dthrel = 0.0894.

mean energy density = 0.017249.

error: med(err) = 0.00309 wavefront, = 0.00516 interpolation, = 0.00792 Fourier fit.

The two significant non-spherical features are the relative shell thickness variation dthrel(dR₁) and the relative inner ellipticity dthrel(r_i). The first depends roughly linearly on the outer ellipticity ΔR₁, plus the value at ΔR₁ = 0 (dthrel(dR₁) = 0.0241), which is results from rotation. The second, dthrel(r_i), is almost equal to the relative outer ellipticity dthrel(r_i), plus the small amount at ΔR₁ = 0 (dthrel(R₁) = 0.00123).

The density distribution is shown in Figures 32-34.

The density distribution increases in θ-direction (θ = 0 is equatorial, θ = π/2 axial).

The mass distribution is shown in Figure 35, Figure 36.

The physical mass distribution ends at the inner boundary at r_i = 16.7, where the density jumps to ρ = 0.

Figure 32. Typical rotating stellar shell star: density over x = radius r₁, y = angle th.

Figure 33. Typical rotating stellar shell star: density over angle th at the inner boundary r₁ = 16.704.

Figure 34. Typical rotating stellar shell star: the fitted density profile for th = 0.1 (equatorial).

Figure 35. Typical rotating stellar shell star: effective inner radius r_i over angle th.

A remarkable result, distinct from the case of the neutron star, is the shape with rotation. The energy-minimal stellar shell-star behaves like a ball of neutron gas (negligible interaction) and decreases slightly its equatorial radius, so that, speaking naively, the increased gravitation counteracts the centrifugal force, the shell-star becomes “cigar-like”, with the shell thickness approximately constant.

Figure 36. Typical rotating stellar shell star: Fourier-fitted total mass in dependence of radius.

The stellar shell-star has all its mass concentrated within a thin shell

(dR₁ = 0.916) which is situated outside its Schwarzschild-radius M₀ = 15.69, where the minimum distance from the horizon is

min(R₁(θ)) − M₀ = 1.01,

therefore the light energy loss is approximately (in Newtonian approximation)

M₀/min(R₁(θ)) = 0.9395 and the attenuation factor 1/(1 − M₀/min(R₁(θ))) = 16.53, it means that visible green light of 0.514 μm is shifted to 8.50 μm into middle-range infrared.

Typical rotating galactic shell-star

This is modelled (approximately) on the central black-hole in the Milky Way with mass M₀ = 4.36 mega-sun-masses (MM_s), radius R₁ = 4.38 mega-sun-Schwarzschild-radii (13.14 × 10⁶ km, Mr_ss) [21].

The underlying calculation is the Mathematica-notebook [22], the results in [23].

The outer Kerr-horizon is r₊ = 4.26Mr_ss.

In order to maintain numerical performance, we are using for mass and distance 10⁶ (mega) units 10⁶ M_s and 10⁶ r_ss and for density 10⁻¹² (mega⁻²) unit 10⁻¹² ρ_s.

Like in the case of the stellar shell-star, the outer ellipticity ΔR₁ is at first a free parameter and calculated from a case-study of minimal mean energy density to ΔR₁ = −2 dTOV, where dTOV is the shell thickness of the spherical shell-star dTOV = 0.057.

The full parameters are:

{alpha1 = 0.670047, omega1 = 0.05239, k1 = 0.0243, k2 = 0.067, R1 = 4.38, gam =, 5-/3,

gam1 =, 5/-3, gam2 =, 4-/3, M0 = 4.36731, infac =, 2/-3., epsi = 0.0024, rilow = 4.3, rhobcx = 4.915}

The r-backward solution is first calculated with the lattice {n_x = 17, n_y = 9} for the rotation-free TOV-case as the initial function.

Then a case study with the parameter ellipticity ΔR₁ is carried out in order to find the minimal mean energy density, on the set of values ΔR₁ = {0.1, −0.3, −1, -2}*dTOV.

The case study yields a minimum at ΔR₁ = −2 dTOV = −0.114 (pancake-like outer boundary), with a mean energy density = 1.734.

This energy-minimalsolution yields the values:

outer radius $R_1\left(\theta \right)=\left\{4.494,\cdots ,4.38\right\}$, mean = 4.43654, dthrel = 0.025367.

inner radius $r_i\left(\theta \right)=\left\{4.46456,\cdots ,4.34109\right\}$, mean = 4.40029, dthrel = 0.02765.

total mass M02eff = 4.3673, inner boundary max(ρ_bc) = 4.96075.

shell thickness dR₁: mean = 0.035256, dthrel = 0.2507.

mean energy density = 1.73479.

error: med(err) = 0.0122 wavefront, = 0.000404 interpolation, = 0.000437 Fourier fit.

In comparison, the spherical-outer-boundary solution with ΔR₁ = 0 yields the values:

inner radius r_i(θ) = {4.34861, …, 4.3423}, mean = 4.343, dthrel = 0.00145.

total mass M02eff = 4.36731, inner boundary max(ρ_bc) = 4.68324.

shell thickness dR₁: mean = 0.0370, dthrel = 0.1736.

mean energy density = 1.75914.

error: med(err) = 0.0101 wavefront, = 0.000390 interpolation, = 0.000324 Fourier fit.

In contrast to the stellar shell-star, here the relative variation of the shell thickness for the spherical-outer-boundarysolution is smaller by a factor of 20 as compared to the minimal solution with a high outer ellipticity, so here there is a dependence of the shell thickness on the ellipticity.

The density distribution is shown in Figures 37-39.

The density distribution increases in th-direction.

The mass distribution is shown in Figure 40, Figure 41.

The physical mass distribution ends at the inner boundary at r_i = 4.46456, where the density jumps to ρ = 0. The fit extrapolates it to lower r-values.

The maximum distance from the horizon is max(r_02e) − r₊ = 0.125, therefore the minimal light energy attenuation is roughly 4.262/0.125 = 34, it means that visible green light of 0.514 μm is shifted to 17 μm into far-infrared.

Figure 37. Typical rotating galactic shell star: density over x = radius r₁, y = angle th, radius in mega-r_s(sun)=3 10⁶ km, th in radian.

Figure 38. Typical rotating galactic shell star: density over angle th at the inner boundary r₁ = 4.464, radius in mega-r_s(sun)=3 10⁶ km, th in radian.

Figure 39. Typical rotating galactic shell star: fitted density profile for th = 0.1 (equatorial), radius in mega-r_s(sun)=3 10⁶ km, th in radian.

Figure 40. Typical rotating galactic shell star: inner radius r_i over angle th, radius in mega-r_s(sun)=3 10⁶ km, th in radian.

The galactic shell-star has all its mass concentrated within a thin shell

(dR₁ = 0.0362) which has its inner radius inside and its outer radius outside its Schwarzschild-radius M₀ = 4.367, where the minimum distance from the horizon is

min(R₁(θ)) − M₀ = 0.0127,

therefore the light energy loss is approximately (in Newtonian approximation)

Figure 41. Typical rotating galactic shell star: total mass M₁(r₁) over radius, radius in r_s(sun)=3 km, th in radian.

Figure 42. Mass distribution of stellar black-holes and neutron stars measured by LIGO.

M₀/min(R₁(θ)) = 0.9971 and the attenuation factor 1/(1 − M₀/min(R₁(θ))) = 345, it means that x-ray-radiation from in-falling matter from the accretion disc with an energy of 5 keV and λ = 0.2 nm is shifted to λ = 69 nm, that is into hard UV-radiation.

9. Experimental Evidence with Recent LIGO and X-Ray Measurements

In November 2018, the LIGO cooperation published the newest statistics of neutron stars and black holes, based on gravitational waves and x-ray measurements [24].

The resulting mass distribution for black-holes and neutron stars is shown in Figure 42 [24].

From these results, we can deduce a confirmed mass range for neutron stars of $0.9M_{sun}\le M_{ns}\le 2.9M_{sun}$ and for stellar black holes $5M_{sun}\le M_{bh}\le 80M_{sun}$.

This mass range confirms our results from chapter 7 for compact neutron stars in the range $M\le 3.04M_{sun}$ , but not for the neutron shell stars in the range $3.04M_{sun}\le M\le 4.91M_{sun}$.

It means that the fluid (non-compact) shell stars are not stable, and only Fermi-gas shell stars are stable.

The mass range for stellar shell star black holes deduced in chapter 7:

$5.5M_{sun}\le M_{bh}\le 81.3M_{sun}$ is confirmed, considering that the uncertainty of mass measurement via x-rays is easily $\delta M_{bh}\ge 0.5M_{sun}$.

10. Conclusions

We introduce in chap. 6 an eos for the nucleon-fluid in the density range ρ_c ≤ ρ ≤ ρ_m, where ρ_c = 0.0417ρ_sand ρ_m = 0.0544ρ_s (sun units r_ss = 3 km, ρ_s = 1.76 × 10¹⁶ g/cm³), which is based on measurement data for the nucleon-nucleon-potential. This suggests, that there is a phase transition at ρ = ρ_c from the (interacting) nucleon fluid to the (weakly interacting) nucleon Fermi-gas.

Based on these 2 eos’s the results for the TOV-equation in chap. 7 are as follows.

Neutron stars obey the nucleon fluid eos and there are compact neutron stars in the range (M₀, R₁) = (0.14M_sun, 1.49r_ss), …, (3.04M_sun, 3.95r_ss), the R-M-relation follows approximately a cubic-root-law: R~M^1/3.

Neutron shell-stars in the range (M₀, R₁) = (3.04M_sun, 3.95r_ss), …, (4.91M_sun, 4.92r_ss) are not stable.

Stellar shell-stars exist in the range of (M₀, R₁) = (5.5M_sun, 9.1r_ss), …, (81.3M_sun, 91.2r_ss).

The underlying equation-of-state is the Fermi-gas of nucleons with the eos $P\left(\rho \right)=K_1{\rho }^{5/3}$. The resulting R-M-relation is practically linear and has a maximum mass value of M_max = 81.3M_sun. The light attenuation factor (redshift) is roughly {1.7, …, 50}. Taken the redshift and the small relative shell thickness of around 0.02, these stellar shell-stars have approximately the properties expected of a genuine black-hole, when measured from a distance $r\gg R_1$. Furthermore, the phase space volume of a thin spherical shell is proportional to its surface A, which approximates the Bekenstein-Hawking black-hole entropy formula $S=\left(k_B/L_P^2\right)A/4$.

The galactic (supermassive) shell-stars have the density scale and the eos of a white-dwarf-star, i.e. of an electron Fermi-gas. The R-M-relation is almost linear and goes from 1MM_sun up to 50MM_sun (MM_sun = 10⁶M_sun, Mr_ss = 10⁶r_ss). dR_rel = (R₁ – r_i)/M₀ is the relative thickness, and shows, that the shells are very thin indeed, with a minimum of 0.001. The relative Schwarzschild-distance dR_srel = (R₁ – M₀)/M₀ has a minimum at {M₀, dR_srel} = {7MM_sun, 0.00142857}, the redshift is around 700. So the overall result is, that the supermassive shell-stars become ever thinner shells, while the distance from the Schwarzschild-horizon is increasing.

In chap. 8 we present numerical results for rotating stars of the 3 types compact neutron star, stellar shell-star and galactic shell-star.

The angular velocity ω was chosen at ω = 0.36ω_max, i.e. about 1/3 of the maximum.

The compact neutron star with M₀ = 0.932M_sun, R_1y = 2.8372r_ss = 8.51 km, R_1x = 2.8391r_ss, ω = 0.1087, has the relative ellipticity of dthrel = 0.00118. The neutron star behaves like a fluid because of its “viscosity”, that is, its nuclear interaction, and becomes slightly “pumpkin-like”.

The stellar shell-star with M₀ = 15.74M_sun, R_1mean = 17.74r_ss, ω = 0.0126, has maximum density ρ_bc = 0.03595ρ_s, outer radii R_1y = 17.89r_ss = 53.67 km, R_1x = 17.59r_ss, inner radii r_iy = 16.98r_ss, r_ix = 16.70r_ss, inner rel. ellipticity dthrel = 0.0162. The redshift is 16.53.

The stellar shell-star behaves like a ball of neutron gas (negligible interaction) and decreases slightly its equatorial radius, so that, speaking naively, the increased gravitation counteracts the centrifugal force, the shell-star becomes “cigar-like”, with the shell thickness approximately constant.

The galactic shell-staris modelled (approximately) on the central black-hole in the Milky Way with mass M₀ = 4.368MM_sun(MM_sun = 10⁶M_sun, Mr_ss = 10⁶r_ss), radius R₁ = 4.38Mr_ss (=13.14 × 10⁶ km), ω = 0.05239.

It has maximum density ρ_bc = 4.961 × 10⁻¹²ρ_s, outer radii R_1y = 4.38Mr_ss, R_1x = 4.494Mr_ss, inner radii r_iy = 4.341Mr_ss, r_ix = 4.464Mr_ss, inner rel. ellipticity dthrel = 0.0276.

The redshift is roughly 345. The galactic shell-star is a shell object with a thin mass shell (ΔR = 0.0352Mr_ss) situated close above its outer Kerr horizon r₊ = 4.26Mr_ss. The polar radius is smaller than the equatorial radius, so the outer shape and the inner shape are both pancake-like.

The overall result is, that the introduction of numerical shell-star solutions of the TOV- and Kerr-Einstein-equations creates shell-star star models, which mimic closely the behaviour of abstract black holes and satisfy the Bekenstein_Hawking entropy formula, but have finite redshifts and escape velocity v < c, no singularity, no information loss paradox, and are classical objects, which need no recourse to quantum gravity to explain their behaviour.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Fliessbach, T. (1990) Allgemeine Relativitätstheorie. Bibliographisches Institut, Leipzig.
[2]	Lattimer, J.M. and Prakash, M. (2000) Neutron Star Structure.
[3]	Visser, M. (2008) The Kerr Spacetime: A Brief Introduction.
[4]	Kokkotas, K.D. and Vavoulidis, M. (2005) Rotating Relativistic Stars. Journal of Physics, Conference Series, 8, 71-80. https://doi.org/10.1088/1742-6596/8/1/009
[5]	Wasserman, A.L. (2011) Thermal Physics. Oregon State University, Corvallis. https://doi.org/10.1017/CBO9780511902611
[6]	Helm, J. (2022) Mathematica-Notebook GRSchwarzTOVOrig3.nb. https://www.researchgate.net/publication/358270854_Star_models_as_solutions_of_TOV-equation, calculation for non-rotating star-models
[7]	Papenbrock, Th. (2008) Physics of Nuclei. National Nuclear Physics Summer School.
[8]	Steiner, A.W., Hempel, M. and Fischer, T. (2013) Core-Collapse Supernova Equations of State Based on Neutron Star Observations. The Astrophysical Journal, 774, Article No. 17. https://doi.org/10.1088/0004-637X/774/1/17
[9]	Hempel, M. (2014) Materie am Limit. Physik in unserer Zeit, 45, 12-20. https://doi.org/10.1002/piuz.201301352
[10]	Hjorth-Jensen, M. (2007) Models for Nuclear Interactions. arxiv-nucl-th9811101.
[11]	Douchin, F. and Haensel, P. (2001) A Unified Equation of State of Dense Matter and Neutron Star Structure. Astronomy and Astrophysics, 380, 151-167.
[12]	Rosenfield, P. (2007) Properties of Rotating Neutron Stars. PhD Thesis, University of Washington, Seattle.
[13]	Urbanec, M. (2010) Equations of State and Structure of Neutron Stars. PhD Thesis, Silesian University in Opava, Opava.
[14]	Stergioulas, N. (2003) Rotating Stars in Relativity. Living Reviews in Relativity, 6, Article No. 3. https://doi.org/10.12942/lrr-2003-3
[15]	Dain, S. (2012) Geometric Inequalities for Axially Symmetric Black Holes. Classical and Quantum Gravity, 29, Article ID: 073001.
[16]	Typel, S., et al. (2010) Composition and Thermodynamics of Nuclear Matter with Light Clusters. Physical Review C, 81, Article ID: 015803. https://doi.org/10.1103/PhysRevC.81.015803
[17]	Chandrasekhar, S. (1992) The Mathematical Theory of Black Holes. Oxford University Press, Oxford.
[18]	Helm, J. (2022) Mathematica-Notebook KerrBLS2e.nb. https://www.researchgate.net/publication/358270862_Star_models_as_solutions_of_Kerr-spacetime_with_matter, calculation for rotating compact neutron-star
[19]	Helm, J. (2022) word-doc KerrBLS_Results.doc. https://www.researchgate.net/publication/358270862_Star_models_as_solutions_of_Kerr-spacetime_with_matter, concise results for rotating star-models
[20]	Helm, J. (2022) Mathematica-Notebook KerrBLS0e.nb. https://www.researchgate.net/publication/358270862_Star_models_as_solutions_of_Kerr-spacetime_with_matter, calculation for rotating shell-star
[21]	Ferrarese, L. and Merritt, D. (2002) Supermassive Black Holes. Physics World, 15, 41-46.
[22]	Helm, J. (2022) Mathematica-Notebook KerrBLS1e.nb. https://www.researchgate.net/publication/358270862_Star_models_as_solutions_of_Kerr-spacetime_with_matter, calculation for rotating galactic shell-star
[23]	Helm, J. (2022) Mathematica-Notebook KerrBLS1e.nb.
[24]	LIGO Cooperation (2018, November). https://www.ligo.caltech.edu/image/ligo20171016a