The Post–Quasi–Static Approximation: An Analytical Approach to Gravitational Collapse

1. Introduction

In the study of self–gravitating systems there are three possible regimes of evolution. The simplest one is the static (stationary when rotations are allowed) regime, which is characterized by the existence of a time–like Killing vector forming a vorticity–free congruence (in the stationary case the congruence is not vorticity–free ). In the coordinate system adapted to this congruence the metric and the physical variables are invariant with respect to translations along the time axe.

Next, we have the quasi–static regime (QSR), which applies when the system is assumed to evolve, but slowly enough, so that it can be considered to be in equilibrium at each moment (the TOV equation is satisfied at all times). This means that the sphere changes slowly, on a time scale that is very long as compared to the typical time in which the sphere reacts to a slight perturbation of hydrostatic equilibrium, this typical time scale is called hydrostatic time scale [5,6] (sometimes this time scale is also referred to as dynamical time scale, e.g. [7]). Thus, in this regime the system is always very close to hydrostatic equilibrium and its evolution may be regarded as a sequence of static models, where the time between any two states of equilibrium is neglected (see [8–10] for applications).

It is worth mentioning that the assumption above is very sensible because the hydrostatic time scale is very small for many phases of the life of the star [6]. It is of the order of 27 minutes for the Sun, $4.5$ seconds for a white dwarf and ${10}^{-4}$ seconds for a neutron star of one solar mass and 10 Km radius. Any of the stellar configurations mentioned above, generally (but not always), changes on a time scale that is very long compared to their respective hydrostatic time scales.

Finally, we have the dynamic regime where the system is out of equilibrium, meaning that the TOV equation is not satisfied. The system changes on a time scale which is smaller than the hydrostatic time scale.

At this point the question arises: can we approach the non–equilibrium by means of successive approximations? Or, equivalently: Is there life between quasi–equilibrium and non–equilibrium?

As it has been proved in the past (see [1–4] and references therein) the answer to the above questions is affirmative (is some cases at least), the corresponding regime is called post–quasi–static (PQSR), and can be regarded as the closest, non–equilibrium, regime to QSR. Before proceeding farther, some important remarks are in order

• First of all it should be stressed that the main motivation to consider the PQSR is to have the possibility to study, in the simplest possible way, those aspects of the object directly related to the non–equilibrium situation, which for obvious reasons cannot be described within the QSR.
• Since we are assuming the fact that we can approach the non–equilibrium by means of successive approximations, it goes without saying that not any self–gravitating fluid will satisfy this requirement. In particular is meaningless, from the physical point of view, to consider geodesic fluids in PQSR, since these fluids are always in the full dynamic regime (the only interaction in this case being the gravitational one).
• It also should be clear that unlike the two precedent regimes, there is not a unique definition for PQSR. Here we shall assume the definition proposed in [1–4]

Let us now elaborate on the main motivation of our endeavor with this work.

One of the most outstanding problems in relativistic astrophysics and gravitation theory today is to provide an accurate description of the gravitational collapse of a supermassive star. The final fate of such process (naked singularities, black holes, anything else), the mechanism behind a type II supernova event [11–17] or the structure and evolution of the compact object resulting from such a process [18–20], stand among the most interesting questions associated to that problem.

We have available three approaches to describe the gravitational collapse in the context of general relativity. On the one hand, one may resort to numerical methods [21–24], which allow for considering more realistic equations of state. However, the obtained results, in general, are restricted and highly model dependents. Also, specific difficulties, associated to numerical solutions of partial differential equations in presence of shocks may complicate further the problem.

On the other hand, one may use analytical solutions to Einstein equations, which are more suitable for a general discussion, and may be very useful in the study of the structure and evolution of self–gravitating systems, since they may be relatively simple to analyze but still may contain some of the essential features of a realistic situation (see for example [25–35] an references therein). However, often they are found, either for too simplistic equations of state and/or under additional heuristic assumptions whose justification is usually uncertain.

Between the two aforementioned approaches, we have semi-numerical techniques, which may be regarded as a “compromise” between the analytical and numerical approaches. These techniques are based on the PQSR approximation mentioned above, and were developed in [1–4] (see also [36,37]).

This third approach, starting from any interior (analytical) static spherically symmetric (“seed”) solution to Einstein equations, leads to a system of ordinary differential equations (referred to as surface equations) for quantities evaluated at the boundary surface of the fluid distribution, whose solution (numerical), allows for modeling, dynamic, self-gravitating spheres, whose static limit (whenever it exists) is the original “seed” solution.

The approach is based on the introduction of a set of conveniently defined “effective” variables (effective pressure and energy density) and an heuristic ansatz on the later, whose rationale and justification become intelligible within the context of the PQSR.

So far, the above mentioned approach has been used, by solving numerically the surface equations. In this work we complement the approach with a sensible physical condition, allowing us to avoid numerical integration, resorting exclusively to analytical methods. Such a condition appears to be the vanishing of the complexity factor, as defined in [38,39]. Other plausible conditions such as the homologous [39] and the quasi–homologous [40] conditions have been considered but were dismissed due to the facts that they, within the PQSR, lead to geodesic fluids.

Besides the vanishing complexity factor condition, we have to resort to additional sources of information in order to obtain a full description of the collapsing system. The number of possible strategies for doing that is very large. Here we emphasize, on the one hand, on conditions suggested by observables such as the luminosity profile and the gravitational redshift. On the other hand we propose some heuristic mathematical constrains, justified by previous experience on finding time–dependent solutions to Einstein equations, or, simply, by the fact that they allow a simple analytical integration.

The organization of the manuscript is as follows. In the next section we introduce the basic variables and definitions, as well as the Einstein and the transport equations. In Section III, we detail the junction conditions with the exterior spacetime, which is Vaidya. The complexity factor and the homologous and quasi–homologous evolution are defined in Section IV. A review of the approach is outlined in Section V, and some examples are analyzed in Section VI. Finally we include a discussion of the results and some concluding remarks in the last section.

2. Basic variables and equations

2.1. The metric

We consider a spherically symmetric distribution of collapsing fluid, bounded by a spherical surface $\Sigma$. The fluid is assumed to be locally anisotropic (principal stresses unequal) and undergoing dissipation in the form of heat flow (to model dissipation in the diffusion approximation). Physical arguments to consider such fluid distributions in the study of gravitational collapse may be found in [41–44] and references therein.

Using comoving coordinates, we write the line element in the form

$$d{s}^2=-A^{2}d{t}^2+B^{2}d{r}^2+R^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2),$$

(1)

where A, B and R are functions of t and r and are assumed positive. We number the coordinates $x^0=t$, $x^1=r$, $x^2=\theta$ and $x^3=\varphi$.

2.2. Energy–momentum tensor

The matter energy–momentum tensor $T_{\alpha \beta }$ inside $\Sigma$ has the form

$$\begin{array}{ccc}\hfill T_{\alpha \beta }& =& (\mu +P_{\perp })V_{\alpha }{V}_{\beta }+P_{\perp }{g}_{\alpha \beta }+(P_r-P_{\perp })K_{\alpha }{K}_{\beta }\hfill \\ & +& q_{\alpha }{V}_{\beta }+V_{\alpha }{q}_{\beta },\hfill \end{array}$$

(2)

where $\mu$ is the energy density, $P_r$ the radial pressure, $P_{\perp }$ the tangential pressure, $q^{\alpha }$ the heat flux, $V^{\alpha }$ the four-velocity of the fluid, and $K^{\alpha }$ a unit four-vector along the radial direction. These quantities satisfy

$$\begin{array}{ccc}\hfill V^{\alpha }{V}_{\alpha }& =& -1,V^{\alpha }{q}_{\alpha }=0,K^{\alpha }{K}_{\alpha }=1,\hfill \\ \hfill K^{\alpha }{V}_{\alpha }& =& 0.\hfill \end{array}$$

