The post--quasi-static approximation: An analytical approach to gravitational collapse

A semi--numerical 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.


I. 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 [5]- [6] (sometimes this time scale is also referred to as dynamical time scale, e.g.[7]).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 [8][9][10] 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 [6], e.g. it is of the order of 27 minutes, 4.5 seconds and 10 −4 seconds, 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 which 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 [1][2][3][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 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 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 [1][2][3][4] 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]- [17] or the structure and evolution of the compact object resulting from such a process [18]- [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]- [24], which allow for including more realistic equations of state.Nevertheless, the obtained results, in general, may be highly model dependent.Besides, difficulties, associated to 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] and references therein).However, often they resort to heuristic assumptions whose justification is unclear.
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][2][3][4] (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 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.

A. The metric
We consider a spherically symmetric distribution of collapsing fluid, bounded by a spherical surface Σ.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 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 = θ and x 3 = φ.

B. Energy-momentum tensor
The matter energy-momentum tensor T αβ inside Σ has the form where µ is the energy density, P r the radial pressure, P ⊥ the tangential pressure, q α the heat flux, V α the fourvelocity of the fluid, and K α a unit four-vector along the radial direction.These quantities satisfy Since we assume the metric (1) comoving then where q is a function of t and r.

C. Kinematical variables
The four-acceleration a α and the expansion Θ of the fluid are given by and its shear σ αβ by where 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 ⊥ , of the collapsing fluid.
From (4) with (3) we have for the four-acceleration and its scalar a, where a α = aK α , and for the expansion 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 and its scalar where Then, the shear tensor can be written as

D. 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 where k, T and τ 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 α , it reads

E. Field equations
The Einstein field equations for the interior spacetime (1) can be written as 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 µ, q, P r , and P ⊥ which, in such a case, can be obtained without further information.

F. Mass and areal velocity
Following Misner and Sharp [48], let us now introduce the mass function m(t, r) (see also [49]), defined by It is useful to introduce the proper time derivative D T given by and the proper radial derivative D R , where R defines the areal radius of a spherical surface inside Σ (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.
Then ( 18) can be rewritten as Using ( 14)-( 16) with ( 19) and ( 20) we obtain from ( 18) and Next, the three-acceleration D T U of an in-falling particle inside Σ can be obtained by using ( 16), ( 18) and ( 22), producing or Finally, from the Bianchi identities we obtain 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.

III. THE EXTERIOR SPACETIME AND JUNCTION CONDITIONS
Outside Σ we assume we have the Vaidya spacetime (i.e.we assume all outgoing radiation is massless), described by where M (v) 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 Σ = constant, was discussed in [50].Now, from the continuity of the first differential form it follows (see [50] for details), and where τ denotes the proper time measured on Σ.
The continuity of the second differential form produces = means that both sides of the equation are evaluated on Σ (observe a misprint in eq.( 40) in [50] and a slight difference in notation).
Comparing (33) with (15) and ( 16) one obtains Thus the matching of ( 1) and ( 28) on Σ implies (32) and (34).Also, we have where L Σ denotes the total luminosity of the sphere as measured on its surface and is given by and where is the total luminosity measured by an observer at rest at infinity.The boundary redshift z Σ is given by with Therefore the time of formation of the black hole is given by Also observe than from (31), (36) and ( 39) it follows and from ( 21), ( 22), ( 31) and ( 39) IV. 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 T F which defines the trace-free part of the electric Riemann tensor (see [51] for details).
Thus, let us define tensor Y αβ by which may be expressed in terms of two scalar functions Y T , Y T F , as Then after lengthy but simple calculations, using field equations, we obtain (see [39,40] for details) In terms of the metric functions the scalar Y T F reads A. 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 and where R I and R II denote the areal radii of two concentric shells (I, II) described by r = r I = constant, and r = r II = constant, respectively.These relationships are reminiscent of the homologous evolution in Newtonian hydrodynamics [5][6][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 which implies (48) only if the fluid is geodesic (A = 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) 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 and combining with condition (50), we obtain whereas, using ( 7) and (10) we get But in the PQSR we have (see equation (65) in section 5.3 below) R = κ(t)r where κ is an arbitrary function of t, producing at once that 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.

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

A. Static regime
In this case all time derivatives vanish, implying: Since B = B(r); A = A(r); R = R(r), reparametrizing r, we may write the line element in the form: Thus, the "dynamic equation ( 27) becomes the well known TOV equation of hydrostatic equilibrium for an anisotropic fluid The Einstein equations in this case read: Also, for the mass function we have or and for the metric function A, we have from ( 26) The important point to keep in mind is that if the radial dependence of µ and P r is known, the metric functions are determined from (61-63).
B. 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).
This implies for U , the metric and the kinematical functions that • The areal velocity U and the kinematical variables are small, (of order O(ǫ), with |ǫ| << 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 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).

C. 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 eq.(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 where κ is an arbitrary (dimensionless) function of t, to be determined later.Taking into account ( 22) and (65), we rewrite the metric as follows Next, defining the effective mass as we obtain Then, equations ( 24) and ( 26) can be written as with where we have followed the terminology used in [2][3][4] and call µ ef f and P ef f 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 From the above, we see at once that if R = κ(t)r and µ ef f have the same radial dependence as µ in the static case, then the radial dependence of m ef f will be the same as in the static case.
On the other hand, if besides the assumption above, we assume that P ef f 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.
1. Take an interior ("'seed") solution to Einstein equations, representing a fluid distribution of matter in equilibrium, with a given µ st = µ st (r); P rst = P rst (r).
2. Assume that the r dependence of the effective density is the same as that of µ st , and R = rκ(t).
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 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 Σ ).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).

VI. 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 T F .Using ( 14)-( 17), ( 46) and (65), we obtain and Let us first consider the q = 0 case, which using (76) produces Since at r = 0, A is different from zero, we must impose and Since the geodesic case in the PQSR should be dismissed by reasons exposed before, we shall consider exclusively dissipative systems.
Then since q = 0, it follows from (76) that B is separable here β is an arbitrary dimensionless function of r, and It is worth stressing that using (83) in ( 10) it follows at once that σ = 0. Thus all our models will be shear-free.
Next, assuming Y T F = 0 we obtain from (79) whose solution reads where f is arbitrary function of integration, and, by reparametrizying t, another function of integration has been put equal to Then eqs.( 75)-( 78) take the form where A is given by (86).Also, from (89) and ( 90) Using ( 65) and (83) we can write 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.

A. 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 µ ef f = F (t), where F (t) is an arbitrary function with units [1/r 2 ], we obtain from (73), and with ( 22) and (68) we have then where c is a constant, with the same units as F (t), given by

With this we have for
and for the field equations On the surface Σ, from (33) or (34) we obtain Redefining α as equations ( 98)-(102) become Introducing the new variable Next, using ( 30), ( 35) and ( 106) we obtain for the luminosity on the surface or using (41), we obtain for the luminosity at infinity (112) Also, observe that using (38) for this model, we obtain for the redshift at the boundary and the time for the formation of a horizon is determined by the equation Thus, the model is completely determined up to two functions of t (f and κ).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 κ.Thus, we shall next consider the case where κ has the linear form where κ 0 and κ 1 are arbitrary functions.Then, introducing (115) in (110) we obtain whose solution is where C is a constant and b 1 and b 2 have the following values with 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 In general for the particular solution (117) the physical variables read whereas for the luminosity we obtain 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 where and now dot denotes differentiation with respect to the dimensionless variable t/r Σ .By introducing the variable the equation above becomes whose solution reads where and γ 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.Our next model is inspired in the well known Tolman V I solution [52], whose equation of state for large values of µ approaches that for a highly compressed Fermi gas.
Thus we assume where g is an arbitrary (dimensionless) function of t.Using the above expression in (69) it follows and replacing (133) into (68) we obtain 1 where c and β are dimensionless constants.
Then using (65), ( 83), ( 86), (134), and redefining the constant α as the metric variables for this model read and the expressions for the physical variables are 8πµ whereas the junction condition, the luminosity and the redshift read implying that the time for the formation of a horizon is determined by the equation It would be convenient to write (143) in terms of the dimensionless variable t where now dots denote derivatives with respect to 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 κ.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 (148) in two equations, as follows Equation ( 149) may be integrated producing where ω and γ are two integration constants and b ≡ 4/β.Solving the above transcendental equation for κ and feeding the result back into (150) 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.

VII. 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][2][3][4], 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 which 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 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 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 (152) 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 σ = 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.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 endeavor.

B
. A model obtained from Tolman VI as seed solution