Non-static fluid spheres admitting a conformal Killing vector: Exact solutions

We carry on a general study on non--static spherically symmetric fluids admitting a conformal Killing vector (CKV). Several families of exact analytical solutions are found for different choices of the CKV, in both, the dissipative and the adiabatic regime. To specify the solutions, besides the fulfillment of the junction conditions on the boundary of the fluid distribution, different conditions are imposed, such as vanishing complexity factor and quasi--homologous evolution. A detailed analysis of the obtained solutions, its prospective applications to astrophysical scenarios, as well as alternative approaches to obtain new solutions, are discussed.


I. INTRODUCTION
The purpose of this work is twofold. On the one hand we want to delve deeper into the physical consequences derived from the assumption that a given space-time admits a CKV. This interest in its turn is motivated by the relevance of such kind of symmetry in hydrodynamics.
Indeed, in general relativity, self-similar solutions are related to the existence of a homothetic Killing vector field (HKV), a generalization of which is a conformal Killing vector field (CKV). The physical interest of systems admitting a CKV is then suggested by the important role played by self-similarity in classical hydrodynamics.
Thus, in Newtonian hydrodynamics, self-similar solutions are those described by means of physical quantities which are functions depending on dimensionless variables x/l(t), where x and t are independent space and time variables and l is a time dependent scale. Therefore the spatial distribution of the characteristics of motion remains similar to itself at all times [1]. In other words, self-similarity is to be expected whenever the system under consideration possesses no characteristic length scale. The above comments suggest that self-similarity plays an important role in the study of systems close to the critical point, where the correlation length becomes infinite, in which case different phases of the fluid (e.g. liquid-vapor) may coexist, the phase boundaries vanish and density fluctuations occur at all length scales. This process may be observed in the critical opalescence.
Besides, examples of self-similar fluids may be found in the study of strong explosions and thermal waves [2][3][4][5].
Motivated by the above arguments many authors, since the pioneering work by Cahill and Taub [6], have focused their interest in the problem of self-similarity in selfgravitating systems. Some of them are restricted to general relativity, with especial emphasis on the ensuing consequences from the existence of HKV or CKV, and possible solutions to the Einstein equations (see for example  and references therein). Also, a great deal of work has been done in the context of other theories of gravitation (see for example [49][50][51][52][53][54][55][56][57][58][59][60][61] and references therein). Finally, it is worth mentioning the interest of this kind of symmetry related to the modeling of wormholes (see [62][63][64][65][66][67][68] and references therein).
On the other hand, the problem of general relativistic gravitational collapse has attracted the attention of researchers since the seminal paper by Oppenheimer and Snyder. The origin of such interest resides in the fact that the gravitational collapse of massive stars represents one of the few observable phenomena where general relativity is expected to play a relevant role. To tackle such a problem there are two different approaches: Numerical methods or analytical exact solutions to Einstein equations. Numerical methods enable researchers to investigate systems that are extremely difficult to handle analytically. However, purely numerical solutions usually hinder the investigation of general, qualitative, aspects of the process. On the other hand, analytical solutions although are generally found either for too simplistic equations of state and/or under additional heuristic assumptions whose justification is usually uncertain, are more suitable for a general discussion and seem to be useful to study non-static models which are relatively simple to analyze but still contain some of the essential features of a realistic situation.
In this manuscript we endeavor to find exact, analytical, non-static solutions admitting a CKV, including dissipative processes. The source will be represented by an anisotropic fluid dissipating energy in the diffusion approximation. In order to find the solutions we shall specialize the CKV to be either space-like (orthogonal to the four-velocity), or time-like (parallel to the fourvelocity). In each case we shall consider separately the dissipative and non-dissipative regime. Also, in order to specify the models, we will assume specific restrictions on the mode of the evolution, (e.g. the quasi-homologous condition), and on the complexity factor, among other conditions. A fundamental role in finding our models is played by the equations ensuing from the junction conditions on the boundary of the fluid distribution, whose integration provides one of the functions defining the metric tensor.
Several families of solutions are found and discussed in detail. A summary of the obtained results and a discussion on the physical relevance of these solutions are presented in last section. Finally several appendices are included containing useful formulae.