Since we assume the metric (1) comoving then

$$\begin{array}{c}\hfill V^{\alpha }=A^{-1}{\delta }_0^{\alpha },q^{\alpha }=q{B}^{-1}{\delta }_1^{\alpha },K^{\alpha }=B^{-1}{\delta }_1^{\alpha },\end{array}$$

(3)

where q is a function of t and r.

2.3. Kinematical variables

The four–acceleration $a_{\alpha }$ and the expansion $\Theta$ of the fluid are given by

$$a_{\alpha }=V_{\alpha ;\beta }{V}^{\beta },\Theta ={V^{\alpha }}_{;\alpha },$$

(4)

and its shear ${\sigma }_{\alpha \beta }$ by

$${\sigma }_{\alpha \beta }=V_{(\alpha ;\beta )}+a_{(\alpha }{V}_{\beta )}-\frac{1}{3}\Theta h_{\alpha \beta },$$

(5)

where $h_{\alpha \beta }=g_{\alpha \beta }+V_{\alpha }{V}_{\beta }.$

We do not explicitly add bulk viscosity to the system because it can be absorbed into the radial and tangential pressures, $P_r$ and $P_{\perp }$, of the collapsing fluid.

From (4) with (3) we have for the four–acceleration and its scalar a,

$$a_1=\frac{A^{\prime }}{A},a^2=a^{\alpha }{a}_{\alpha }={\left(\frac{A^{\prime }}{AB}\right)}^2,$$

(6)

where $a^{\alpha }=a{K}^{\alpha }$, and for the expansion

$$\Theta =\frac{1}{A}\left(\frac{\dot{B}}{B}+2\frac{\dot{R}}{R}\right),$$

(7)

where the prime stands for r differentiation and the dot stands for differentiation with respect to t. With (3) we obtain for the shear (5) its non zero components

$${\sigma }_{11}=\frac{2}{3}{B}^2\sigma ,{\sigma }_{22}=\frac{{\sigma }_{33}}{{sin}^2\theta }=-\frac{1}{3}{R}^2\sigma ,$$

(8)

and its scalar

$${\sigma }^{\alpha \beta }{\sigma }_{\alpha \beta }=\frac{2}{3}{\sigma }^2,$$

(9)

where

$$\sigma =\frac{1}{A}\left(\frac{\dot{B}}{B}-\frac{\dot{R}}{R}\right).$$

(10)

Then, the shear tensor can be written as

$${\sigma }_{\alpha \beta }=\sigma \left(K_{\alpha }{K}_{\beta }-\frac{1}{3}{h}_{\alpha \beta }\right).$$

(11)

2.4. Transport equations

In the dissipative case we shall need a transport equation in order to find the temperature distribution and its evolution. Assuming a causal dissipative theory (e.g. the Israel– Stewart theory [45–47]) the transport equation for the heat flux reads

$$\begin{array}{ccc}\hfill \tau h^{\alpha \beta }{V}^{\gamma }{q}_{\beta ;\gamma }& +& q^{\alpha }=-k{h}^{\alpha \beta }\left(T_{,\beta }+T{a}_{\beta }\right)\hfill \\ & -& \frac{1}{2}k{T}^2{\left(\frac{\tau V^{\beta }}{\kappa T^2}\right)}_{;\beta }{q}^{\alpha },\hfill \end{array}$$

(12)

where k, T and $\tau$ denote thermal conductivity, temperature and relaxation time respectively.

In the spherically symmetric case under consideration, the transport equation has only one independent component which may be obtained from (12) by contracting with the unit spacelike vector $K^{\alpha }$, it reads

$$\tau V^{\alpha }{q}_{,\alpha }+q=-k\left(K^{\alpha }{T}_{,\alpha }+Ta\right)-\frac{1}{2}k{T}^2{\left(\frac{\tau V^{\alpha }}{\kappa T^2}\right)}_{;\alpha }q.$$

(13)

2.5. Field equations

The Einstein field equations for the interior spacetime (1) can be written as

$$\begin{array}{c}\hfill 8\pi \mu A^2=\left(2\frac{\dot{B}}{B}+\frac{\dot{R}}{R}\right)\frac{\dot{R}}{R}-{\left(\frac{A}{B}\right)}^2\left[2\frac{R^{″}}{R}+{\left(\frac{R^{\prime }}{R}\right)}^2-2\frac{B^{\prime }}{B}\frac{R^{\prime }}{R}-{\left(\frac{B}{R}\right)}^2\right],\end{array}$$

(14)

$$4\pi qAB=\left(\frac{{\dot{R}}^{\prime }}{R}-\frac{\dot{B}}{B}\frac{R^{\prime }}{R}-\frac{\dot{R}}{R}\frac{A^{\prime }}{A}\right),$$

(15)

$$\begin{array}{c}\hfill 8\pi P_r{B}^2=-{\left(\frac{B}{A}\right)}^2\left[2\frac{\ddot{R}}{R}-\left(2\frac{\dot{A}}{A}-\frac{\dot{R}}{R}\right)\frac{\dot{R}}{R}\right]+\left(2\frac{A^{\prime }}{A}+\frac{R^{\prime }}{R}\right)\frac{R^{\prime }}{R}-{\left(\frac{B}{R}\right)}^2,\end{array}$$

(16)

$$\begin{array}{c}\hfill 8\pi P_{\perp }{R}^2=-{\left(\frac{R}{A}\right)}^2\left[\frac{\ddot{B}}{B}+\frac{\ddot{R}}{R}-\frac{\dot{A}}{A}\left(\frac{\dot{B}}{B}+\frac{\dot{R}}{R}\right)+\frac{\dot{B}}{B}\frac{\dot{R}}{R}\right]\\ \hfill +{\left(\frac{R}{B}\right)}^2\left[\frac{A^{″}}{A}+\frac{R^{″}}{R}-\frac{A^{\prime }}{A}\frac{B^{\prime }}{B}+\left(\frac{A^{\prime }}{A}-\frac{B^{\prime }}{B}\right)\frac{R^{\prime }}{R}\right].\end{array}$$

(17)

Observe that if functions $A(t,r)$, $B(t,r)$ and $R(t,r)$ are completely determined, the system above becomes an algebraic system of four equations for the four unknown functions $\mu$, q, $P_r$, and $P_{\perp }$ which can be obtained without further information.

2.6. Mass and areal velocity

Following Misner and Sharp [48], let us now introduce the mass function $m(t,r)$ (see also [49]), defined by

$$m=\frac{R^3}{2}{R_{23}}^{23}=\frac{R}{2}\left[{\left(\frac{\dot{R}}{A}\right)}^2-{\left(\frac{R^{\prime }}{B}\right)}^2+1\right].$$

(18)

It is useful to introduce the proper time derivative $D_T$ given by

$$D_T=\frac{1}{A}\frac{\partial }{\partial t},$$

(19)

and the proper radial derivative $D_R$,

$$D_R=\frac{1}{R^{\prime }}\frac{\partial }{\partial r},$$

(20)

where R defines the areal radius of a spherical surface inside $\Sigma$ (as measured from its area).

Using (19) we can define the velocity U of the collapsing fluid as the variation of the areal radius with respect to proper time, i.e.

$$U=D_{T}R.$$

(21)

Then (18) can be rewritten as

$$E\equiv \frac{R^{\prime }}{B}={\left(1+U^2-\frac{2m}{R}\right)}^{1/2}.$$

(22)

Using (14)-(16) with (19) and (20) we obtain from (18)

$$\begin{array}{c}\hfill D_{T}m=-4\pi \left(P_{r}U+qE\right)R^2,\end{array}$$

(23)

and

$$\begin{array}{c}\hfill D_{R}m=4\pi \left(\mu +q\frac{U}{E}\right)R^2.\end{array}$$

(24)

