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

Abstract

A seminumerical approach proposed many years ago for describing gravitational collapse in the post-quasi-static approximation is modified in order to avoid the numerical integration of the basic differential equations the approach is based upon. For doing that we have to impose some restrictions on the fluid distribution. More specifically, we shall assume the vanishing complexity factor condition, which allows for analytical integration of the pertinent differential equations and leads to physically interesting models. Instead, we show that neither the homologous nor the quasi-homologous evolution are acceptable since they lead to geodesic fluids, which are unsuitable for being described in the post-quasi-static approximation. Also, we prove that, within this approximation, adiabatic evolution also leads to geodesic fluids, and therefore, we shall consider exclusively dissipative systems. Besides the vanishing complexity factor condition, additional information is required for a full description of models. We shall propose different strategies for obtaining such an information, which are based on observables quantities (e.g., luminosity and redshift), and/or heuristic mathematical ansatz. To illustrate the method, we present two models. One model is inspired in the well-known Schwarzschild interior solution, and another one is inspired in Tolman VI solution.

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), in which case 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 implies that the fluid distribution changes on a time scale that is very long as compared to the hydrostatic time scale [1,2] (sometimes, this time scale is also referred to as the dynamical time scale, e.g., [3]). Thus, in this regime, the evolution of the fluid may be regarded as a sequence of static models, where the time between any two states of equilibrium is neglected (see [4,5,6] for applications).

The QSR applies to a large variety of scenarios due to the fact that the hydrostatic time scale is very small during many phases of the life of the star [2], e.g., it is of the order of 27 min, $4.5$ s and ${10}^{-4}$ s, respectively, for the Sun, a white dwarf and a neutron star of one solar mass and 10 Km radius.

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 that is smaller than the hydrostatic time scale.

All this having been said, the following question is in order: 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 [7,8,9,10] and references therein), the answer to the above questions is affirmative (in 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.

1.

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.

2.

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, it 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).

3.

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 [7,8,9,10].

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

To provide an accurate description of the gravitational collapse of a supermassive star, including the final fate of such process (naked singularities, black holes, anything else), the mechanism behind a type II supernova event [11,12,13,14,15,16,17] or the structure and evolution of the compact object resulting from such a process [18,19,20], is a task of utmost relevance.

We have available three approaches to study the gravitational collapse in the context of general relativity. On the one hand, numerical methods [21,22,23,24] allow for including more realistic equations of state. Nevertheless, the obtained results, in general, may be highly model-dependent. Moreover, difficulties associated with numerical solutions of partial differential equations in presence of shocks may complicate further the problem.

Alternatively, one may resort to analytical solutions to Einstein equations, which are more suitable for general discussions, and may be relatively simple to analyze, still containing some of the essential features of a realistic situation (see for example [25,26,27,28,29,30,31,32,33,34,35] an references therein). However, often they resort to heuristic assumptions, whose justification is unclear.

Between the two aforementioned approaches, we have seminumerical 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 [7,8,9,10] (see also [36,37]).

This third approach allows to reduce the initial system of partial differential equations into a system of ordinary differential equations (referred to as surface equations) for quantities evaluated at the boundary surface of the fluid distribution.

The approach relies on a set of conveniently defined variables (referred to as “effective” variables) plus an heuristic ansatz on the latter, 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 carrying that out 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 3, 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 4. A review of the approach is outlined in Section 5, and some examples are analyzed in Section 6. 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,42,43,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$ is the radial pressure, $P_{\perp }$ is the tangential pressure, $q^{\alpha }$ is the heat flux, $V^{\alpha }$ is the four-velocity of the fluid and $K^{\alpha }$ is 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) to be 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,46,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)

At this point, the following remark is in order: The knowledge of $A(t,r)$, $B(t,r)$ and $R(t,r)$ casts the system above in an algebraic system of four equations for the four unknown functions $\mu$, q, $P_r$ and $P_{\perp }$ which, in such a case, 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 to say, at this point, the first term on the right-hand side describes the gravitational force term.

3. The Exterior Spacetime and Junction Conditions

Outside $\Sigma$, we assume that we have the Vaidya spacetime (i.e., we assume that 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 that (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 Equation (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 that from (31), (36) and (39), it follows that

$$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 that has been proposed in order 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)

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,\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 [1,2,3].

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 nonrelativistic 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 geodesic 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 obtain

$${\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 the following:

$$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 following 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 “dynamic 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 as follows:

$$\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 that it can be considered to be in equilibrium at each moment (Equation (57) is satisfied).

This implies the following for U, the metric and the kinematical functions:

• The areal velocity U and the 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 the second-order 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 is that these variables depend upon time according to Equation (15).

5.3. Post-Quasi-Static Regime (PQSR)

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 Equation (57)? Such a situation is described by what we call PQSR.

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 [8,9,10] 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, we see at once that if $R=\kappa \left(t\right)r$ and ${\mu }_{eff}$ have the same radial dependence as $\mu$ in the static case, then the radial dependence of $m_{eff}$ will be the same as in the static case. 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 [10]. 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 [10].

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

1.

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).$$

2.

Assume that the r dependence of the effective density is the same as that of ${\mu }_{st}$, and $R=r\kappa \left(t\right)$.

3.

Impose the vanishing complexity factor condition.

4.

From the two conditions above, we are able to determine the metric functions up to two arbitrary functions of t.

5.

For these functions of t, one has the junction condition (33).

6.

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 for 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 (33) 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)

where $\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\left(t\right),$$

(86)

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

Then, Equations (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=\mathrm{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, the 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.

6.2. 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 that

$$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) and (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 \\ \hfill 8\pi P_r& =& -\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 \end{array}$$

(141)

$$\begin{array}{ccc}& +& \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, the 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=\mathrm{constant}$ case, or the assumption of the linearity of $\kappa$. In both cases, the procedure is very similar to the preceding case. Instead, we shall propose a different approach here.

Specifically, we shall split (149) into 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.

7. 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 [7,8,9,10], we have assumed the vanishing complexity factor as the cornerstone of the proposed method. As far as we are aware, this is the first approach for modeling gravitational collapse that includes both the PQSR and the vanishing complexity conditions. 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. In this work, and with the sole purpose to illustrate the method, we have resorted to heuristic mathematical restrictions. It must be clear that the full potential of the approach may only be deployed when the missing information is provided by either of the observables mentioned above. Although this last issue remains one the most important pending question regarding our approach, it is out of the scope of this manuscript.

Invoking the vanishing complexity factor as the main assumption behind the proposed approach is not arbitrary, and its rationale becomes intelligible when we remember 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 of describing gravitational collapse in the simplest possible way.

There is an additional argument reinforcing the assumption of a 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 nonvanishing 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 Equations (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. We would like to emphasize once again that the optimal path to display the power of the presented method would be to supply such information through physical data obtained from astrophysical observations, among which the luminosity and the gravitational redshift appear to be the most relevant. We harbor the hope that some of our colleagues will be able to succeed in such an endeavor.