II. THE METRIC, THE SOURCE AND RELEVANT EQUATIONS AND VARIABLES
In what follows we shall briefly summarize the definitions and main equations required for describing spherically symmetric dissipative fluids. We shall heavily rely on [69], therefore we shall omit many steps in the calculations, details of which the reader may find in [69].
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 (diffusion approximation).
The justification to consider anisotropic fluids is provided by the fact that pressure anisotropy is produced by many different physical phenomena of the kind expected in gravitational collapse scenario (see [70] and references therein). Furthermore we expect that the final stages of stellar evolution should be accompanied by intense dissipative processes, which, as shown in [71], should produce pressure anisotropy.
Choosing comoving coordinates, the general interior metric can be written as 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 = φ. Observe that A and B are dimensionless, whereas R has the same dimension as r.
The energy momentum tensor in the canonical form, reads where µ is the energy density, P r the radial pressure, P ⊥ the tangential pressure, q α = qK α the heat flux, V α the four-velocity of the fluid, and K α a unit four-vector along the radial direction. Since we are considering comoving observers, we have These quantities satisfy It is worth noticing that we do not explicitly add bulk or shear viscosity to the system because they can be trivially absorbed into the radial and tangential pressures, P r and P ⊥ , of the collapsing fluid (in Π). Also we do not explicitly introduce dissipation in the free streaming approximation since it can be absorbed in µ, P r and q.
The Einstein equations for (1) and (2), are explicitly written in Appendix A.
The acceleration a α and the expansion Θ of the fluid are given by and its shear σ αβ by From (5) we have for the four-acceleration and its scalar a, and for the expansion where the prime stands for r differentiation and the dot stands for differentiation with respect to t. We obtain for the shear (6) its non zero components and its scalar where Next, the mass function m(t, r) reads Introducing the proper time derivative D T given by we can define the velocity U of the collapsing fluid as the variation of the areal radius with respect to proper time, i.e.
where R defines the areal radius of a spherical surface inside the fluid distribution (as measured from its area). Then (12) can be rewritten as Using (15) we can express (A6) as where D R denotes the proper radial derivative, Using (A2)-(A4) with (13) and (17) we obtain from (12) and which implies satisfying the regular condition m(t, 0) = 0.
Integrating (20) we find A. The Weyl tensor and the complexity factor Some of the solutions exhibited in the next section are obtained from the condition of vanishing complexity factor. This is a scalar function intended to measure the degree of complexity of a given fluid distribution [72,73], and is related to the so called structure scalars [74].
In the spherically symmetric case the magnetic part of the Weyl tensor (C ρ αβµ ) vanishes, accordingly it is defined by its "electric" part E γν , defined by whose non trivial components are where Observe that the electric part of the Weyl tensor, may be written as As shown in [72,73] the complexity factor is identified with the scalar function Y T F which defines the trace-free part of the electric Riemann tensor (see [74] 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 [75] for details) Next, using (A2), (A4), (A5) with (12) and (24) we obtain 3m which combined with (21) and (28) produces It is worth noticing that due to a different signature, the sign of Y T F in the above equation differs from the sign of the Y T F used in [72] for the static case.
Thus the scalar Y T F may be expressed through the Weyl tensor and the anisotropy of pressure or in terms of the anisotropy of pressure, the density inhomogeneity and the dissipative variables.
In terms of the metric functions the scalar Y T F reads B. The exterior spacetime and junction conditions Since we are considering bounded fluid distributions then we still have to satisfy the junction (Darmois) conditions. Thus, outside Σ we assume we have the Vaidya spacetime (i.e. we assume all outgoing radiation is mass-less), described by where M (v) denotes the total mass, and v is the retarded time.
The matching of the full nonadiabatic sphere to the Vaidya spacetime, on the surface r = r Σ = constant, requires the continuity of the first and second fundamental forms across Σ (see [76] and references therein for details), which implies where Σ = means that both sides of the equation are evaluated on Σ.
Finally, the total luminosity (L ∞ ) for an observer at rest at infinity is defined by