Next, the three–acceleration $D_{T}U$ of an in-falling particle inside $\Sigma$ can be obtained by using (16), (18) and (22), producing

$$D_{T}U=-\frac{m}{R^2}-4\pi P_{r}R+E\frac{A^{\prime }}{AB},$$

(25)

or

$$\frac{A^{\prime }}{A}=\frac{4\pi RB}{E}\left(\frac{D_{T}U}{4\pi R}+\frac{m}{4\pi R^3}+P_r\right).$$

(26)

Finally, from the Bianchi identities we obtain

$$\begin{array}{ccc}\hfill (\mu +P_r)D_{T}U& =& -(\mu +P_r)\left(\frac{m}{R^2}+4\pi P_{r}R\right)-E^2\left[D_R{P}_r+\frac{2}{R}(P_r-P_{\perp })\right]\hfill \\ & -& E\left[D_{T}q+2q\left(\frac{2U}{R}+\sigma \right)\right].\hfill \end{array}$$

(27)

The physical meaning of different terms in (27) has been discussed in detail in [43,44]. Suffice is to say in this point that the first term on the right hand side describes the gravitational force term.

3. The exterior spacetime and junction conditions

Outside $\Sigma$ we assume we have the Vaidya spacetime (i.e. we assume all outgoing radiation is massless), described by

$$d{s}^2=-\left[1-\frac{2M\left(v\right)}{\rho }\right]d{v}^2-2d\rho dv+{\rho }^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2),$$

(28)

where $M\left(v\right)$ denotes the total mass, and v is the retarded time.

The matching of the full non–adiabatic sphere (including viscosity) to the Vaidya spacetime, on the surface $r=r_{\Sigma }=$ constant, was discussed in [50].

Now, from the continuity of the first differential form it follows (see [50] for details),

$$Adt\stackrel{\Sigma }{=}dv\left(1-\frac{2M\left(v\right)}{\rho }\right)\stackrel{\Sigma }{=}d\tau ,$$

(29)

$$R\stackrel{\Sigma }{=}\rho \left(v\right),$$

(30)

and

$${\left(\frac{dv}{d\tau }\right)}^{-2}\stackrel{\Sigma }{=}\left(1-\frac{2m}{\rho }+2\frac{d\rho }{dv}\right),$$

(31)

where $\tau$ denotes the proper time measured on $\Sigma$.

The continuity of the second differential form produces

$$m(t,r)\stackrel{\Sigma }{=}M\left(v\right),$$

(32)

and

$$\begin{array}{c}\hfill 2\left(\frac{{\dot{R}}^{\prime }}{R}-\frac{\dot{B}}{B}\frac{R^{\prime }}{R}-\frac{\dot{R}}{R}\frac{A^{\prime }}{A}\right)\stackrel{\Sigma }{=}-\frac{B}{A}\left[2\frac{\ddot{R}}{R}-\left(2\frac{\dot{A}}{A}-\frac{\dot{R}}{R}\right)\frac{\dot{R}}{R}\right]+\frac{A}{B}\left[\left(2\frac{A^{\prime }}{A}+\frac{R^{\prime }}{R}\right)\frac{R^{\prime }}{R}-{\left(\frac{B}{R}\right)}^2\right],\end{array}$$

(33)

where $\stackrel{\Sigma }{=}$ means that both sides of the equation are evaluated on $\Sigma$ (observe a misprint in eq.(40) in [50] and a slight difference in notation).

Comparing (33) with (15) and (16) one obtains

$$q\stackrel{\Sigma }{=}{P}_r.$$

(34)

Thus the matching of (1) and (28) on $\Sigma$ implies (32) and (34).

Also, we have

$$q\stackrel{\Sigma }{=}\frac{L}{4\pi {\rho }^2},$$

(35)

where $L_{\Sigma }$ denotes the total luminosity of the sphere as measured on its surface and is given by

$$L\stackrel{\Sigma }{=}{L}_{\infty }{\left(1-\frac{2m}{\rho }+2\frac{d\rho }{dv}\right)}^{-1},$$

(36)

and where

$$L_{\infty }=-\frac{dM}{dv}\stackrel{\Sigma }{=}-\left[\frac{dm}{dt}\frac{dt}{d\tau }{\left(\frac{dv}{d\tau }\right)}^{-1}\right],$$

(37)

is the total luminosity measured by an observer at rest at infinity.

The boundary redshift $z_{\Sigma }$ is given by

$$\frac{dv}{d\tau }\stackrel{\Sigma }{=}1+z,$$

(38)

with

$$\frac{dv}{d\tau }\stackrel{\Sigma }{=}{\left(\frac{R^{\prime }}{B}+\frac{\dot{R}}{A}\right)}^{-1}.$$

(39)

Therefore the time of formation of the black hole is given by

$$\left(\frac{R^{\prime }}{B}+\frac{\dot{R}}{A}\right)\stackrel{\Sigma }{=}E+U\stackrel{\Sigma }{=}0.$$

(40)

Also observe than from (31), (36) and (39) it follows

$$L\stackrel{\Sigma }{=}\frac{L_{\infty }}{{(E+U)}^2},$$

(41)

and from (21), (22), (31) and (39)

$$\frac{d\rho }{dv}\stackrel{\Sigma }{=}U(U+E).$$

(42)

4. The complexity factor

The condition we shall impose on our system in order to integrate analytically the ensuing differential equations, is the vanishing of the complexity factor. This is a scalar function intended to measure the degree of complexity of a given fluid distribution [38,39], and is related to the so called structure scalars [51].

As shown in [38,39] the complexity factor is identified with the scalar function $Y_{TF}$ which defines the trace–free part of the electric Riemann tensor (see [51] for details).

Thus, let us define tensor $Y_{\alpha \beta }$ by

$$Y_{\alpha \beta }=R_{\alpha \gamma \beta \delta }{V}^{\gamma }{V}^{\delta },$$

(43)

which may be expressed in terms of two scalar functions $Y_T,Y_{TF}$, as

$$\begin{array}{c}\hfill Y_{\alpha \beta }=\frac{1}{3}{Y}_T{h}_{\alpha \beta }+Y_{TF}\left(K_{\alpha }{K}_{\beta }-\frac{1}{3}{h}_{\alpha \beta }\right).\end{array}$$

(44)

Then after lengthy but simple calculations, using field equations, we obtain (see [39,40] for details)

$$Y_{TF}=-8\pi \Pi +\frac{4\pi }{R^3}{\int }_0^r{R}^3\left(D_R\mu -3q\frac{U}{RE}\right)R^{\prime }d\tilde{r}.$$

(45)

In terms of the metric functions the scalar $Y_{TF}$ reads

$$\begin{array}{c}\hfill Y_{TF}=\frac{1}{A^2}\left[\frac{\ddot{R}}{R}-\frac{\ddot{B}}{B}+\frac{\dot{A}}{A}\left(\frac{\dot{B}}{B}-\frac{\dot{R}}{R}\right)\right]+\frac{1}{B^2}\left[\frac{A^{″}}{A}-\frac{A^{\prime }}{A}\left(\frac{B^{\prime }}{B}+\frac{R^{\prime }}{R}\right)\right].\end{array}$$

(46)

4.1. The homologous and quasi–homologous evolution

Another set of possible conditions, which might be considered in order to avoid numerical integration, are conditions on the pattern of evolution.

One of these conditions is represented by the homologous evolution (H). In [39] it was assumed that the H evolution describes the simplest mode of evolution of the fluid distribution. Such a condition is defined by

$$U=\tilde{a}\left(t\right)R,\hspace{1em}\hspace{1em}\tilde{a}\equiv \frac{U_{\Sigma }}{R_{\Sigma }},$$

(47)

and

$$\frac{R_I}{R_{II}}=\mathrm{constant},$$

(48)

where $R_I$ and $R_{II}$ denote the areal radii of two concentric shells ($I,II$) described by $r=r_I=\mathrm{constant}$, and $r=r_{II}=\mathrm{constant}$, respectively.

