New Solutions of Gravitational Collapse in General Relativity and in the Newtonian Limit

Jan Helm

Technical University of Berlin, Berlin, Germany.

DOI: 10.4236/jhepgc.2022.82034                

Abstract

We discuss the Oppenheimer-Snyder-Datt (OSD) solution from a new perspective, introduce a completely new formulation of the problem exclusively in external Schwarzschild space-time (ESM) and present a new treatment of the singularities in this new formulation. We also give a new Newtonian approximation of the problem. Furthermore, we present new numerical solutions of the modified OSD-model and of the ball-to-ball-collapse with 4 different numerical methods.

Keywords

General Relativity, Gravitational Collapse, External Schwarzschild Space-Time

1. Introduction

The gravitational collapse (GC) is, together with the Robertson-Walker-Friedmann-Lemaitre (RWFL) cosmological models, the most important dynamical model in General Relativity. In its OSD form it has, apart from RWFL, the only closed analytic solution in this area. In astrophysics and cosmology it is of utmost importance, because it is considered to be the valid model for the formation of stars from dust and gas clouds and for catastrophic events like star collapse to a neutron star or a Black Hole. However, the OSD formulation has severe drawbacks: its assumptions of homogeneous density and zero pressure are completely unrealistic and especially the latter can be even regarded as unphysical (Mitra [1] ). Furthermore, the formation of the Black Hole asserted in OSD happens at infinite time for an external observer in the corresponding external Schwarzschild or Vaidya space-time (ESM) (see e.g. [2] ). Finally, OSD uses co-moving coordinate frame as the formulation basis and introduces junction conditions (continuity of space-time function and derivative) to the external observer frame in ESM, which are difficult to calculate symbolically and to implement numerically.

Recent publications on this subject can be divided into 7 categories:

· observational astrophysics with application of GC-formalism Kotake [3];

· review Naidu [4] Joshi [5] Lasky [6];

· gravitational collapse of fluid Lasky [6];

· gravitational collapse in star formation in Newtonian formulation Girichidis [7];

· extensions of OSD in COF (co-moving frame) with respect to heat-flux Herrera ( [8] [9] ) Goswami [10], radiation Sharma [11], neutrino-emission Nakazato [12], equation-of-state Sanwe [13] Joshi [14];

· singularity freedom in gravitational collapse Marshall [15], Mitra [1];

· gravitational collapse quantization from quantum-theoretical point of view Hajicek [16], Corda [17].

In this article we introduce the essential feature system equations for GC exclusively in ESM, which adds 1 differential equation to the 4 of OSD, and radial velocity as a new, fifth variable function.

This makes the system more complicated, particularly a separation ansatz in r and t like in OSD is not possible anymore, but there are no junction conditions and all results are directly observable for an external observer.

The article can be subdivided into three parts:

· fundamentals and the system equations for the different models in chapters 2, 3, 4;

· the particular models and discussion of singularities in chapters 5, 6, 7;

· numerical methods and results in chapters 8, 9.

2. GR Fundamentals

The most general spherically symmetric line element, in spherical coordinates (x^μ) = (t, r, θ, f), can be written as

$\text{d}{s}^2=-A{\left(t,r\right)}^2{c}^2\text{d}{t}^2+B{\left(t,r\right)}^2\text{d}{r}^2+Y{\left(t,r\right)}^2\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)$ (1)

where A, B and Y are functions of the coordinates t and r.

The Einstein field equations with the above line element are:

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

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)$ (3)

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 }$ (4)

The Einstein tensor is given by

$G_{\mu \nu }=R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }{R}_0$ (5)

we obtain from (3) and (4) the non-vanishing Einstein tensor components [4] [8],