III. THE TRANSPORT EQUATION
In the dissipative case we shall need a transport equation in order to find the temperature distribution and evolution. Assuming a causal dissipative theory (e.g.the Israel-Stewart theory [77][78][79] ) the transport equation for the heat flux reads where κ denotes the thermal conductivity, and T and τ denote 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 (37) by contracting with the unit spacelike vector K α , it reads (38) Sometimes it is possible to simplify the equation above, in the so called truncated transport equation, when the last term in (37) may be neglected [80], producing

IV. THE HOMOLOGOUS AND QUASI-HOMOLOGOUS CONDITIONS
As mentioned before, in order to specify some of our models we shall impose the condition of vanishing complexity factor. However, for time dependent systems, it is not enough to define the complexity of the fluid distribution. We need also to elucidate what is the simplest pattern of evolution of the system.
In [73] the concept of homologous evolution was introduced, in analogy with the same concept in classical astrophysics, as to represent the simplest mode of evolution of the fluid distribution.
Thus, the field equation (A3) written as can be easily integrated to obtain whereã is an integration function, or If the integral in the above equations vanishes we have from (41) or (42) that This relationship is characteristic of the homologous evolution in Newtonian hydrodynamics [81][82][83]. In our case, this may occur if the fluid is shear-free and non dissipative, or if the two terms in the integral cancel each other.
In [73], the term "homologous evolution" was used to characterize relativistic systems satisfying, besides (43), the condition where R 1 and R 2 denote the areal radii of two concentric shells (1, 2) described by r = r 1 = constant, and r = r 2 = constant, respectively. The important point that we want to stress here is that (43) does not imply (44). Indeed, (43) implies that for the two shells of fluids 1, 2 we have that implies (44) only if A = A(t), which by a simple coordinate transformation becomes A = constant. Thus in the non-relativistic regime, (44) always follows from the condition that the radial velocity is proportional to the radial distance, whereas in the relativistic regime the condition (43) implies (44), only if the fluid is geodesic.
In [69] the homologous condition was relaxed, leading to what was defined as quasi-homologous evolution, restricted only by condition (43), implying