These relationships are reminiscent of the homologous evolution in Newtonian hydrodynamics [5–7].

The important point that we want to stress here is that, in the relativistic regime, (47) does not imply (48).

Indeed, (47) implies that for two comoving shells of fluids $I,II$ we have

$$\frac{U_I}{U_{II}}=\frac{A_{II}{\dot{R}}_I}{A_I{\dot{R}}_{II}}=\frac{R_I}{R_{II}},$$

(49)

which implies (48) only if the fluid is geodesic ($A=\mathrm{constant}$). However, in the non–relativistic regime, (48) always follows from the condition that the radial velocity is proportional to the radial distance.

Another possible condition (less restrictive) could be represented by the so called “quasi–homologous” regime ($QH$), characterized by condition (47) alone, which implies (see [40] for details)

$$\frac{4\pi }{R^{\prime }}Bq+\frac{\sigma }{R}=0.$$

(50)

Thus the H condition implies (48) and (50), while the $QH$ condition only requires (50).

However both conditions lead (within the PQSR) to geodesics fluids, which, as already mentioned, are physically without interest.

Indeed, writing (15) as

$$4\pi qB=\frac{1}{3}{(\Theta -\sigma )}^{\prime }-\sigma \frac{R^{\prime }}{R},$$

(51)

and combining with condition (50), we obtain

$${(\Theta -\sigma )}^{\prime }=0,$$

(52)

whereas, using (7) and (10) we get

$${\left(\Theta -\sigma \right)}^{\prime }={\left(\frac{3}{A}\frac{\dot{R}}{R}\right)}^{\prime }=0.$$

(53)

But in the PQSR we have (see equation (65) in Section 5.3 below) $R=\kappa \left(t\right)r$ where $\kappa$ is an arbitrary function of t, producing at once that

$$A^{\prime }=0,$$

(54)

implying that the fluid is geodesic, as it follows from (6).

Thus from physical considerations we must exclude the H or the $QH$ conditions for the mode of evolution.

We shall next, define mathematically the three regimes of evolution mentioned in the Introduction, in order to understand the rationale behind the proposed approach.

5. Evolution regimes

Let us now express the three possible regimes of evolution, in terms of the metric and physical variables.

5.1. Static regime

In this case all time derivatives vanish, implying:

$$q=U=\Theta =\sigma =0.$$

(55)

Since $B=B\left(r\right);A=A\left(r\right);R=R\left(r\right)$, reparametrizing r, we may write the line element in the form:

$$d{s}^2=-A^{2}d{t}^2+B^{2}d{r}^2+r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2).$$

(56)

Thus, the “Euler” equation (27) becomes the well known TOV equation of hydrostatic equilibrium for an anisotropic fluid

$$\begin{array}{c}\hfill P_r^{\prime }+\frac{2}{r}(P_r-P_{\perp })=-\frac{(\mu +P_r)}{r(r-2m)}(m+4\pi P_r{r}^3).\end{array}$$

(57)

The Einstein equations in this case read:

$$\begin{array}{c}\hfill 8\pi \mu A^2=-{\left(\frac{A}{B}\right)}^2\left[{\left(\frac{1}{r}\right)}^2-2\frac{B^{\prime }}{Br}-{\left(\frac{B}{r}\right)}^2\right],\end{array}$$

(58)

$$\begin{array}{c}\hfill 8\pi P_r{B}^2=\left(2\frac{A^{\prime }}{A}+\frac{1}{r}\right)\frac{1}{r}-{\left(\frac{B}{r}\right)}^2,\end{array}$$

(59)

$$\begin{array}{c}\hfill 8\pi P_{\perp }{r}^2={\left(\frac{r}{B}\right)}^2\left[\frac{A^{″}}{A}-\frac{A^{\prime }}{A}\frac{B^{\prime }}{B}+\left(\frac{A^{\prime }}{A}-\frac{B^{\prime }}{B}\right)\frac{1}{r}\right].\end{array}$$

(60)

Also, for the mass function we have

$$m=\frac{r}{2}\left(1-\frac{1}{B^2}\right)\Rightarrow B^2={\left(1-\frac{2m}{r}\right)}^{-1},$$

(61)

or

$$m=4\pi {\int }_0^r\mu r^{2}dr,$$

(62)

and for the metric function A, we have from (26)

$$ln\left(\frac{A}{A_{\Sigma }}\right)={\int }_{r_{\Sigma }}^r\frac{(m+4\pi r^3{P}_r)}{r(r-2m)}dr.$$

(63)

The important point to keep in mind is that if the radial dependence of $\mu$ and $P_r$ is known, the metric functions are determined from (61–63).

5.2. Quasi-static regime (QSR)

As mentioned before, in this regime the system is assumed to evolve, but sufficiently slow, so that it can be considered to be in equilibrium at each moment (Eq. (57) is satisfied).

Let us now translate this assumption in conditions to U, metric and kinematical functions.

The QSR implies that

• The areal velocity U as well as other kinematical variables are small, (of order $O\left(ϵ\right)$, with $\left|ϵ\right|<<1$) which in turn implies that dissipative variables and all first order time derivatives of metric functions are also small, implying that we shall neglect terms of order ${ϵ}^2$ and higher.
• From the above and the fact that the system always satisfies the equation of hydrostatic equilibrium, it follows from (27) that second time derivatives of metric functions can be neglected.

Thus in QSR we have

$$O\left(U^2\right)={\dot{A}}^2={\dot{B}}^2=\dot{A}\dot{B}=\ddot{R}=\ddot{B}\approx 0$$

(64)

and the radial dependence of the metric functions as well as that of physical variables is the same as in the static case. The only difference with the latter case being that these variables depend upon time according to equation (15).

5.3. Post–quasi–static regime (PQSR)

In the two regimes considered above the system is always in (or very close to) hydrostatic equilibrium. Let us now move one step forward into non–equilibrium and let us assume that (57) is not satisfied.

Then the question arises: What is the closest situation to QSR not satisfying eq. (57)? Such a situation is described by what we call PQSR.

Now, since in both, the static and QSR regimes, the radial dependence of metric variables is the same, we shall keep that radial dependence as much as possible, but of course the time dependence of those variables is such that now (64) is not satisfied.

Then from the above we write

$$R=r\kappa \left(t\right),$$

(65)

where $\kappa$ is an arbitrary (dimensionless) function of t, to be determined later.

Taking into account (22) and (65), we rewrite the metric as follows

$$d{s}^2=-A^{2}d{t}^2+{\kappa }^2[E^{-2}d{r}^2+r^2(d{\theta }^2+sin{\theta }^{2}d{\varphi }^2)].$$

(66)

Next, defining the effective mass as

$$m_{eff}\equiv m-\frac{1}{2}R{U}^2,$$

(67)

we obtain

$$E^2=1-\frac{2m_{eff}}{R}.$$

(68)

Then, equations (24) and (26) can be written as

$$\begin{array}{ccc}\hfill \frac{1}{\kappa }{m}_{eff}^{\prime }& =& 4\pi R^2{\mu }_{eff},\hfill \\ \hfill \end{array}$$

(69)

$$\begin{array}{ccc}\hfill \frac{1}{\kappa }{(lnA)}^{\prime }& =& \frac{4\pi R^2{P}_{eff}+m_{eff}/R}{R-2m_{eff}},\hfill \end{array}$$

(70)

with

$$\begin{array}{ccc}\hfill {\mu }_{eff}& =& \mu +\frac{qU}{E}-\frac{U{D}_{R}U}{4\pi R}-\frac{U^2}{8\pi R^2},\hfill \\ \hfill \end{array}$$

(71)

$$\begin{array}{ccc}\hfill P_{eff}& =& P_r+\frac{D_{T}U}{4\pi R}+\frac{U^2}{8\pi R^2},\hfill \end{array}$$