$G_{00}=\left(2\frac{B°}{B}\frac{Y°}{Y}+\frac{{\left(Y°\right)}^2}{Y^2}\right)-\frac{A^2}{B^2}\left(-2\frac{B\text{'}}{B}\frac{Y\text{'}}{Y}+\frac{{\left(Y\text{'}\right)}^2}{Y^2}+2\frac{Y\text{'}\text{'}}{Y}\right)+\frac{A^2}{Y^2}$ (6)

$G_{11}=\frac{B^2}{Y^2}\left(2\frac{A°}{A}\frac{Y°}{Y}-\frac{{\left(Y°\right)}^2}{Y^2}-2\frac{Y°°}{Y}\right)+\left(2\frac{A\text{'}}{A}\frac{Y\text{'}}{Y}+\frac{{\left(Y\text{'}\right)}^2}{Y^2}\right)-\frac{B^2}{Y^2}$ (7)

$G_{01}=2\left(\frac{B°}{B}\frac{Y\text{'}}{Y}+\frac{A\text{'}}{A}\frac{Y°}{Y}-\frac{Y°\text{'}}{Y}\right)$ (8)

$\begin{array}{l}{G}_{22}=-\frac{Y^2}{A^2}\left(\frac{B°°}{B}-\frac{A°}{A}\frac{B°}{B}+\frac{B°}{B}\frac{Y°}{Y}-\frac{A°}{A}\frac{Y°}{Y}+\frac{Y°°}{Y}\right)\\ -\frac{Y^2}{B^2}\left(\frac{A\text{'}\text{'}}{A}-\frac{A\text{'}}{A}\frac{B\text{'}}{B}+\frac{A\text{'}}{A}\frac{Y\text{'}}{Y}-\frac{B\text{'}}{B}\frac{Y\text{'}}{Y}+\frac{Y\text{'}\text{'}}{Y}\right)\end{array}$ (9)

$G_{33}={\sin }^2\theta G_{22}$ (10)

The dot $°$ represents ${\partial }_t$ and the prime $\text{'}$ represents ${\partial }_r$.

The energy-momentum tensor T for a perfect fluid with 4-velocity u^μ, heat-flow q^μ, density ρ, and pressure P is

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

We will consider here only adiabatic systems without heat-flow to outside, so in the following q = 0.

For the gravitation collapse problem two coordinate systems are used.

In the co-moving coordinate system (COF) the spatial part of the velocity is zero:

$u^{\mu }=\frac{1}{A}{\delta }_0^{\mu }$

and the acceleration 4-vector is

$u{°}^{\mu }=\frac{A\text{'}}{A}{\delta }_1^{\mu }$

In the external Schwarzschild space-time (ESM) with Schwarzschild-radius r_s

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

we take the coordinate system of the non-moving observer at infinity.

In order to make the variables dimensionless, we introduce “sun units”, like in [18], where r_s means the Schwarzschild-radius.

In the following, we make all equations and variables dimensionless, by using the “proper units” of the gravitating system $r_s=r_{ss}\frac{M}{M_{sun}}$ Schwarzschild-radius for r, M_sun for M, ${\rho }_s={\rho }_{ss}/{\left(\frac{M}{M_{sun}}\right)}^2$ for ρ, $P_s={\rho }_s{c}^2$ for P.

3. System Equations for COF

We insert the above expression for u^μ into T and get the Einstein-equations as the system equations:

eq00: $\begin{array}{c}{c}_0\rho =\exp \left(-2fA\right)\left(2fB°fY°+{\left(fY°\right)}^2\right)\\ -\exp \left(-2fB\right)\left(2fY\text{'}\text{'}+3{\left(fY\text{'}\right)}^2-2fB\text{'}fY\text{'}\right)+\exp \left(-2fY\right)\end{array}$ (13)

eq11: $\begin{array}{c}{c}_{0}P=\exp \left(-2fA\right)\left(2fA°fY°-3{\left(fY°\right)}^2-fY°°\right)\\ -\exp \left(-2fB\right)\left(2fA\text{'}fY\text{'}+{\left(fY\text{'}\right)}^2\right)-\exp \left(-2fY\right)\end{array}$ (14)

eq01: $-c_{0}q=0=2\exp \left(-2fB-fA\right)\left(fB°fY\text{'}+fA\text{'}fY°-fY°\text{'}+fY°fY\text{'}\right)$ (15)

eq22: $\begin{array}{c}{c}_{0}P=-\exp \left(-2fA\right)\left(fB°°+{\left(fB°\right)}^2-fA°fB°+fB°fY°+fY°°+{\left(fY°\right)}^2\right)\\ +\exp \left(-2fB\right)\left(fA\text{'}\text{'}+{\left(fA\text{'}\right)}^2-fA\text{'}fB\text{'}+fA\text{'}fY\text{'}-fB\text{'}fY\text{'}+fY\text{'}\text{'}+{\left(fY\text{'}\right)}^2\right)\end{array}$ (16)

where $A\left(t,r\right)=\exp \left(fA\left(t,r\right)\right)$, $B\left(t,r\right)=\exp \left(fB\left(t,r\right)\right)$ , $Y\left(t,r\right)=\exp \left(fY\left(t,r\right)\right)$ .

The constant $c_0=\frac{8\pi G}{c^2}=\frac{3}{{\left(r_{ss}\right)}^2{\rho }_s}$, explicitly and in sun-units.

In the dimensionless representation and M = M_sun:c₀= 3.

When calculating in COF, one uses in parallel an external space-time valid outside the gravitating mass system boundary, the usual choice is the Vaidya space-time for a radiating spherically symmetric mass system with total mass M(t)and initial massM₀ in new coordinates r₂, t₂:

$\text{d}{s}^2=-\left(1-\frac{r_s}{r_2}\frac{M\left(t_2\right)}{M_0}\right)c^2\text{d}{t}_2^2-2\text{d}{t}_2\text{d}{r}_2+r_2^2\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)$ (17)

When the system is adiabatic, i.e. heat-flow = 0 and M(t)= M₀, the Vaidya space-time becomes the Schwarzschild space-time in Eddington-Finkelstein coordinates:

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

The line element functions of the inner and the external space-time must be continuous, and so must be their derivatives normal to the surface boundary, this generates junction conditions on the surface [8]: with the boundary surface space-time $\text{d}{s}^2=-c^2\text{d}{\tau }^2+Y^2\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)$

$A\text{d}t=\text{d}{t}_2\left(1-\frac{r_s}{r_2}\frac{M\left(t_2\right)}{M_0}\right)=\text{d}\tau$ (19)

$Y=r_2\left(t_2\right)$ (20)

$A\frac{\text{d}{t}_2}{\text{d}\tau }=1-\frac{r_s}{r_2}\frac{M\left(t_2\right)}{M_0}+2\frac{\text{d}{r}_2}{\text{d}{t}_2}$ (21)

$\begin{array}{l}-2\left(\frac{B°}{B}\frac{Y\text{'}}{Y}+\frac{A\text{'}}{A}\frac{Y°}{Y}-\frac{Y°\text{'}}{Y}\right)\\ =-\frac{B}{A}\left(2\frac{Y°°}{Y}-2\frac{A°}{A}\frac{Y°}{Y}+2{\left(\frac{Y°}{Y}\right)}^2\right)+\frac{A}{B}\left(2\frac{A\text{'}}{A}\frac{Y\text{'}}{Y}+{\left(\frac{Y\text{'}}{Y}\right)}^2-\frac{B}{Y^2}\right)\end{array}$ (22)

In the last equation, the dot $°$ represents ${\partial }_t$ and the prime $\text{'}$ represents ${\partial }_r$.

4. System Equations for Schwarzschild Space-Time ESM

Here we give up the COF, and calculate exclusively in ESM. Now the radial velocity u(t,r) becomes a new variable function and we derive an additional system equation for it from the GR orbit equations, where τ is the proper time:

$\frac{{\text{d}}^2{x}^{\kappa }}{\text{d}{\tau }^2}=-{\Gamma }_{\mu \nu }^{\kappa }\frac{\text{d}{x}^{\mu }}{\text{d}\tau }\frac{\text{d}{x}^{\nu }}{\text{d}\tau }$ (23)

and the relativistic velocity normalization condition (dimensionless, c = 1)

$-1=g_{\mu \nu }\frac{\text{d}{x}^{\mu }}{\text{d}\tau }\frac{\text{d}{x}^{\nu }}{\text{d}\tau }$ (24)

we get the following differential equations for t = x⁰ and r = x¹ and their derivative $\text{'}$ for τ:

$t\text{'}\text{'}=-{\Gamma }_{00}^{0}t{\text{'}}^2-2{\Gamma }_{01}^{0}t\text{'}r\text{'}-{\Gamma }_{11}^{0}r{\text{'}}^2=-\frac{A^{\left(1,0\right)}}{A}t{\text{'}}^2-2\frac{A^{\left(0,1\right)}}{A}t\text{'}r\text{'}-\frac{B{B}^{\left(1,0\right)}}{A^2}t{\text{'}}^2$

$r\text{'}\text{'}=-{\Gamma }_{00}^{1}t{\text{'}}^2-2{\Gamma }_{01}^{1}t\text{'}r\text{'}-{\Gamma }_{11}^{1}r{\text{'}}^2=-\frac{A\text{}{A}^{\left(0.1\right)}}{B^2}t{\text{'}}^2-2\frac{B^{\left(1,0\right)}}{B}t\text{'}r\text{'}-\frac{B^{\left(0,1\right)}}{B}t{\text{'}}^2$

$t\text{'}=\frac{\sqrt{1+r{\text{'}}^2{B}^2}}{A}$

And from that for the radial velocity u_r= r':

$u_r\text{'}=-\frac{\text{}{A}^{\left(0,1\right)}}{A}\left(u_r^2+\frac{1}{B^2}\right)-2\frac{B^{\left(1,0\right)}}{A}{u}_r\sqrt{u_r^2+\frac{1}{B^2}}-\frac{B^{\left(0,1\right)}}{B}{u}_r^2$

We get for the energy-momentum tensor

$u_0=A\sqrt{u_r^2{B}^2+1}$

$u_1=u_r{B}^2$

$T_{00}=\left(\rho +P\right)\left(u_r^2{B}^2+1\right)A^2$

$T_{11}=-\left(\rho +P\right)B^4{u}_r^2+P{B}^2$

$T_{22}=P{Y}^2$

$T_{33}=P{Y}^2{\sin }^2\theta$

From this and the expressions for G_μν (6…10) get the system equations for ESM, where u(t,r)= u_r

eq00: $\begin{array}{l}{c}_0\left(\rho \left(1+\exp \left(fu+fB\right)\right)+P\left(\exp \left(fu+fB-fA\right)\right)\right)\\ =\exp \left(-2fA\right)\left(2fB°fY°+{\left(fY°\right)}^2\right)\\ -\exp \left(-2fB\right)\left(2fY\text{'}\text{'}+3{\left(fY\text{'}\right)}^2-2fB\text{'}fY\text{'}\right)+\exp \left(-2fY\right)\end{array}$ (25)

eq11: $\begin{array}{l}{c}_0\left(P\left(1-\exp \left(fu+fB\right)\right)-\rho \exp \left(fu+fB\right)\right)\\ =\exp \left(-2fA\right)\left(2fA°fY°-3{\left(fY°\right)}^2-fY°°\right)\\ -\exp \left(-2fB\right)\left(2fA\text{'}fY\text{'}+{\left(fY\text{'}\right)}^2\right)-\exp \left(-2fY\right)\end{array}$ (26)

eq01: $-c_{0}q=0=2\exp \left(-2fB-fA\right)\left(fB°fY\text{'}+fA\text{'}fY°-fY°\text{'}+fY°fY\text{'}\right)$ (27)

eq22: $\begin{array}{l}{c}_{0}P=-\exp \left(-2fA\right)\left(fB°°+{\left(fB°\right)}^2-fA°fB°+fB°fY°+fY°°+{\left(fY°\right)}^2\right)\\ +\exp \left(-2fB\right)\left(fA\text{'}\text{'}+{\left(fA\text{'}\right)}^2-fA\text{'}fB\text{'}+fA\text{'}fY\text{'}-fB\text{'}fY\text{'}+fY\text{'}\text{'}+{\left(fY\text{'}\right)}^2\right)\end{array}$ (28)

eq44: $\begin{array}{l}fu°\exp \left(fA\right)+fu\text{'}\exp \left(fB\right)\\ =-fA\text{'}\exp \left(fu\right)\left(1+\exp \left(-2fB-2fu\right)\right)\\ -2fB°\exp \left(fB-fA+fu\right)\sqrt{\left(1+\exp \left(-2fB-2fu\right)\right)}-fB\text{'}\exp \left(fu\right)\end{array}$ (29)

where $A\left(t,r\right)=\exp \left(fA\left(t,r\right)\right)$, $B\left(t,r\right)=\exp \left(fB\left(t,r\right)\right)$ , $Y\left(t,r\right)=\exp \left(fY\left(t,r\right)\right)$, $u\left(t,r\right)=\exp \left(fu\left(t,r\right)\right)$.

5. The Gravitational Dust Cloud Collapse Model of Oppenheimer-Snyder-Datt

The Oppenheimer-Snyder-Datt (OSD) model of a gravitational collapse was the first exact GR-solution of a dynamic GR system, Joshi [5]. With the (unrealistic) assumption of a no-pressure dust cloud with homogeneous density and calculating in COF one can use a separation-ansatz in the system equations and gets an explicit analytic expression for the function variables.

In the following we use the terminology in [2] and [5].

In COF the solution space-time becomes (dimensionless)

$\text{d}{s}^2=-\text{d}{t}^2+\frac{v{\left(t\right)}^2}{1-k{r}^2}\text{d}{r}^2+v{\left(t\right)}^2{r}^2\left(\text{d}{\theta }^2+{\sin }^2\theta \text{d}\phi \right)$ (30)

where the edge of the cloud is in COF

$R_1\left(t\right)=r_{0}v(\; t\; )$

and $k=1/r_0^2$, and r₀ is the initial edge radius, v(0)= 1, density $\rho \left(t\right)={\rho }_0/v{\left(t\right)}^3$.

The function v(t) = 0 at the dimensionless time $T_{col}=\frac{\pi }{2}{r}_0^{3/2}$ (see [2] ), where the radius of the cloud R₁(T_col)= 0 and there is a singularity with density blowing-up.

In ESF, with coordinates t₂ and r₂, the (dimensionless) equation for the edge is ( [2] 45.10):

$\frac{\text{d}{t}_2}{\text{d}{r}_2}=-\frac{1}{\sqrt{\frac{1}{r_2}}\left(1-\frac{1}{r_2}\right)}$

which gives the relationship

$t_2=t_{20}-\left(2\sqrt{r_2}\left(1+\frac{r_2}{3}\right)+\log \left(\frac{\sqrt{r_2}-1}{\sqrt{r_2}+1}\right)\right)$ (31)

$t_2\to \infty ,r_2=1+const\ast \exp \left(-t_2\right)$

so in the limit of infinite time t₂ the edge in ESM goes exponentially slowly (with the characteristic time t_s= r_s/c) to the Schwarzschild-radius r_s: there is no singularity, as the Black Hole forms at infinity.

This result can be found in [2] [15] and [19].

In [19] the author gives a proof for non-existence of singularities in gravitational collapse in COF under certain conditions.

6. Gravitational Collapse in the Newtonian Approximation

In the Newtonian approximation, we consider the total energy density of pressure, kinetic energy and gravitational energy,

$E_d=P_s\left(\rho \right)+\frac{\rho v^2}{2}-\frac{GM\left(r\right)\rho }{r}$ (32)

with the equation-of-state (eos) of the ideal gas $P_s\left(\rho \right)=\rho T_r$ which is conserved and therefore stationary under ${\partial }_t$ and ${\partial }_r$, which gives 2 differential equations for the 2 function variables ρ(t,r) and v(t,r).

$M\left(r\right)=4\pi {\displaystyle \underset{0}{\overset{r}{\int }}\rho \left(r_x\right)r_x^2\text{d}{r}_x}$ is the mass within the radius r.

In proper units, i.e. dimensionless the total energy density becomes

$E_d=\rho T_r+\frac{\rho v^2}{2}-\frac{3\rho }{2r}{\displaystyle \underset{0}{\overset{r}{\int }}\rho \left(r_x\right)r_x^2\text{d}{r}_x}$

Into the eos enters the relative temperature $T_r=\frac{kT}{m_n{c}^2}$, relative to the nucleon mass.

The resulting (dimensionless) equations are of order 1 in ρ(t,r) and v(t,r).

equt: ${\partial }_t{E}_d=0$

$\frac{\rho °}{\rho }=\frac{\left(v°v-\frac{3\rho }{2r}I\left[r^2\rho °\right]\right)}{\frac{3\rho }{2r}I\left[r^2\rho \right]-\left(T_r+\frac{v^2}{2}\right)}$ (33)

equr: ${\partial }_r{E}_d=0$

$\frac{\rho \text{'}}{\rho }=\frac{\left(v\text{'}v+\frac{3\rho }{2r^2}I\left[r^2\rho \right]-\frac{3{\rho }^{2}r}{2}\right)}{\frac{3\rho }{2r}I\left[r^2\rho \right]-\left(T_r+\frac{v^2}{2}\right)}$ (34)

with the abbreviation

$I\left[r^2\rho \right]={\displaystyle \underset{0}{\overset{r}{\int }}\rho \left(r_x\right)r_x^2\text{d}{r}_x}$ and $I\left[r^2\rho °\right]={\displaystyle \underset{0}{\overset{r}{\int }}\rho °\left(r_x\right)r_x^2\text{d}{r}_x}$

and mean density ${\rho }_m=\frac{I\left[r^2\rho \right]}{r^2}.$

These equations can be simplified by using the integrability condition:

${\partial }_t$ equr = ${\partial }_r$ equt

$\frac{3\rho }{2r}I\left[r^2\rho \right]-\left(T_r+\frac{v^2}{2}\right)=\frac{3{\rho }^2}{2r\rho °}\left(\frac{I\left[r^2\rho °\right]}{r^2}-1\right)$ (35)

With this the Equations (33), (34) become

$\rho {°}_m\left(1+r^2\right)=1+\frac{1}{3}\frac{r}{\rho }\left(v^2\right)°$ (33a)

$\rho \text{'}\left(1-\rho {°}_m\right)=\rho °\left(\rho r^2-r{\rho }_m-\frac{1}{2\rho }\left(v^2\right)\text{'}\right)$ (34a)

Using (33a) and $v^2<1$ one can show that ρ cannot blow-up at r = 0, so there is no singularity.

The numerical results for the ball-to-ball model are given in 9.6 below.

7. Singularity in ESM Gravitational Collapse

In this section we examine the singularity behavior in Schwarzschild space-time gravitational collapse.

We assume that there is a singularity in one or several of the function variable at the time t = t₀, i.e. A(t₀) or B(t₀) or Y(t₀)= 0. For simplicity, by a coordinate shift we can make t0= 0.

We make the ansatz $B\left(t,r\right)=t^k\exp \left(fBr\left(r\right)\right),t\to 0$, $Y\left(t,r\right)=t^k\exp \left(fYr\left(r\right)\right),t\to 0$,

$A\left(t,r\right)=\exp \left(fYr\left(r\right)\right),t\to 0$

that is, B and Y vanish at t = 0 like t^k(k ≥2) and A remains non-zero, this behavior is modeled at the OSD-COF singularity (where k = 2), but any other choice of vanishing space-time function with k = 2 will bring the same result.

For the logarithmic function variables it means $fB\left(t,r\right)=k\log \left(t\right)+fBr\left(r\right),t\to 0$, $fY\left(t,r\right)=k\log \left(t\right)+fYr\left(r\right),t\to 0$ both B and Y have a genuine logarithmic singularity at 0, where they become negative infinite.

Before the singularity is reached, all physical feasibility conditions must be met, i.e. in particular $0\le u\left[t,r\right]<1$, $v_{sound}=\sqrt{\frac{\partial P}{\partial \rho }}=k_1\gamma {\rho }^{\gamma -1}<1$, mass condition $M_0={\displaystyle \underset{0}{\overset{r_0}{\int }}\rho \left[t,r\right]4\pi r^2\text{d}r}$, so with pressure non-zero, there is a physical limit for the density.

Under these assumptions and taking only terms, which have the strongest divergence in t, the ESM Equations (25)…(29) take the form

eq00: $0=-c_0\rho -\exp \left(-2fB\right)\left(2fYr\text{'}\text{'}+3{\left(fYr\text{'}\right)}^2-2fBr\text{'}fYr\text{'}-\exp \left(2fBr-2fYr\right)\right)$ (25s)

eq11: $0=-c_0{k}_1{\rho }^{\gamma }-\exp \left(-2fB\right)\left(2fAr\text{'}fYr\text{'}+{\left(fYr\text{'}\right)}^2-\exp \left(2fBr-2fYr\right)\right)$ (26s)

eq01: $0=fBt°\left(fYr\text{'}+fAr\text{'}-fYr\text{'}\right)=k{t}^{k-1}\left(fYr\text{'}+fAr\text{'}-fYr\text{'}\right)$ (27s)

eq22: $\begin{array}{l}0=-c_0{k}_1{\rho }^{\gamma }+\exp \left(-2fB\right)(fAr\text{'}\text{'}+{\left(fAr\text{'}\right)}^2-fAr\text{'}fBr\text{'}+fAr\text{'}fYr\text{'}\\ -fBr\text{'}fYr\text{'}+fYr\text{'}\text{'}+{\left(fYr\text{'}\right)}^2)\end{array}$ (28s)

eq44: $0=-fu°\exp \left(fA\right)-2fB°\exp \left(-fA\right)$ (29s)

where $A\left(t,r\right)=\exp \left(fAr\left(r\right)\right)$, $B\left(t,r\right)=\exp \left(k\ast \log \left(t\right)+fBr\left(r\right)\right)$ , $Y\left(t,r\right)=\exp \left(k\ast \log \left(t\right)+frY\left(r\right)\right)$, $u\left(t,r\right)=\exp \left(fut\left(t\right)+fur\left(r\right)\right)$ , $\rho \left(t,r\right)=rhot\left(t\right)\ast rhor\left(r\right)$ : with t-limit behavior a product ansatz in t and r is always possible, where ρ is the density and u is the velocity.

Now, as the density ρ must not blow-up (this would mean an unphysical solution before the Black Hole generation), the brackets with $\exp \left(-2fB\right)$ must vanish identically.

This gives $fAr\text{'}=0$, then $fBr\text{'}fYr\text{'}=0$, and finally $fBr\text{'}=0=fYr\text{'}$, and from (29s) also $fur=const$ so that the space-time and the velocity is spatially constant, which is impossible.

One can show, that other choices of divergent space-time functions lead to a blow-up of the density or to unphysical results.

In essence, this confirms in ESM the results of [19], that (in COF) the gravitational collapse singularity can arise only with zero pressure, which is physically untenable.

8. The Differential Equations and Their Numerical Form

The 5 differential equations for the gravitational collapse in SMF (Schwarzschild metric frame) have the form (25…29).

The first 4 are the Einstein-equations for R⁰⁰, R¹¹, R⁰¹, R²², and the last is the equation for the radial velocity u from the GR equations-of-motion. The dot represents ${\partial }_{t}Ak$ and the prime’ represents ${\partial }_{r}Ak$. The constant

$c_0=\frac{8\pi G}{c^2}=\frac{3}{{\left(r_{ss}\right)}^2{\rho }_s}$, explicitly and in sun-units, and P is replaced by the equation-of-state $P\left(\rho \right)=f_{eos}\left(\rho \right)$, the total mass M₀= 1(M_sun).

The chosen space-time range is with edge radius r₀ and final time T₁: t = 0...50(r_s/c),r = 0…10(r_s), r_s = Schwarzschild-radius.

The variable functions are $Ak=\left\{fA\left(t,r\right),fB\left(t,r\right),fY\left(t,r\right),\rho \left(t,r\right),fu\left(t,r\right)\right\}$, where $A\left(t,r\right)=\exp \left(fA\left(t,r\right)\right)$, $B\left(t,r\right)=\exp \left(fB\left(t,r\right)\right)$ , $Y\left(t,r\right)=\exp \left(fY\left(t,r\right)\right)$, $u\left(t,r\right)=\exp \left(fu\left(t,r\right)\right).$

It is important to note the highest derivatives of the variable functions, they are:

$fA°,fA\text{'}\text{'},fB°°,fB\text{'},fY°°,fY\text{'}\text{'},\rho ,fu°,fu\text{'}$

For the variables with the differential degree 2 in t, the t-boundary-condition becomes

$t=0$ : $Ak=fbx\left[r\right]$, ${\partial }_{t}Ak=fDbx\left[r\right]$,

for degree 1 it becomes

$t=0$ : $Ak=fbx\left[r\right]$,

and correspondingly in r.

The density ρ is the only algebraic variable function, so no boundary condition is imposed, its boundary value is calculated from the equations at the boundary, together with the other highest derivatives.

All numerical solution methods used here operate on an equidistant 2-dimensional lattice {t_i, r_j}with step size h₁= r₀/ndimx in t, and step size h₂= T₁/ndimy in r, ndim= 8…32, depending on the required execution time.

For Ritz-Galerkin test functions and for fits a 2-dimensional trigonometric (Fourier) expansion ${\displaystyle \underset{k,l=0}{\overset{ngrad}{\sum }}{a}_{k,l}\exp \left(Ik{\omega }_{t}t\right)\exp \left(Il{\omega }_{r}r\right)}$, with basic angular frequencies ω_t ω_r and $I=\sqrt{-1}$, of degree ngrad is used.

Before the solution procedure starts, an appropriate “seed function” w_varfit is chosen with the desired boundary conditions, and then corrected to solve the equations at the boundary, the boundary conditions are corrected appropriately. Its error is evaluated and a global trigonometric t-r-fit on the lattice {t_i,r_j}is calculated.

9. Numerical Solution

9.1. Numerical Solution Methods

Table 1 Ritz-Galerkin global minimization with a trigonometric test function and the parameter vector qv2t operates on the differential-equation-set with the variable functions replaced by the test function, see Helm [19].

Table 1. Ritz-Galerkin global minimization.

The goal function is the sum of the absolute error values on the lattice {t_i,r_j}.

First, an appropriate seed function w_varfit is chosen, which fulfils the equations at the boundary.

Its error is evaluated and a global t-r-fit on the lattice {t_i,r_j}to the Ritz-Galerkin-test function is performed with a resulting parameter value vector qv20t.

The boundary conditions

$t=0$ : $Ak=fbx\left[r\right]$, ${\partial }_{t}Ak=fDbx\left[r\right]$

$r=r_0$ : $Ak=fby\left[t\right]$, ${\partial }_{r}Ak=fDby\left[t\right]$

are enforced via a penalty added to the goal function.

The physical feasibility conditions

$0\le \rho \left[t,r\right]\le {\rho }_{\max }$, $0\le u\left[t,r\right]\le 1$, $M_0={\displaystyle \underset{0}{\overset{r_0}{\int }}\rho \left[t,r\right]4\pi r^2\text{d}r}$

are used as minimization precondition for the minimization procedure, the start vector for qv2t is qv20t, the initial fit parameter vector.

The result of the procedure is the global parameter vector qv2t, parameter replacement in the trigonometric test functions gives the solution functions.

The variable vector qv2t has the length 5 *ngrad².

Table 2 Finite-difference-minimization is based on the finite-difference discretization of the differential equations, see Helm [19].

In the finite-difference scheme, the variables are the values Ak(t_i,r_j) on the lattice {t_i,r_j}.

The differences in the scheme are t-forward and r-backward, because it starts with the boundary condition: the start-values in {t_i} and {r_j} are replaced by the boundary values.

The goal function is the sum of the absolute error values on the lattice {t_i,r_j}.

The start vector are the values {Ak(t_i,r_j)(w_varfit)} of the seed-function.

The physical feasibility conditions are imposed via a penalty function added to the goal function.

Finally, a global minimization with the variable-vector {Ak(t_i,r_j)} and the start vector is performed, yielding the array {Ak(t_i,r_j)}₀ as the solution.

The variable vector {Ak(t_i,r_j)} has the length 5*ndimx*ndimy, much more than global Ritz-Galerkin, but the minimization has no precondition and is therefore much faster.

Table 3 R-profile-wavefront for fixed t = t_i, calculates r-profiles for a fixed t_i, starting with the boundary t = t₀= 0. For each profile, the highest derivatives are

Table 2. Finite-difference-minimization.

Table 3. R-profile-wavefront.

calculated from the differential equations, the values of the lower derivatives are calculated from the predecessor profile via Euler step and inserted into the equations, see Helm [19].

The Euler-step is:

$A{k}^{\left(d,0\right)}\left(t_{i+1},r_j\right)=A{k}^{\left(d,0\right)}\left(t_i,r_j\right)+A{k}^{\left(d+1,0\right)}\left(t_i,r_j\right)\ast hstep$,

where $A{k}^{\left(d,0\right)}$ is the d-th derivative in t.

After the insertion of the lower derivatives the equations become algebraic equations for the highest-derivative-profiles $A{k}^{\left(d,0\right)}\left(r,t_i\right)$.

These functional algebraic equations in r can be solved point-wise consecutively starting with the r-boundary r = r_ndimy = r₀ and obeying continuity requirements (discrete algebraic profile calculation).

Alternatively, they can be solved by 1-dimensional trigonometric Ritz-Galerkin minimization starting with a 1-dimensional trigonometric r-fit to the predecessor profile as start parameter-vector, and the r-boundary values as precondition (trigonometric Ritz-Galerkin r-profile calculation).

Here, we used the trigonometric Ritz-Galerkin r-profile calculation, because it is faster and the continuity is guaranteed automatically.

The initial step i = 0t = t₀ uses the boundary r-profile for t = t₀ as the predecessor.

The result of the wavefront procedure is the array of the solution function values on the lattice vla (value local array).

For each profile, the variable vector qvt has the length ngrad*5, so the total number of variables is ndimx*ngrad*5.

Finite-difference-wavefront for fixed t = t_i, calculates discrete r-profiles for a fixed t_i, starting with the boundary t = t₀= 0. Here, the differential equations are discretized via a finite-difference-scheme with variables Ak(t_i,r_j) on the lattice {t_i,r_j}. For each profile t = t_i, the discrete variables {Ak(t_i,r_j)}_j are calculated by minimization of the discretized equations, with start vector from the predecessor profile.

The initial step i = 0t = t₀ uses the boundary r-profile for t = t₀ as the predecessor.

The result of the wavefront procedure is the array of the solution function values on the lattice vla (value local array).

For each profile, the variable vector ${\left\{A{k}^{\left(d,0\right)}\left(t_i,r_j\right)\right\}}_j$ has the length ndimy*5, so the total number of variables is ndimx*ndimy*5.

9.2. Numerical Result Statistics

All calculations were made in Mathematica-code Helm [19], which is available for download and can be viewed either with Mathematica or with the freely available Wolfram CDF player.

We present the results for the function variables as 3D-plots and for the rest in form of a table, which contains the essential numerical parameters: lattice size, execution time on a conventional 3 GHz desktop and the error statistics.

For the minimization methods Ritz-Galerkin and finite-difference the numerical parameters are: lattice size, execution time in sec, error = median error of the solution, error-testfunc = median error of the seed functions, total-error = total minimization error, for comparison testf-error = total error of the seed functions, and the mass = numerical mass of the solution over time (should be 1). The minimization is carried out globally (global-min) and with start-vector = seed functions (local-min).

For the wavefront methods r-profile and finite-difference the numerical parameters are: lattice size, execution time in sec, error-fit = median error of the fit (spline, Fourier or polynomial interpolation) to the discrete solution values on the lattice, error-testfunc = median error of the seed functions, error-discrete = median error of the discrete solution array vla, and the mass = numerical mass of the solution over time (should be 1).

In the 3D-plots the first coordinate is time t = 0...50(r_s/c) and the second coordinate is radius r = 0…10(r_s), r_s = Schwarzschild-radius.

9.3. The Equation-of-State

The equation-of-state, which we are using here for the GR-gravitational-collapse is $P=k1{\rho }^{\gamma }$ for the non-interacting nucleon Fermi gas, with γ = 5/3 near the critical nucleon density.

In the limit of the Newtonian approximation, we use, as customary for sun-like stars, simply the dimensionless ideal-gas-equation (c = 1) $P=\rho T_r$, where $T_r=\frac{kT}{m_n{c}^2}$ is the relative temperature, relative to the nucleon mass.

9.4. Extended Oppenheimer-Snyder-Datt Model

The original OSD postulates dust without pressure (P = 0) and a homogeneous density (ρ = ρ(t)) in COF(co-moving frame), which results in a collapsing ball with a diminishing edge R₁(t), in ESM(external Schwarzschild metric) $R_1\left(t\right)=r_s\left(1+const\exp \left(-ct/r_s\right)\right),t\to \infty$, i.e. the edge reaches the Schwarzschild radius r_s in an infinite time.

We extend here OSD to incorporate $P=k1{\rho }^{\gamma }$ and the density initially constant with a smooth (exponentially zero) edge and we calculate exclusively in ESM.

We chose as the seed-functions the solution of OSD in COF transformed into ESM and smoothed. Still, the derivatives of the seed-functions are numerically not very well behaved, therefore the results of the 32 × 32-lattice are modest (error Ritz-Galerkin minimization 0.50 (Table 1) and for r-profile-wavefront 0.64 (Table 3).

The best results achieve the finite-difference minimization with an acceptable error 0.035 (Table 2).

The most interesting result is the density: it goes from an approximately constant distribution in r through a “build-up” at the center to an approximately constant distribution again at the final time, there is no blow-up.

Plot of numerical results:

time-radius distribution of the spacetime functions: fA(r,t) = radial scale function (dr²-factor) as shown in Figure 1, fB(r,t) = time scale function (dt²-factor) as shown in Figure 2, fY(r,t) = angular scale function (dθ² - dφ² factor) as shown in Figure 3, ρ(r,t) is the density as shown in Figure 4, u(r,t) is the velocity as shown in Figure 5, r and t are in Schwarzschild-units.

9.5. Ball-to-Ball Gravitational Collapse

We present the numerical results for the ball-to-ball gravitational collapse. This models the (fully GR) transition of an initially non-moving gas ball with radius r₀= R₁ to a smaller ball with radius R₂ = R₁/2 at the time t = T₁= 50, transition characteristic time T_cr= 0.2T₁.

For this purpose, we choose a corresponding ball-to-ball seed-function for ρ, Schwarzschild-expressions for fAfB fY, Newtonian form for the gravitational

acceleration axt(t,r) and calculating the velocity u(t,r) via the relativistic addition-of-velocity-theorem $u+\text{d}u=\frac{u+axt\ast \text{d}t}{1+u\ast axt\ast \text{d}t}$, which gives the formula $u\left(t,r\right)=\frac{\exp \left(2Iaxt\left(t,r\right)\right)-1}{\exp \left(2Iaxt\left(t,r\right)\right)+1}$, where $Iaxt\left(t,r\right)={\displaystyle {\int }_0^{t}axt\left(t_x,r\right)\text{d}{t}_x}$ is the integral

Figure 1. Radial scale function fA(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 2. Time scale function fB(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 3. Angular scale function fY(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

of the acceleration axt(t,r). This expression for u(t,r) fulfills automatically the physical feasibility condition u <1 (for dimensionless equations c = 1).

In a second step, we correct the boundary values (and the seed-functions) by solving the GR-equations on the boundary via Ritz-Galerkin-minimization as shown in Table 4.

Figure 4. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 5. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Table 4. Ritz-Galerkin global minimization.

The best results with error = 0.012 give the finite-difference minimization (Table 5), but the other methods yield very similar results (Table 6, Table 7).

The space-time function Y(r,t) is in normal Schwarzschild space-time Y(r,t)= r, here it reaches at the edge Y(r = 10,t)= 2.5, so area is scaled down by a factor of 4²= 16.

The density builds-up at the center to a dense inner ball with approximate radius r = 3 at the time t = 10= T_cr, and stays then approximately constant.

The radial velocity builds-up at t = 10= T_cr to about 0.2c, then rises at t = 30 and decreases sharply at final time.

Plot of numerical results:

time-radius distribution of the spacetime functions: fA(r,t) = radial scale function (dr²-factor) as shown in Figure 6, fB(r,t) = time scale function (dt²-factor) as shown in Figure 7, fY(r,t) = angular scale function (dθ² - dφ² factor) as shown in Figure 8, ρ(r,t) is the density as shown in Figure 9, u(r,t) is the velocity as shown in Figure 10, r and t are in Schwarzschild-units.

Plot of numerical results:

Time-radius distribution of the spacetime functions: fA(r,t) = radial scale function (dr²-factor) as shown in Figure 11, fB(r,t) = time scale function (dt²-factor)

Table 5. Finite-difference-minimization.

Figure 6. Radial scale function fA(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 7. Time scale function fB(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 8. Angular scale function fY(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

as shown in Figure 12, fY(r,t) = angular scale function (dθ² − dφ² factor) as shown in Figure 13, ρ(r,t) is the density as shown in Figure 14, u(r,t) is the velocity as shown in Figure 15, r and t are in Schwarzschild-units.

Plot of numerical results:

Time-radius distribution of the spacetime functions: fA(r,t) = radial scale

Figure 9. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 10. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s)

Table 6. R-profile-wavefront.

Figure 11. Radial scale function fA(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 12. Time scale function fB(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Table 7. Finite-difference-wavefront.

Figure 13. Angular scale function fY(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 14. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

function (dr²-factor) as shown in Figure 16, fB(r,t) = time scale function (dt²-factor) as shown in Figure 17, fY(r,t) = angular scale function (dθ² - dφ² factor) as shown in Figure 18, ρ(r,t) is the density as shown in Figure 19, u(r,t) is the velocity as shown in Figure 20, r and t are in Schwarzschild-units.

Figure 15. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 16. Radial scale function fA(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

9.6. Newtonian Gravitational Collapse

We present the numerical results for the ball-to-ball gravitational collapse in the Newtonian limit. This models the transition of an initially non-moving gas ball with radius r₀= R₁ to a smaller ball with radius R₂ = R₁/2 at the time t = T₁= 50, transition characteristic time T_cr= 0.2T₁, using simple Newtonian gravity and

Figure 17. Time scale function fB(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

Figure 18. Angular scale function fY(r,t)t = 0...50(r_s/c),r = 0…10(r_s).

the ideal-gas as the equation-of-state.

Here, there are only 2 differential equations derived from the energy conservation, for the 2 variable functions u(r,t), ρ(r,t).

Here, the density behaves essentially like in the full-GR calculation, the velocity rises to about 0.2c in the medium radius range, and stays low in the central

Figure 19. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 20. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Table 8. Ritz-Galerkin global minimization.

“bulb” area, where most of the mass is concentrated.

The Newtonian collapse was calculated with Ritz-Galerkin minimization (Table 8), finite-difference minimization (Table 9), and R-profile wavefront (Table 10).

Plot of numerical results:

Time-radius distribution ρ(r,t) the density as shown in Figure 21, u(r,t) the velocity as shown in Figure 22, r and t are in Schwarzschild-units.

Plot of numerical results:

Figure 21. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 22. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Time-radius distribution ρ(r,t) the density as shown in Figure 23,u(r,t) the velocity as shown in Figure 24, r and t are in Schwarzschild-units.

Table 9. Finite-difference-minimization.

Figure 23. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 24. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Table 10. R-profile-wavefront.

Figure 25. Density ρ(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Figure 26. Velocity u(r,t) t = 0...50(r_s/c),r = 0…10(r_s).

Plot of numerical results:

Time-radius distribution ρ(r,t) the density as shown in Figure 25, u(r,t) the velocity as shown in Figure 26, r and t are in Schwarzschild-units.

Conflicts of Interest

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

References

[1]	Mitra, A. (2000) Non-Occurrence of Trapped Surfaces and Black Holes in Spherical Gravitational Collapse. Foundations of Physics Letters, 13, 543-579, arXiv: gr-qc/9810038. https://doi.org/10.1023/A:1007810414531
[2]	Fliessbach, T. (1990) Allgemeine Relativitatstheorie. Bibliographisches Institut, Leipzig.
[3]	Kotake, K. (2013) Multiple Elements to Determine the Gravitational Wave Signatures of Core-Collapse Supernovae. Comptes Rendus Physique, 14, 318-315. arXiv: 1110.5107. https://doi.org/10.1016/j.crhy.2013.01.008
[4]	Naidu, N.F. (2008) Dynamics of Dissipative Gravitational Collapse. M.S. Thesis, University of Durban, Durban.
[5]	Joshi, P.S. and Malafarina, D. (2011) Recent Developments in Gravitational Collapse. International Journal of Modern Physics D, 20, 2641-2729. arXiv: 1201.3660. https://doi.org/10.1142/S0218271811020792
[6]	Lasky, P.D. (2007) A Unified Treatment of Gravitational Collapse in General Relativity. Ph.D. Thesis, Monash University, Melbourne.
[7]	Girichidis, Ph. (2012) Importance of the Initial Conditions for Star Formation. Ph.D. Thesis, University of Heidelberg, Heidelberg.
[8]	Herrera, L. and Barreto, W. (2011) Relativistic Gravitational Collapse in Comoving Coordinates. International Journal of Modern Physics D, 20, 1265-1288, arXiv: 1010.4087. https://doi.org/10.1142/S0218271811019426
[9]	Herrera, L., Di Prisco, A., Fuenmayor, E. and Troconis, O. (2009) Dynamics of Viscous Dissipative Gravitational Collapse. International Journal of Modern Physics D, 18, 129-145. arXiv: 0804.3584. https://doi.org/10.1142/S0218271809014285
[10]	Goswami, R. (2007) Gravitational Collapse of Dustlike Matter with Heat Flux. arXiv: 0707.1122.
[11]	Sharma, R. and Das, S. (2013) Collapse of a Relativistic Self-Gravitating Star with Radial Heat Flux. Journal of Gravity, 2013, Article ID: 659605. arXiv: 1304.47765. https://doi.org/10.1155/2013/659605
[12]	Nakazato, K. (2008) Gravitational Collapse and Neutrino Emission of Massive Stars. Ph.D. Thesis, Waseda University, Tokyo.
[13]	Sanwe, S., Saraykar, R.V. and Joshi, P.S. (2015) Gravitational Collapse with Equation of State. arXiv: 1207.3200.
[14]	Joshi, P.S., Malafarina, D. and Narayan, R. (2011) Equilibrium Configurations from Gravitational Collapse. Classical and Quantum Gravity, 28, Article ID: 235018. https://doi.org/10.1088/0264-9381/28/23/235018
[15]	Marshall, T.W. (2009) The Gravitational Collapse of a Dust Ball. arXiv: 0907.2339.
[16]	Hajicek, P. and Kiefer, C. (2001) Embedding Variables in the Canonical Theory of Gravitating Shells. Nuclear Physics B, 603, 531-554. arXiv: hep-th/0007004. https://doi.org/10.1016/S0550-3213(01)00141-9
[17]	Corda, C., Feleppa, F. and Tamburini, F. (2021) Quantum Oscillations in the Black Hole Horizon. arXiv: gr-qc/2104.05451.
[18]	Helm, J. (2019) New Solutions of the Tolman-Oppenheimer-Volkov Equation and of Kerr Space-Time with Matter. arXiv: 1434.0037. http://www.researchgate.net/publication/330737764
[19]	Helm, J. (2022) Gravitational Collapse. http://www.janhelm-works.com https://www.researchgate.net/publication/358270663