V. CONFORMAL MOTIONS: EXACT SOLUTIONS
We shall consider spacetimes whose line element is defined by (1), admitting a CKV, i.e. satisfying the equation where L χ denotes the Lie derivative with respect to the vector field χ, which unless specified otherwise, has the general form and ψ in principle is a function of t, r. The case ψ = constant corresponds to a HKV. Our goal consists in finding exact solutions admitting a one parameter group of conformal motions, expressed in terms of elementary functions.
Two different families of solutions will be obtained depending on the choice of χ α . One of these families corresponds to the case with χ α orthogonal to V α , while the other corresponds to the case with χ α parallel to V α . For both families we shall consider separately the nondissipative (q = 0) and the dissipative (q =) case.
For the non-dissipative case of the family of solutions with χ α orthogonal to V α , we shall obtain from the matching conditions and specific values of the relevant parameters, solutions I, II, III, and for the particular case M = 0 we shall obtain solutions IV, V . For the dissipative case of this family, imposing the vanishing complexity factor condition and the shear-free condition we shall obtain solution V I.
For the non-dissipative case of the family of solutions with χ α parallel to V α , we shall obtain from the matching conditions and the vanishing complexity factor condition solution V II, whereas from specific values of relevant parameters we shall obtain solution V III. Also imposing the condition M = 0 we shall obtain in this case solutions IX, X.
Finally for the dissipative case of this family, imposing the complexity factor condition, we shall obtain solution XI.
Let us start by considering the case χ α orthogonal to V α and q = 0.
Then from and From (50) and (52) it follows where h is an arbitrary function of t, which without loos of generality may be put equal to 1 by reparametrizing t.
Thus we may write where α is a unit constant with dimensions of 1 length . Next, taking the time derivative of (51) and (52) and using (53) we obtain where G(r) is an arbitrary function of r which may be put equal to 1 by a a reparametrization of r, and F 1 is an arbitrary dimensionless function of t.
Thus we have and Then, feeding back (55) and (57) into (A3) with q = 0, one obtains where f and g are two arbitrary functions of their arguments and F (t) ≡ 1 F1(t) . So far we see that any model is determined up to three arbitrary functions F (t), f (t), g(r).
Then the field equations read Using the results above the matching conditions (33) and (35) on the surface r = r Σ = constant reaḋ and with ω ≡ g ′ (r Σ ) 2 . It is a simple matter to check that (63) is just the first integral of (64), therefore we only need to consider the former equation.
It would be useful to write (63) in the forṁ with orṼ Obviously all solutions have to satisfy the conditions ω > V (R Σ ). Among them we have: In this case we may have solutions evolving between the singularity and some value of R Σ in the interval [2M, 3M ] (region A in figure 1), and solutions with R Σ in the interval [3M, ∞] (region B in figure 1). figure 1).
In general, we may write from (65) from which we may obtain R Σ expressed in terms of elliptic functions. However, in some cases analytical solutions may be found in terms of elementary functions. For doing that we shall proceed as follows. Let us introduce the variable r = x + b in the polynomial which allows us to write or where b 1,2,3 and b are solutions of the following equations Then the integration of (68) produces where To obtain explicit solutions expressed through elementary functions, we shall assume b = 0, thus in our notation we have Imposing b = 0, we are led to two sub-cases, b 2 = 0, or b 3 = 0, in both sub-cases M = 1 3 √ 3ω . Using (79)- (81) in (76), we obtain for both sub-cases the same solutions, namely and In the first case, the areal radius of the boundary (R (I) Σ ) expands from 0 (the singularity) approaching 3M asymptotically as t → ∞, thereby representing a white hole scenario.
In the second case, the areal radius of the boundary (R Thus the fulfillment of the matching conditions provides one of the arbitrary functions of time describing our metric. In order to specify further our model we shall impose the quasi-homologous evolution and the vanishing complexity factor condition. As we can see from (42), in the non-dissipative case the quasi-homologous condition implies that the fluid is shear-free (σ = 0), implying in its turṅ Thus the metric functions become (85) Therefore our models are now specified up to an arbitrary function of r (g(r)). In order to fix this function we shall further impose the vanishing complexity factor condition. Then feeding back (85) into (31) we obtain Using (85) in (86), it follows at once with c 1 ≡ ± √ ω, and c 2 is another integration constant, we shall choose the negative sign in c 1 in order to ensure that R ′ > 0. However it should be noticed that the regularity conditions, necessary to ensure elementary flatness in the vicinity of the axis of symmetry, and in particular at the center (see [84], [85], [86]), are not satisfied.
Therefore after the imposition of the two conditions above (quasi-homologous evolution and vanishing complexity factor) we have all the metric functions completely specified for any of the above solutions to (65).
Thus in the case b = 0 we obtain from (82) from which the physical variables are easily found to be From (90) it follows at once that (P r ) Σ = 0. It is worth noticing that the expansion scalar for this model reads Thus the expansion is homogeneous and positive, diverging at t = t 0 , and tending to zero as t → ∞. The fast braking of the expansion for t > t 0 is produced by the negative initially large (diverging at t = t 0 ) value of D T U . This can be checked from (B6), where the negative gravitational term proportional to m/R 2 provides the leading term in the equation ( m R 2 ∼ equilibrium is reached asymptotically, but not, as usual, by the balance between the gravitational term (the first term on the right of (B6)) and the hydrodynamic terms (the second term on the right of (B6)). Instead, both terms cancel independently. Indeed, as t → ∞, the gravitational term vanishes due to the fact that the inertial mass density (the "passive gravitational mass density") µ + P r → 0, and the hydrodynamic term vanishes because, as it can be easily checked, the radial pressure gradient cancels the anisotropic factor, as t → ∞. Next, if we take (83) we obtain for f (t) whereas the expressions for the physical variables read In the limit t → ∞ the two above solutions converge to the same static distribution whose physical variables are 8πµ = 1 where the constant F 0 has been chosen F 0 = rΣα √ 3 . It is worth noticing that the ensuing equation of state for the static limit is the Chaplygin-type equation µ = −P r .
In the case ω = 0 the expression for R Σ is given by and from(87) Then the following expressions may be obtained for the metric functions and the physical variables read It is worth stressing the presence of important topological pathologies in this solution (e.g. R ′ = 0), implying the appearance of shell crossing singularities.
Before closing this subsection we would like to call the attention to a very peculiar solution that may be obtained by assuming that the space-time outside the boundary surface delimiting the fluid is Minkowski. This implies M = 0, and then the solutions to (63) read and Assuming further that the evolution is quasihomologous and the complexity factor vanishes, we obtain for the functions f (t) and The corresponding physical variables for f (IV ) read whereas for f (V ) they are (114) In the above the constants F 0 , α, r Σ have been chosen such that F0 αrΣ = 1. This kind of configurations have been considered in [11,87].
Let us now consider the general dissipative case when the vector χ α is orthogonal to the four-velocity.
Then from (49) we obtain, following the same procedure as in the non-dissipative case where α is a unit constant with dimensions of 1 length , where F (t) is and arbitrary function of t, and Then, feeding back (115) and (116) into (A3) with q = 0, one obtainsḂ The equation above may be formally integrated, to obtain where f and g are two arbitrary functions of their arguments.
In order to find a specific solutions we shall impose next the vanishing complexity factor condition (Y T F = 0).
Then from the above expressions and (31), the condition Y T F = 0 reads In order to find a solution to the above equation we shall assume that and The integration of (123) produces where β, γ are arbitrary functions of t. It is worth noticing that β has dimensions of 1 length , and γ is dimensionless.
Next, taking the r derivative of (122) we obtain γ = β α . Then we may write Next, combining (122) with (125) we obtain where c 3 , c 4 are arbitrary constants. From the above expression it follows at oncė On the other hand, (119) with (125), imply 4π qAdtdr = β(t)r, g(r) = g 1 r + g 0 , where g 1 , g 0 are constant. We may now write the physical variables in terms of the function β(t), they read 8πµ = β 2 (αr + 1) 2 α 2 c 2 3 e 2c4 βdt α 2 − 4c 4β + (132) The function β may be found, in principle, from the junction condition (35), however since this is in practice quite difficult at this level of generality, we shall first impose further constraints on our fluid distribution in order to obtain a simpler model, and afterwards we shall use the junction conditions. We shall start by imposing the quasi-homologous condition (46). Then using (119) and (120) in (46) we get Using (133) with (125)-(128) one obtains where c 5 is a constant with dimensions of length. So the metric functions may be written as It is worth noticing that the areal radius is independent on time (U = 0), solutions of this kind have been found in [69] Next, instead of quasi-homologous condition we shall impose the shear-free condition. Then assuming σ = 0 it follows at once thatḞ = 0 implying c 4 = 0. Then the metric functions become from which we can write the physical variables as Now we can find β from the junction condition (35), which using (139) and (140) with In order to integrate the above equation, let us introduce the variable s =β β , which casts (142) into the Ricatti equation whose solution is producing for β where α 2 is a negative constant of integration with the same dimensions as β.
Using the truncated version of the transport equation (39), we obtain for the temperature where c 3 and T 0 (t) are arbitrary constant and function of integration, respectively. The model described by equations (137)-(141) and (146), (147) will be named as model V I.
We shall next analyze the case when the vector χ is parallel to the four-velocity vector. We start by considering the non-dissipative case. In this case the equation (49) produces where h(r) is an arbitrary function of its argument. It is worth noticing that in this case the fluid is necessarily shear-free.
Thus the line element may be written as (149) Next, using (148) in (A3), the condition q = 0 readṡ whose solution is implying where g, f are two arbitrary functions of their argument.
Thus the metric is defined up to three arbitrary functions (g(r), f (t), h(r)).
The function f (t) will be obtained from the junction conditions (33), (35).
Indeed, evaluating the mass function at the boundary surface Σ we obtain from (33) and (151) where hΣ . On the other hand, from (35), using (151) we obtain To specify a model we have to obtain f (t) from the solution to the above equations.
In the special case a 1 = 1 (153) becomeṡ which has exactly the same form as (65) and therefore admits the same kind of solutions, and (155) reads a first integral of which, as it can be easily shown, is (156), therefore we only need to satisfy (156). In order to determine the functions g(r), h(r) we shall assume the vanishing complexity factor condition Y T F = 0.
Using (151) in (31) the condition Y T F = 0 reads or producing where c 4 , c 5 are arbitrary constants. If we choose implying a 1 = 1, then we obtain from (161) where c 6 , c 7 are constants. Thus, let us consider the following model. The time dependence described by f (t) is obtained from the solution to (156) given by with ǫ = 1 27M 2 , and the radial dependence of the model is given by the functions g(r), h(r) given by (162) (163).
The physical variables corresponding to this model read where the following relationships between the constants has been used α 2 ≡ c 2 In the limit t → ∞ the above model tends to a static fluid distribution described by satisfying the equation of state P r = −µ.
Another case which allows integration in terms of elementary function may be obtained from the conditions ǫ = 0 and a 1 = 1/2. Then (153) readṡ The above equation may be easily integrated, producing witht ≡ √ 3α 2 (t − t 0 ). Next, in order to specify further the model, we shall impose the vanishing complexity factor condition. In this case (a 1 = 1/2), the general solution to (159) reads however since ǫ = 0 the constant c 2 must vanish. The physical variables for this model read This solution represents fluid distribution oscillating between R Σ = 0 and R Σ = 8M 3 . It is worth noticing that the energy density is always positive, whereas the radial pressure is not.
Finally, we shall present two solutions describing a "ghost" compact object, of the kind already discussed in the previous section.
Thus assuming M = 0, equation (153) becomeṡ Solutions to the above equation in terms of elementary functions may be obtained by assuming a 1 = 1, in which case the two possible solutions to (177) are and Imposing further the vanishing complexity factor condition, then functions h(r), g(r) are given by (162) and (163). The physical variables corresponding to (178) and (179) read respectively and Finally, we shall consider the case where the CKV is parallel to the four-velocity, and the system is dissipative. As result of the admittance of the CKV the metric functions read as (148). Then, feeding this back into (A3) producesḂ which may be formally integrated, to obtain implying and where f (t) and g(r) are arbitrary functions of their arguments. To specify a model we shall impose the vanishing complexity factor condition. Thus, using (187)-(189) in (31) the condition Y T F = 0 reads a formal integration of which produces where γ(t) is an arbitrary function. Also, taking the t-derivative of (190) we obtain Using the above expressions we may write the metric functions (187)-(189) as implying A = 1 and Further restrictions on functions f (t), γ(t) will be obtained from the junction condition (P r = q) Σ .
Indeed, using (192)-(195) and (A4), the condition (P r = q) Σ , reads where In order to solve the above equation we shall assume and where a 1 ≡ h ′ Σ rΣ hΣ . From (199) it follows at once that producing Using (200) and (201) in (198), this last equation becomes In order to integrate (202), let us introduce the variable y =Ẋ Σ XΣ , in terms of which (202) readṡ This is a Ricatti equation, a particular solution of which is Then, in order to find the general solution to (203) let us introduce the variable z = y − y 0 , producinġ whose solution reads where b is an arbitrary constant of integration and δ ≡ α(1 − a 1 ) − y 0 . With this result, we can easily find X Σ , whose expression reads where c is a constant of integration.
Using (207) in (200) we obtain the explicit form of γ(t), and using this expression and (207) in (197) we obtain the explicit from of f (t). Thus, the model is completely determined up to a single function of r (h(r)).
In terms of X(t, r) and h(r), the physical variables read In order to obtain a specific model we shall assume a 1 = 2, which implies y 0 = −α and δ = 0, then feeding back these values in (207), the expression for X Σ becomes withc ≡ c (1+b) 2 . Next, we shall assume for h(r) the form where c 2 is a constant with dimensions 1/[length] 2 , producing Using (213) and (214) we obtain for the radial dependence of X rdr h(r) 2 From (200) we obtain at once for γ(t) and from (197) and (215) the expression for f (t) reads Finally the expression for X(t, r) reads Thus, the physical variables for this model XI (including the total mass and the temperature) read 8πµ = 3c 2 α 2 e −2αt 4r 4 5r 4 + 2r 2 r 2 Σ + r 4 Σ , this last expression was obtained using the truncated transport equation (39).
It is worth noticing that this model is intrinsically isotropic in pressure, the energy density is positive and larger than the pressure, and the matching condition q Σ = P r is obviously satisfied. However the physical variables are singular at the center.

VI. DISCUSSION
We have seen so far that the admittance of CKV leads to a wealth of solutions to the Einstein equations for a general spherically symmetric fluid distributions, which could be applied to a variety of astrophysical problems, or serve as testbeds for discussions about theoretical issues such as wormholes and white holes.
In order to find solutions expressed in terms of elementary functions we have imposed further constrains on the fluid distribution. Some of which are endowed with a distinct physical meaning (e.g. the vanishing complexity factor, or the quasi-homologous condition), while others have been imposed just to produce models described by elementary functions.
We started by considering non-dissipative fluids admitting a CKV orthogonal to the four-velocity. In this case the assumed symmetry reduces the metric variables to three functions (two functions of t and one function of r). Then, the matching conditions reduce to a single differential equation (65) whose solution provides one of the three functions describing the metric. In order to obtain a solution expressed in terms of elementary functions we have assumed specific values of the parameters entering into the equation.
The first choice (ω = 1 27M 2 ) leads to two expressions for the areal radius of the boundary ( (82) and (83)). The first one describes a fluid distribution whose boundary areal radius expands from 0 to 3M , while the second one describes a contraction of the boundary areal radius from infinity to 3M . To find the remaining two functions to determine the metric we have assumed the quasihomologous condition and the vanishing complexity factor condition. In this way we are lead to our models I and II, both of which have positive energy densities and the physical variables are singular free, except the model I for t = t 0 .
As t → ∞ both solutions tend to the same static solution (97) satisfying a Chaplygin-type equation of state µ = −P r . The way of reaching this static limit deserves some comments. Usually the hydrostatic equilibrium is reached when the "gravitational force term" (the first term on the right of (B6)) cancels the "hydrodynamic force term" (the second term on the right of (B6)). However here the situation is different, the equilibrium is reached because as t → ∞ both terms tend to zero.
In spite of the good behavior of these two models, it should be mentioned that regularity conditions are not satisfied by the resulting function R on the center of the distribution. Accordingly for the modeling of any specific scenario, the central region should be excluded.
Next, we have considered the case ω = 0, which together with the vanishing complexity factor condition produces the model III. In this model the boundary areal radius oscillates between 0 and 2M . The energy density and the tangential pressure of this model are positive and homogeneous, while the radial pressure vanishes identically. As in the previous two models this solution does not satisfy the regularity condition at the center.
As an additional example of analytical solution we have considered the case M = 0. The two models for this kind of solution are the models IV and V . They represent a kind of "ghost" stars, formed by a fluid distribution not producing gravitational effects outside the boundary surface. They present pathologies, both physical and topological, and therefore their physical applications are dubious. However since this kind of distributions have been considered in the past (see for example [87]) we present them here.
Next we have considered the subcase where the CKV is orthogonal to the four velocity and the fluid is dissipative. For this case we have found a model satisfying the vanishing complexity factor and the quasi-homologous condition, which together with the fulfillment of the matching conditions determine all the metric functions. This model (model V I) is described by expressions (137)-(141), and the expression (147) for the temperature, which has been calculated using the truncated version of the transport equation. It contains contribution from the transient regime (proportional to τ ) as well as from the stationary regime. As previous models, this solution does not satisfy the regularity conditions at the center.
The other family of solutions corresponds to the case when the CKV is parallel to the four-velocity. In the non-dissipative case, as consequence of this symmetry, the metric functions are determined up to three functions (two functions of r and one function of t). Besides, the fluid is necessarily shear-free, a result which was already known [88,89]. The function of t is obtained from the fulfillment of the matching conditions (Eqs. (153), (155)). These equations have been integrated for different values of the parameters entering into them. Thus, for a 1 = 1 and ǫ = 1 27M 2 , together with the vanishing complexity factor condition and h(r) = c 6 r, we have found model V II. The boundary areal radius of this model expands from zero to 3M , and the physical variables are given by (165)-(167). In the limit t → ∞ the model tends to a static sphere whose equation of state is P r = −µ. The energy density is positive, and presents a singularity only at t = t 0 , however regularity conditions are not satisfied at the center.
The integration of the matching conditions for ǫ = 0 and a 1 = 1 2 together with the vanishing complexity factor, produce the model V III. The boundary areal radius of this model oscillates between zero and 8M 3 . The energy density is positive and larger than the radial pressure, but the fluid distribution is singular at r = 0.
For M = 0 and a 1 = 1 we obtain models IX and X they describe the kind of "ghost stars" mentioned before. However they are plagued with, both, physical and topological pathologies which renders them unviable for physical modeling. We include them just for sake of completeness.
Finally, we have considered the dissipative case for the CKV parallel to the four-velocity. The metric variables for this case take the form (187)-(189), which after im-posing the vanishing complexity factor condition become (193)- (195). Thus the metric is determined up to three functions (two functions of t and one function of r). The two functions of t will be obtained from the integration of the matching conditions, while the function of r is assumed as (213). The model is further specified with the choice a 1 = 2. This produce the model XI.
As it follows from (218) the boundary areal radius of the model tends to infinity as t → ∞, while in the same limit the total mass m Σ tends to infinity, whereas both q and µ tend to zero. The explanation for this strange result comes about from the fact that R Σ grows exponentially with t, overcompensating the decreasing of µ and q in (20). It is also worth noticing the negative sign of q, implying an inward heat flux driving the expansion of the fluid distribution.
Overall, we believe that the eleven models exhibited (or at least some of them) could be useful to describe some stages of some regions of self-gravitating fluid, in the evolution of compact objects. Each specific scenario imposing specific values on the relevant parameters. It should be reminded that in any realistic collapsing scenario we do not expect the same equation of state to be valid all along the evolution and for the whole fluid configuration.
Before concluding, some general comments are in order.
1. The analytical integration of the equations derived from the matching conditions have been carried out by imposing specific values on the parameters entering into those equations, also the models have been specified by using some conditions such as the quasi-homologous condition. Of course the number of available options is huge. Among them we would like to mention the prescription of the total luminosity measured by an observer at rest at infinity (36). Let us recall that this is one of the few observables in the process of stellar evolution. Equivalently one could propose a specific evolution of the total mass with time.
2. In some cases, when the topological pathologies are not "severe", the time interval of viability of the solution may be restricted by the condition that U ≤ 1 (e.g. for solutions I and II). In other cases however, due to topological defects, the interpretation of U as a velocity becomes dubious and therefore it is not clear that U should satisfy the above mentioned condition.
3. Model XI is dissipative and intrinsically isotropic in pressure. However as shown in [71] dissipation produces pressure anisotropy, unless a highly unlikely cancellation of the four terms on the right of equation (28) in [71] occurs. This happens in model XI, which renders this solution a very remarkable one.
4. For reasons explained in the Introduction we have focused on the obtention of analytical solutions expressed through elementary functions. However it should be clear that for specific astrophysical scenarios, a numerical approach for solving the matching conditions, could be more appropriate.

ACKNOWLEDGMENTS
This work was partially supported by the Spanish Ministerio de Ciencia e Innovación under Research Projects No. FIS2015-65140-P (MINECO/FEDER). ADP ac-knowledges hospitality from the Physics Department of the Universitat de les Illes Balears.
This last equation may be further transformed as fol-lows, the acceleration D T U of an infalling particle can be obtained by using (7), (A4), (12) and (15), producing and then, substituting a from (B5) into (B4), we obtain (µ + P r ) D T U = − (µ + P r ) m R 2 + 4πP r R − E 2 D R P r + 2(P r − P ⊥ )