(72)

where we have followed the terminology used in [2–4] and call ${\mu }_{eff}$ and $P_{eff}$ the “effective density” and the “effective pressure”, respectively. The meaning of these variables will become clear in the discussion below, however we remark at this point that in the static and QSR cases, the effective variables coincide with the corresponding physical variables. (in what concerns their radial dependence).

Next, from (69)–(72), with (65) we may write

$$\frac{1}{{\kappa }^3}{m}_{eff}={\int }_0^{r}4\pi r^2{\mu }_{eff}dr,$$

(73)

$$\frac{1}{\kappa }ln\left(\frac{A}{A_{\Sigma }}\right)={\int }_{r_{\Sigma }}^r\left[\frac{4\pi R^3{P}_{eff}+m_{eff}}{R(R-2m_{eff})}\right]dr.$$

(74)

From the above, it follows at once that if $R=\kappa \left(t\right)r$ and ${\mu }_{eff}$ shares the same radial dependence as $\mu$ in the static case, then obviously the radial dependence of $m_{eff}$ will be the same as in the static case. The inverse is true of course, if the radial dependence of $m_{eff}$ is the same as in the static case, then ${\mu }_{eff}$ shares the same radial dependence as $\mu$ static.

On the other hand, if besides the assumption above, we assume that $P_{eff}$ shares the same radial dependence as $P_r$ static, then it follows from (74) that A shares the same radial dependence as in the static case.

All these considerations provided the rationale for the algorithm as exposed in [4]. Thus, the proposed method, starting from any interior (analytical) static spherically symmetric (“seed”) solution to Einstein equations, leads to a system of ordinary differential equations for quantities evaluated at the boundary surface of the fluid distribution, whose solution (numerical), allows for modeling, dynamic self-gravitating spheres, whose static limit is the original “seed” solution.

In this work, motivated by our interest in resorting to purely analytical methods we shall modify the algorithm described in [4].

Specifically, the main steps of the formalism we propose may be summarized as follows.

• Take an interior (“‘seed”) solution to Einstein equations, representing a fluid distribution of matter in equilibrium, with a given

  $${\mu }_{st}={\mu }_{st}\left(r\right);P_{rst}=P_{rst}\left(r\right).$$
• Assume that the r dependence of the effective density is the same as that of ${\mu }_{st}$, and $R=r\kappa \left(t\right)$.
• Impose the vanishing complexity factor condition.
• From the two conditions above we are able to determine the metric functions up to two arbitrary functions of t.
• For these functions of t one has the junction condition (33).
• In order to determine the remaining function and to integrate analytically (33) we have a large number of possible strategies. Here we shall mention some of them, which may be based on the information obtained from the observables of the collapsing star. Such observables are the luminosity and the redshift. Alternatively we may assume additional heuristic constraints on some other physical variables, or ad hoc mathematical conditions based in previous works on gravitational collapse, or simply justified by the fact that it allows a simple integration of (33). We list below some possible strategies of the kind mentioned above.
  • Assuming a specific luminosity profile obtained from observations and using (36) or (37) we obtain a relationship between the two arbitrary functions of t mentioned above, thereby reducing (38) to an ordinary differential equation for one variable.
  • Assuming a specific form for the evolution of the redshift we obtain again a relationship between the two arbitrary functions of t
  • We may consider a specific pattern evolution of the areal radius of the star, or equivalently of its velocity ($U_{\Sigma }$). This could be useful if for example we want to check the possibility of a bouncing of the boundary surface.
  • Assuming different profiles of either one of the two arbitrary functions of t, we can look for conditions allowing the formation (or not) of a horizon, according to (40).

6. Modeling

We shall now proceed to implement the approach for modeling that we propose, and illustrate it by means of two examples.

Let us first write the general expressions for the field equations and $Y_{TF}$. Using (14)–(17), (46) and (65), we obtain

$$8\pi \mu =\frac{1}{A^2}\left(\frac{2\dot{B}}{B}+\frac{\dot{\kappa }}{\kappa }\right)\frac{\dot{\kappa }}{\kappa }-\frac{1}{B^2}\left(\frac{1}{r}-\frac{2B^{\prime }}{B}\right)\frac{1}{r}+\frac{1}{r^2{\kappa }^2},$$

(75)

$$4\pi q=\frac{1}{AB}\left(\frac{\dot{\kappa }}{r\kappa }-\frac{\dot{B}}{rB}-\frac{A^{\prime }\dot{\kappa }}{A\kappa }\right),$$

(76)

$$\begin{array}{c}\hfill 8\pi P_r=-\frac{1}{A^2}\left[\frac{2\ddot{\kappa }}{\kappa }-\left(\frac{2\dot{A}}{A}-\frac{\dot{\kappa }}{\kappa }\right)\frac{\dot{\kappa }}{\kappa }\right]\\ \hfill +\frac{1}{B^2}\left(\frac{2A^{\prime }}{A}+\frac{1}{r}\right)\frac{1}{r}-\frac{1}{r^2{\kappa }^2},\end{array}$$

(77)

$$\begin{array}{c}\hfill 8\pi P_{\perp }=-\frac{1}{A^2}\left[\frac{\ddot{B}}{B}+\frac{\ddot{\kappa }}{\kappa }-\frac{\dot{A}}{A}\left(\frac{\dot{B}}{B}+\frac{\dot{\kappa }}{\kappa }\right)+\frac{\dot{B}\dot{\kappa }}{B\kappa }\right]\\ \hfill +\frac{1}{B^2}\left[\frac{A^{″}}{A}-\frac{A^{\prime }{B}^{\prime }}{AB}+\left(\frac{A^{\prime }}{A}-\frac{B^{\prime }}{B}\right)\frac{1}{r}\right],\end{array}$$

(78)

and

$$\begin{array}{c}\hfill Y_{TF}=\frac{1}{A^2}\left[\frac{\ddot{\kappa }}{\kappa }-\frac{\ddot{B}}{B}+\frac{\dot{A}}{A}\left(\frac{\dot{B}}{B}-\frac{\dot{\kappa }}{\kappa }\right)\right]\\ \hfill +\frac{1}{B^2}\left[\frac{A^{″}}{A}-\frac{A^{\prime }}{A}\left(\frac{B^{\prime }}{B}+\frac{1}{r}\right)\right].\end{array}$$

(79)

Let us first consider the $q=0$ case, which using (76) produces

$$\frac{1}{r}\left(\frac{\dot{\kappa }}{\kappa }-\frac{\dot{B}}{B}\right)-\frac{A^{\prime }\dot{\kappa }}{A\kappa }=0.$$

(80)

Since at $r=0$, A is different from zero, we must impose

$$\frac{\dot{\kappa }}{\kappa }=\frac{\dot{B}}{B},\Rightarrow \mathrm{B}\mathrm{separable},$$

(81)

and

$$\frac{A^{\prime }\dot{\kappa }}{A\kappa }=0,\Rightarrow A=A\left(t\right),\left(\mathrm{geodesic}\right).$$

(82)

Since the geodesic case in the PQSR should be dismissed by reasons exposed before, we shall consider exclusively dissipative systems.

Then since $q\ne 0$, it follows from (76) that B is separable

$$B(r,t)=\kappa \left(t\right)\beta \left(r\right),$$

(83)

here $\beta$ is an arbitrary dimensionless function of r, and

$$4\pi q=-\frac{1}{A\kappa \beta }\left(\frac{A^{\prime }\dot{\kappa }}{A\kappa }\right).$$

(84)

It is worth stressing that using (83) in (10) it follows at once that $\sigma =0$. Thus all our models will be shear–free.

Next, assuming $Y_{TF}=0$ we obtain from (79)

$$\frac{A^{″}}{A^{\prime }}=\frac{\beta {\left(r\right)}^{\prime }}{\beta \left(r\right)}+\frac{1}{r},$$

(85)

whose solution reads

$$A=\alpha \int \beta \left(r\right)rdr+f(t),$$

(86)

where f is arbitrary function of integration, and, by reparametrizying t, another function of integration has been put equal to $\alpha =constant=1$, with dimensions $[1/r^2]$.

Then eqs.(75)–(78) take the form

$$8\pi \mu =\frac{1}{A^2}\frac{3{\dot{\kappa }}^2}{{\kappa }^2}-\frac{1}{{\beta }^{2}r{\kappa }^2}\left(\frac{1}{r}-\frac{2{\beta }^{\prime }}{\beta }\right)+\frac{1}{r^2{\kappa }^2},$$

(87)

$$4\pi q=-\frac{\alpha r\dot{\kappa }}{A^2{\kappa }^2},$$

(88)

$$\begin{array}{c}\hfill 8\pi P_r=-\frac{1}{A^2}\left(\frac{2\ddot{\kappa }}{\kappa }-\frac{2\dot{f}\dot{\kappa }}{A\kappa }+\frac{{\dot{\kappa }}^2}{{\kappa }^2}\right)\\ \hfill +\frac{1}{{\beta }^{2}r{\kappa }^2}\left(\frac{2\alpha \beta r}{A}+\frac{1}{r}\right)-\frac{1}{r^2{\kappa }^2},\end{array}$$

(89)

$$\begin{array}{c}\hfill 8\pi P_{\perp }=-\frac{1}{A^2}\left(\frac{2\ddot{\kappa }}{\kappa }-\frac{2\dot{f}\dot{\kappa }}{A\kappa }+\frac{{\dot{\kappa }}^2}{{\kappa }^2}\right)\\ \hfill +\frac{1}{{\beta }^2{\kappa }^2}\left(\frac{2\alpha \beta }{A}-\frac{{\beta }^{\prime }}{r\beta }\right),\end{array}$$

(90)

where A is given by (86).

Also, from (89) and (90)

$$8\pi (P_r-P_{\perp })=\frac{1}{{\beta }^2{\kappa }^{2}r}\left(\frac{1}{r}+\frac{{\beta }^{\prime }}{\beta }\right)-\frac{1}{{\kappa }^2{r}^2}.$$

(91)

Using (65) and (83) we can write

$${\mu }_{eff}=\mu +\frac{qr\beta \dot{\kappa }}{A}-\frac{{\dot{\kappa }}^2}{8\pi A^2{\kappa }^2}\left(3-\frac{2\alpha r^2\beta }{A}\right),$$

(92)

$$P_{eff}=P_r+\frac{1}{4\pi A^2}\left(\frac{\ddot{\kappa }}{\kappa }-\frac{\dot{\kappa }\dot{f}}{\kappa A}\right)+\frac{{\dot{\kappa }}^2}{8\pi A^2{\kappa }^2},$$

(93)

where A is given by (86).

We shall now use the equations above to present some analytical models of collapsing objects. It should be stressed that the obtained models are presented with the sole purpose of illustrating the method, and not to describe any specific astrophysical scenario.

6.1. A model with homogenous effective energy–density

The first model, is obtained by taking as our “seed” solution the well known Schwarzschild interior solution characterized by homogeneous energy-density and isotropic pressure.

Thus, assuming ${\mu }_{eff}=F\left(t\right)$, where $F\left(t\right)$ is an arbitrary function with units $[1/r^2]$, we obtain from (73),

$$m_{eff}=\frac{4\pi r^3{\kappa }^{3}F\left(t\right)}{3},$$

(94)

and with (22) and (68) we have

$$\frac{1}{r^2}\left(1-\frac{1}{{\beta }^2}\right)=\frac{8\pi {\kappa }^{2}F\left(t\right)}{3},$$

(95)

then

$${\beta }^2=\frac{1}{1-c{r}^2},$$

(96)

where c is a constant, with the same units as $F\left(t\right)$, given by

$$c=\frac{8\pi {\kappa }^{2}F\left(t\right)}{3}.$$

(97)

With this we have for A

$$A=f\left(t\right)-\frac{\alpha }{c}\sqrt{1-c{r}^2},$$

(98)

and for the field equations

$$8\pi \mu =\frac{3c^2{\dot{\kappa }}^2}{{\left(cf-\alpha \sqrt{1-c{r}^2}\right)}^2{\kappa }^2}+\frac{3c}{{\kappa }^2},$$

(99)

$$4\pi q=-\frac{\alpha c^{2}r\dot{\kappa }}{{\left(cf-\alpha \sqrt{1-c{r}^2}\right)}^2{\kappa }^2},$$

(100)

$$\begin{array}{ccc}\hfill 8\pi P_r=8\pi P_{\perp }& =& -\frac{c^2}{{\left(cf-\alpha \sqrt{1-c{r}^2}\right)}^2}\left[\frac{2\ddot{\kappa }}{\kappa }-\frac{2c\dot{f}\dot{\kappa }}{\left(cf-\alpha \sqrt{1-c{r}^2}\right)\kappa }+\frac{{\dot{\kappa }}^2}{{\kappa }^2}\right]\hfill \\ & +& \frac{2c\alpha \sqrt{1-c{r}^2}}{\left(cf-\alpha \sqrt{1-c{r}^2}\right){\kappa }^2}-\frac{c}{{\kappa }^2}.\hfill \end{array}$$

(101)

On the surface $\Sigma$, from (33) or (34) we obtain

$$2\kappa \ddot{\kappa }-\frac{2c\dot{f}\kappa \dot{\kappa }}{\left(cf-\alpha \sqrt{1-c{r}^2}\right)}+{\dot{\kappa }}^2-2\alpha r\dot{\kappa }\stackrel{\Sigma }{=}4\alpha f\sqrt{1-c{r}^2}-c{f}^2-\frac{3{\alpha }^2\left(1-c{r}^2\right)}{c}.$$

(102)

Redefining $\alpha$ as

$$\alpha =\frac{c}{\sqrt{1-c{r}_{\Sigma }^2}},$$

(103)

equations (98)–(102) become

$$A=f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}},$$

(104)

$$8\pi \mu =\frac{3{\dot{\kappa }}^2}{{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^2{\kappa }^2}+\frac{3c}{{\kappa }^2},$$

(105)

$$4\pi q=-\frac{cr\dot{\kappa }}{\sqrt{1-c{r}_{\Sigma }^2}{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^2{\kappa }^2},$$

(106)

$$8\pi P_r=8\pi P_{\perp }=-\frac{1}{{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^2}\left[\frac{2\ddot{\kappa }}{\kappa }-\frac{2\dot{f}\dot{\kappa }}{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)\kappa }+\frac{{\dot{\kappa }}^2}{{\kappa }^2}\right]+\frac{2c\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}}{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right){\kappa }^2}-\frac{c}{{\kappa }^2},$$

(107)

and

$$2\kappa \ddot{\kappa }-\frac{2\dot{f}\kappa \dot{\kappa }}{\left(f-1\right)}+{\dot{\kappa }}^2-\frac{2c{r}_{\Sigma }\dot{\kappa }}{\sqrt{1-c{r}_{\Sigma }^2}}=4fc-c{f}^2-3c.$$

(108)

Introducing the new variable

$$X\equiv \sqrt{c}(f-1),$$

(109)

(108) reads

$$2\kappa \ddot{\kappa }-\frac{2\dot{X}\kappa \dot{\kappa }}{X}+{\dot{\kappa }}^2-\frac{2c{r}_{\Sigma }\dot{\kappa }}{\sqrt{1-c{r}_{\Sigma }^2}}=-X^2+2\sqrt{c}X.$$

(110)

Next, using (30), (35) and (106) we obtain for the luminosity on the surface

$$L_{\Sigma }=-\frac{c{r}_{\Sigma }^3\dot{\kappa }}{\sqrt{1-c{r}_{\Sigma }^2}{(f-1)}^2},$$

(111)

or using (41), we obtain for the luminosity at infinity

$$L_{\infty }=-\frac{c{r}_{\Sigma }^3\dot{\kappa }}{\sqrt{1-c{r}_{\Sigma }^2}{(f-1)}^2}{\left(\sqrt{1-c{r}_{\Sigma }^2}+\frac{\dot{\kappa }{r}_{\Sigma }}{f-1}\right)}^2.$$

(112)

Also, observe that using (38) for this model, we obtain for the redshift at the boundary

$$z=\frac{(f-1)({\beta }_{\Sigma }-1)-\dot{\kappa }{r}_{\Sigma }{\beta }_{\Sigma }}{f-1+\dot{\kappa }{r}_{\Sigma }{\beta }_{\Sigma }},$$

(113)

and the time for the formation of a horizon is determined by the equation

$$\frac{\dot{\kappa }}{f-1}=-\frac{1}{{\beta }_{\Sigma }{r}_{\Sigma }}.$$

(114)

Thus, the model is completely determined up to two functions of t (f and $\kappa$). As mentioned before, in order to determine these two functions we have a large number of possible strategies. Here we shall resort to heuristic mathematical conditions, in order to fully determine the system.

As a first example we shall assume a heuristic mathematical condition on $\kappa$. Thus, we shall next consider the case where $\kappa$ has the linear form

$$\kappa ={\kappa }_{0}t+{\kappa }_1,$$

(115)

where ${\kappa }_0$ and ${\kappa }_1$ are arbitrary functions. Then, introducing (115) in (110) we obtain

$$\frac{2\dot{f}{\kappa }_0}{c(f-1)(f+b_1)(f+b_2)}=\frac{1}{{\kappa }_{0}t+{\kappa }_1},$$

(116)

whose solution is

$$\begin{array}{c}\hfill {(f-1)}^{b_1-b_2}{(f+b_1)}^{b_2+1}{(f+b_2)}^{-(b_1+1)}=\\ \hfill C{({\kappa }_{0}t+{\kappa }_1)}^{\frac{c(b_1+1)(b_2+1)(b_1-b_2)}{2{\kappa }_0^2}},\end{array}$$

(117)

where C is a constant and $b_1$ and $b_2$ have the following values

$$b_1=-2\pm \sqrt{1-\frac{{\kappa }_0{\kappa }_2}{c}},$$

(118)

$$b_2=-2\mp \sqrt{1-\frac{{\kappa }_0{\kappa }_2}{c}},$$

(119)

with

$${\kappa }_2\equiv {\kappa }_0-\frac{2c{r}_{\Sigma }}{\sqrt{1-c{r}_{\Sigma }^2}}.$$

(120)

In order to obtain f we have to solve the algebraic equation (117), for any given set of constants.

Thus, for example, for $b_1=0$, which implies $b_2=-4$, equation (117) reads

$$\begin{array}{c}\hfill \frac{{(f-1)}^4}{f^3(f-4)}=C{({\kappa }_{0}t+{\kappa }_1)}^{\frac{-6c}{{\kappa }_0^2}}.\end{array}$$

(121)

In general for the particular solution (117) the physical variables read

$$8\pi \mu =\frac{1}{{({\kappa }_{0}t+{\kappa }_1)}^2}\left[\frac{3{\kappa }_0^2}{{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^2}+3c\right],$$

(122)

$$4\pi q=-\frac{cr{\kappa }_0}{\sqrt{1-c{r}_{\Sigma }^2}{({\kappa }_{0}t+{\kappa }_1)}^2{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^2},$$

(123)

$$\begin{array}{ccc}\hfill 8\pi P_r=8\pi P_{\perp }& =& \frac{2{\kappa }_0\dot{f}}{({\kappa }_{0}t+{\kappa }_1){\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^3}-\frac{{\kappa }_0^2}{{({\kappa }_{0}t+{\kappa }_1)}^2{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}^2}\hfill \end{array}$$

(124)

$$\begin{array}{ccc}& +& \frac{1}{{({\kappa }_{0}t+{\kappa }_1)}^2}\left[\frac{2c\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}}{\left(f-\sqrt{\frac{1-c{r}^2}{1-c{r}_{\Sigma }^2}}\right)}-c\right],\hfill \end{array}$$

(125)

whereas for the luminosity we obtain

$$L_{\Sigma }=-\frac{c{r}_{\Sigma }^3{\kappa }_0}{\sqrt{1-c{r}_{\Sigma }^2}{(f-1)}^2}.$$

(126)

Observe that in this particular case the condition for the formation of the horizon as implied by (114) implies $f=constant$, which obviously contradicts (121). Thus no black hole results from the evolution of such a model.

As a second example we shall next consider the particular case $X=$ constant, for which (110) becomes

$$2\kappa \ddot{\kappa }+{\dot{\kappa }}^2-2ϵ\dot{\kappa }=\xi ,$$

(127)

where

$$ϵ\equiv \frac{c{r}_{\Sigma }^2}{\sqrt{1-c{r}_{\Sigma }^2}},\xi \equiv r_{\Sigma }^2(-X^2+2\sqrt{c}X),$$

(128)

and now dot denotes differentiation with respect to the dimensionless variable $t/r_{\Sigma }$.

By introducing the variable

$$\dot{\kappa }=z\Rightarrow \ddot{\kappa }=\dot{\kappa }\frac{dz}{d\kappa }=z\frac{dz}{d\kappa },$$

(129)

the equation above becomes

$$2\kappa \frac{dz}{d\kappa }+\frac{1}{\kappa }(z^2-2ϵz)=\frac{\xi }{\kappa },$$

(130)

whose solution reads

$$z\equiv \dot{\kappa }=\frac{{\xi }^{1/2}\sqrt{\kappa +h}}{\sqrt{\kappa }},$$

(131)

where

$$h=\frac{2}{\gamma }\left[ln\kappa \pm \sqrt{1+\gamma {\kappa }^2}\mp ln\left|\frac{1+\sqrt{1+\gamma {\kappa }^2}}{\kappa {\gamma }^{1/2}}\right|\right],$$

(132)

and $\gamma$ is an arbitrary constant.

We shall not elaborate further on these models, since the resulting expressions are too cumbersome, and our sole purpose here is to illustrate the way of using the proposed formalism, and not describe any specific astrophysical scenario.

7. A model obtained from Tolman VI as seed solution

Our next model is inspired in the well known Tolman $VI$ solution [52], whose equation of state for large values of $\mu$ approaches that for a highly compressed Fermi gas.

Thus we assume

$${\mu }_{eff}=\frac{g\left(t\right)}{r^2},$$

(133)

where g is an arbitrary (dimensionless) function of t. Using the above expression in (69) it follows

$$m_{eff}=4\pi {\kappa }^{3}g\left(t\right)r,$$

(134)

and replacing (134) into (68) we obtain

$$\frac{1}{{\beta }^2}=1-8\pi {\kappa }^{2}g\left(t\right)=1-c,$$

(135)

where c and $\beta$ are dimensionless constants.

Then using (65), (83), (86), (135), and redefining the constant $\alpha$ as

$$\alpha =\frac{2\sqrt{1-c}}{r_{\Sigma }^2},$$

(136)

the metric variables for this model read

$$\begin{array}{ccc}\hfill A& =& f\left(t\right)+{\left(\frac{r}{r_{\Sigma }}\right)}^2,\hfill \end{array}$$

(137)

$$\begin{array}{ccc}\hfill B& =& \frac{\kappa }{\sqrt{1-c}}=\beta \kappa ,\hfill \end{array}$$

(138)

$$\begin{array}{ccc}\hfill R& =& \kappa \left(t\right)r,\hfill \end{array}$$

(139)

and the expressions for the physical variables are

$$\begin{array}{ccc}\hfill 8\pi \mu & =& \frac{3{\dot{\kappa }}^2}{{\kappa }^2{(\frac{r^2}{r_{\Sigma }^2}+f)}^2}+\frac{{\beta }^2-1}{r^2{\kappa }^2{\beta }^2},\hfill \end{array}$$

(140)

$$\begin{array}{ccc}\hfill 4\pi q& =& -\frac{2r\dot{\kappa }}{{\kappa }^2\beta r_{\Sigma }^2{(\frac{r^2}{r_{\Sigma }^2}+f)}^2},\hfill \end{array}$$

(141)

$$\begin{array}{cc}8\pi P_r\hfill & =-\frac{1}{{(\frac{r^2}{r_{\Sigma }^2}+f)}^2}\left[\frac{2\ddot{\kappa }}{\kappa }-\frac{2\dot{f}\dot{\kappa }}{\kappa (\frac{r^2}{r_{\Sigma }^2}+f)}+\frac{{\dot{\kappa }}^2}{{\kappa }^2}\right]\hfill \\ & +\frac{4}{{\kappa }^2{\beta }^2{r}_{\Sigma }^2(\frac{r^2}{r_{\Sigma }^2}+f)}-\frac{{\beta }^2-1}{{\beta }^2{\kappa }^2{r}^2},\hfill \end{array}$$

(142)

$$\begin{array}{ccc}\hfill 8\pi (P_r-P_{\perp })& =& -\frac{{\beta }^2-1}{{\beta }^2{\kappa }^2{r}^2},\hfill \end{array}$$

(143)

whereas the junction condition, the luminosity and the redshift read

$$2\ddot{\kappa }\kappa -\frac{2\dot{f}\dot{\kappa }\kappa }{(f+1)}+{\dot{\kappa }}^2-4\frac{\dot{\kappa }}{\beta r_{\Sigma }}=\frac{4}{{\beta }^2{r}_{\Sigma }^2}\left(f+1\right)-\frac{{\beta }^2-1}{{\beta }^2{r}_{\Sigma }^2}{\left(f+1\right)}^2$$

(144)

$$L_{\Sigma }=-\frac{2r_{\Sigma }\dot{\kappa }}{\beta {(f+1)}^2},$$

(145)

$$L_{\infty }=-\frac{2r_{\Sigma }\dot{\kappa }{(f+1+\beta r_{\Sigma }\dot{\kappa })}^2}{{\beta }^3{(f+1)}^4},$$

(146)

and

$$z=\frac{(f+1)(\beta -1)-\dot{\kappa }{r}_{\Sigma }\beta }{f+1+\dot{\kappa }{r}_{\Sigma }\beta },$$

(147)

implying that the time for the formation of a horizon is determined by the equation

$$\frac{\dot{\kappa }}{f+1}=-\frac{1}{\beta r_{\Sigma }}.$$

(148)

It would be convenient to write (144) in terms of the dimensionless variable $\overline{t}\equiv t/r_{\Sigma }$, it reads

$$2\ddot{\kappa }\kappa -\frac{2\dot{f}\dot{\kappa }\kappa }{(f+1)}+{\dot{\kappa }}^2-4\frac{\dot{\kappa }}{\beta }=\frac{4}{{\beta }^2}\left(f+1\right)-\frac{({\beta }^2-1)}{{\beta }^2}{\left(f+1\right)}^2,$$

(149)

where now dots denote derivatives with respect to $\overline{t}$.

As in the precedent case we have a large number of possible strategies to obtain the two functions of t determining the whole system. Thus we could consider for example the $f=constant$ case, or the assumption of the linearity of $\kappa$. In both cases the procedure is very similar as in the preceding case. Instead, we shall propose a different approach here.

Specifically we shall split (149) in two equations, as follows

$$2\ddot{\kappa }\kappa +{\dot{\kappa }}^2-4\frac{\dot{\kappa }}{\beta }=0,$$

(150)

$$-\frac{2\dot{f}\dot{\kappa }\kappa }{(f+1)}=\frac{4}{{\beta }^2}\left(f+1\right)-\frac{({\beta }^2-1)}{{\beta }^2}{\left(f+1\right)}^2.$$

(151)

Equation (150) may be integrated producing

$$\frac{-2\omega b\sqrt{\kappa }+b^2\kappa +2{\omega }^{2}ln(\omega +b\sqrt{\kappa })}{b^3}=t+\gamma ,$$

(152)

where $\omega$ and $\gamma$ are two integration constants and $b\equiv 4/\beta$.

Solving the above transcendental equation for $\kappa$ and feeding the result back into (151) we obtain f.

Once the functions of time are determined, we have to resort to a transport equation (e.g. (12)) in order to find the distribution and evolution of the temperature.

As in the previous example, the resulting expressions are too burdensome and not very illuminating, so we shall not elaborate further on them.

8. Discussion and Conclusions

We have proposed an analytical approach to describe spherical collapse within the context of PQSR. To avoid the numerical integration of differential equations appearing in the algorithm put forward in [1–4], we have assumed the vanishing complexity factor as the cornerstone of the proposed method. Doing so, starting with a given “seed” static analytical solution to the Einstein equation, we are led to a situation where the whole system is determined by two arbitrary functions of t. These functions are related through the junction condition (33). For the additional information required to obtain the above mentioned functions, we have presented a list of possible strategies, based on either information obtained from observables such as luminosity and gravitational redshift, or from ad hoc heuristic mathematical conditions imposed on the system. It goes without saying that the presented list is not exhaustive, and much more possibilities can be considered. This last issue remains one the most important pending question regarding our approach.

Invoking the vanishing complexity factor as the main assumption behind the proposed approach is not arbitrary, and its rationale becomes intelligible when we remind that the complexity factor has been shown to be a good measure of the degree of complexity of a fluid distribution. Thus, assuming such a condition we ensure that we are dealing with the “simplest” fluid distributions available within the PQSR, in concord with one of the main goals of our endeavor consisting in describing gravitational collapse in its simplest possible way.

There is an additional argument reinforcing the assumption of vanishing complexity factor within the context of PQRS. Indeed, as we have seen, all models obtained with the approach here presented, are necessarily shear–free. On the other hand, as shown in [53], the shear–free condition is unstable in the presence of pressure anisotropy and/or dissipation. However, writing the complexity factor in terms of kinematical variables as

$$Y_{TF}=\frac{a^{\prime }}{B}-a\frac{R^{\prime }}{RB}+a^2-\frac{\dot{\sigma }}{A}-\frac{{\sigma }^2}{3}-\frac{2}{3}\Theta \sigma ,$$

(153)

it can be shown that the vanishing of the complexity factor implies the stability of the shear–free condition in the geodesic case (seen [53] for details). In the non–geodesic, static, case the combination of the first three terms on the right of (153) must be equal to zero if we assume the vanishing of the complexity factor, implying in its turn that such combination must remain non–vanishing but small (bounded) in the PQSR. In such a case we may safely conclude that the quasi–stability of $\sigma =0$ is ensured (see the discussion between Eqs.(63) and (67) in [53]).

Conditions on the complexity of the pattern of evolution such as H and $QH$, appear to be too strong and have to be excluded since they lead to geodesic fluids, which as mentioned before are physically incompatible with the very idea behind the PQSR.

Also, the adiabatic condition implies that the fluid is geodesic, accordingly we have considered exclusively dissipative fluids.

In order to illustrate the method we have presented two models. One is based on the interior Schwarzschild solution as the “seed” solution, whereas the other is inspired in the well known Tolman VI solution. The purpose of these calculations was to show how the algorithm works. In order to provide the missing information we have resorted to some mathematical ansatz. Of course it would be desirable to supply such information with physical data obtained from astrophysical observations, among which the luminosity and the gravitational redshift appear to be the most relevant. Such a task is out of the scope of this manuscript, but remains as the most important issue related to the discussion here presented.