Einstein-aether models III: conformally static metrics, perfect fluid and scalar fields

The asymptotic properties of conformally static metrics in Einstein-aether theory with a perfect fluid source and a scalar field are analyzed. In case of perfect fluid, some relativistic solutions are recovered such as: Minkowski spacetime, the Kasner solution, a flat FLRW space and static orbits depending on the barotropic parameter $\gamma$. To analyze locally the behavior of the solutions near a sonic line $v^2=\gamma-1$, where $v$ is the tilt, a new"shock"variable is used. Two new equilibrium point on this line are found. These points do not exist in General Relativity when $1<\gamma<2 $. In the limiting case of General Relativity these points represent stiff solutions with extreme tilt. Lines of equilibrium points associated with a change of causality of the homothetic vector field are found in the limit of General Relativity. For non-homogeneous scalar field $\phi(t,x)$ with potential $V(\phi(t,x))$ the symmetry of the conformally static metric restrict the scalar fields to be considered to $ \phi(t,x)=\psi (x)-\lambda t, V(\phi(t,x))= e^{-2 t} U(\psi(x))$, $U(\psi)=U_0 e^{-\frac{2 \psi}{\lambda}}$. An exhaustive analysis (analytical or numerical) of the stability conditions is provided for some particular cases.


Introduction
According to the measurements from type Ia supernovae [64] the Universe is experiencing an accelerated expansion due to an unknown "Dark Energy" source, that was introduced in the standard cosmological model to account for a e-mail: genly.leon@ucn.cl b alfredo.millano@alumnos.ucn.cl c lattaj@mathstat.dal.ca 68% of the energy content of the universe [2]. Measurements of anisotropies of the cosmic microwave background (CMB) from experiments including the WMAP [10] and Planck [3] satellites, have provided strong support for the standard Λ CDM model of cosmology where Λ is a cosmological constant. However, there are some tensions with local measurements of the Hubble expansion rate from supernovae Ia [65] and other cosmological data [1], that settled this model under question mark. A very well-known issue of Λ CDM model is that the energy density comprised in a cosmological constant Λ has to be fine-tuned by ∼ 55 orders of magnitude to account for the present acceleration [52]. Therefore, various attempts to explain the cosmic acceleration within General Relativity (GR) were proposed as alternatives to Λ CDM, as well as several alternatives that abandon GR and modify the Einstein-Hilbert action. Within the last group, a very interesting alternative is the so-called Einstein-aether theory, which is an effective field theory in which the Hilbert action is modified by the introduction of a dynamical timelike unit vector field, u a , the aether, which is covariantly coupled, at Lagrangian level, up to the second order derivatives of the spacetime metric g ab , excluding total derivatives. The unitarity is imposed by introducing a Lagrangian multiplier in action. This theory has some features that make it of interest to mathematicians, physics, and to cosmologists. These are: a) it violates the Lorentz invariance, but preserves locality and covariance; b) it has some imprints on the inflationary scenery; c) it satisfies conditions for linearized stability, positive energy, and vanishing of preferred-frame post-Newtonian parameters; d) for generic values of the coupling constants, the aether and the metric isotropizes (although for large angles or large angle derivatives of the tilt angle there is a runaway behavior in which the anisotropies increases with time, and some singularities may appear); and e) every hypersurface-orthogonal Einstein aether solution is a Hořava solution, etc., see, e.g., arXiv:2010.03033v1 [gr-qc] 6 Oct 2020 [4-7, 11, 13, 14, 18-31, 34-48, 51, 53, 55, 59, 63, 66-70, 72, 75].
Einstein-aether theory has applications in various anisotropic and inhomogeneous contexts. In [19] it was implemented the 1+3 orthonormal frame formalism, adopting the comoving aether gauge, to obtain evolution equations in normalized variables, which are suitable for numerical calculations and for phase space analyzes. Spatially homogeneous Kantowski-Sachs models were studied. e.g., in [5,19,47,56,72]. In [5] the scalar field interacts to both the aether field expansion and shear scalars through the potential. The stability against spatially curvature and anisotropic perturbations was studied. The late-time attractor of the theory is the vacuum de-Sitter expansionary phase. Static metrics for non-tilted a perfect fluid with linear and polytropic equations of state, and with a scalar field with exponential or monomial potentials, were studied in [18,19,48]. Other solutions were examined elsewhere: vacuum Bianchi Type V [67]; Friedmann-Lemaître-Robertson-Walker metric (FLRW) [59,66]; a Locally Rotationally Symmetric (LRS) Bianchi Type III [66]; modified scalar field cosmology with interactions between the scalar field and the aether [55] based on Einstein-aether theories by [45] and [20]. An emphasis was set on the issue of the existence of solutions of the reduced equations, the classification of the singularities, and the stability analysis.
These theories are different from Scalar-tensor theories, and they are similar to the particle creation, bulk viscosity, and varying vacuum theories, or varying-mass dark matter particles theories [8,9,49,50,54,[60][61][62]. In [57] the Einstein-aether theory which incorporates a scalar field nonminimally coupled to the aether through an effective coupling B (φ ) = 6B 0 φ 2 [45] was studied. It was found there are five families of scalar field potentials on the form V A (φ ) = V 0 φ p + V 1 φ r , where p, r are specific constants, which lead to Liouville-integrable systems, and which admit conservation laws quadratic in the momenta. Following an analogous strategy in [58] were determined exact and analytic solutions of the gravitational field equations in Einstein-aether scalar model field with a Bianchi I background space with nonlinear interactions of the scalar field with the aether field. Conservation laws for the field equations for specific forms of the unknown functions such that the field equations are Liouville integrable were derived. Furthermore, the evolution of the anisotropies was studied by determining the equilibrium points and analyzing their stability.
This paper is the third of a series of works devoted to Einstein-aether theory with perfect fluids and scalar fields. In paper I [18], the field equations in the Einstein-aether theory for static spherically symmetric spacetimes and a perfect fluid source, and subsequently with the addition of a scalar field (with an exponential self-interacting potential) were investigated. Appropriate dynamical variables which facilitate the study of the equilibrium points of the resulting dynam-ical system were introduced. In addition, the dynamics at infinity was discussed. The qualitative properties of the solutions are of particular interest, as well as their asymptotic behavior and whether they admit singularities. A number of new solutions were presented. Continuing this line, in paper II [48] the existence of analytic solutions for the field equations in the Einstein-aether theory for a static spherically symmetric spacetime was investigated. A detailed dynamical system analysis of the field equations was provided.
This paper is focused on the study of timelike self-similar (TSS) spherically symmetric models with perfect fluid and/or scalar fields, using the covariant decomposition 1 + 3 [19,33,73,74]. This formalism is well-suited for performing qualitative and numerical analysis. TSS spherically symmetric models are characterized by a 4-dimensional symmetric homothetic group H 4 acting multiply transitively on 3dimensional timelike surfaces.
This paper is organized as follows. In section 2 the 1 + 3 orthonormal frame formalism is summarized. In section 3 the action for the Einstein-aether Theory is presented, following [13,30,36]. In section 4, TSS spherically symmetric models with perfect fluid are studied using the covariant decomposition 1 + 3 [19,74], following the approach given in [33]. In section 5 different stability conditions of the equilibrium solutions of dynamical systems will be established. Numerical methods will be used to support and validate the analytical results. In section 6 a non-homogeneous scalar field φ (t, x) with potencial V (φ (t, x)) which satisfies the symmetry of the conformally static metric [15] is studied. The θ -normalization procedure will be implemented in section 7. An exhaustive analysis (analytical or numerical) of the stability conditions is provided for some particular cases. Section 8 is devoted to conclusions.
The frame rotation Ω αβ is zero. The cosmological constant is chosen to be 0 for simplicity. The quatity n 13 only appears in the equations together with e 2 n 13 through the Gauss spatial curvature of the 2-spheres: 2 K := 2(e 2 − 2n 13 )n 13 , which is simplified to: Hence, the dependence on ϑ does not appears explicitly in the equations. The quantity 2 K is used instead of e 2 2 to write the field equations. The spatial curvatures simplify to: The components of the Weyl curvature are simplified to: with E + given by: Using the following simplifications, 2 K = K,u 1 =u, a 1 = a, the essential variables are: N, e 1 1 , K, H, σ + , a, µ, q 1 , p, π + , (14) and the auxiliary variables are 3 K, 3 S + ,u.
The field equations are written as: The restrictions are the Gauss and Codazzi equations together with the definition of a are the following: where the spatial curvatures are given by: Afterwards, the lapse function N is provided specifying the time gauge and, since there are not evolution equations for p and π + , they should be specified by the fluid model through equations of state for p and the transport equation for π + .

Einstein-aether Gravity
The action for the Einstein-aether Theory is the most general covariant functional involving partial derivatives of order at most two (not including total derivatives) of the space-time metric g ab and a vector field u a , called aether [13,30,36] given by: where: is the Einstein-aether lagrangian [36] with: That action contains the Einstein-Hilbert term 1 2 R, wherein: R denotes the Ricci scalar, g ab denotes the metric tensor, and K ab cd is a tensor of four indices corresponding to the kinetic terms of the aether. It contains four dimensionless constants c i and M is the Lagrange multiplier that forces unitarity of the aether vector, u c u c = −1. That is, u c is a timelike vector [30]. The signature of the metric g ab is (−+++). Physical units are such that c = 1, κ 2 ≡ 8πG = 1, where c is the speed of light.
The field equations of the Einstein-aether theory accounts for [20,41]: -The effects of anisotropy and inhomogeneities (e.g., curvature) on the geometry of the spherically symmetric models under consideration. -The contribution from the energy-momentum tensor T ae ab of the aether field, which depends on the dimensionless parameters c i , i = 1, . . . 4. In General Relativity all c i = 0, hence the Einstein's field equations are generalized.
To study the effects of matter, the values corresponding to General Relativity, or values close to them, can be substituted.
-When studying the phenomenology of theories within a preferred framework, and particularly, in the isotropic and spatially homogeneous universe, it is generally assumed the aether field will be aligned with the cosmic frame (natural resting frame preferred by the CMB) and therefore is related to the expansion rate of the universe. -In principle, in spherically symmetric models the preferred frame determined by the aether can be different (that is, tilted) to the CMB rest frame. This adds additional terms to the energy-momentum tensor of the aether T ae ab , for example, an hyperbolic angle of tilt, v, which measures the aether boost with respect the CMB rest frame [14,45]. In homogeneous but spatially anisotropic models, it is expected that the hyperbolic inclination angle v will decay along with its derivative as t → +∞ [16,17].
All spherically symmetric aether fields are hypersurface-orthogonal. Therefore, all spherically symmetric solutions of the aether theory will also be solutions in the infrared limit of gravity of Hořava. The opposite is not true in general, but it is true for solutions with spherical symmetry with regular center [36]. When spherical symmetry is imposed, the aether is hypersurface-orthogonal, and it has zero twist. Therefore, without loss of generality it is possible to make c 4 zero [36]. After redefining parameters to remove c 4 , the parameter space is 3-dimensional. The c i contributes to the effective Newtonian gravitational constant G. Then a parameter c i can be specified to make 8πG = 1. The remaining two parameters characterize two non-trivial physical quantities, for example, the Schwarzschild mass and radius of a matter distribution. The other restrictions imposed on the c i are summarized in [36] and in equations 43-46 in [6]. The field equations obtained by variation of (19) with respect to g ab , u a , and M are respectively given by [29]: where G ab is the Einstein tensor of the metric g ab , T T OT ab is the total energy-momentum tensor, T T OT ab = T ae ab + T mat ab , where T mat ab is the total contribution of matter. T mat ab will be omitted for the moment (and will be added later for models with perfect fluid and with scalar field), starting with vacuum case (L m = 0). The quantities J a b ,u a and the aether energy-momentum tensor T ae ab are given by: Taking the contraction of (22b) with u b and with the induced metric h bc := g bc + u b u c the following equations are obtained: Equation (24a) is used as the definition of the Lagrange multiplier. The second system of equations give compatibility conditions that the aether vector must satisfy.
4 Timelike self-similar spherically symmetric perfect fluid models In the diagonal homothetic formulation, the line element can be written in diagonal form, where one of the coordinates adapts to the homothetic symmetry [12]: For the conformally static metric with line element given by (25), the following scalars can be defined: where ... denotes the derivative with respect to the spatial variable x. The quantities α and β are respectively the expansion scalar and the shear scalar of the normal congruence to the symmetry surface of the static universe M , ds 2 , conformally related to the physical spacetime M , ds 2 through the homothetic factor e 2t . A non-tilted aether vector u = e −t b 1 ∂ t is considered. Assuming that the matter content of the physical universe M , ds 2 is a perfect fluid, that is specified by a 4-velocity vector field v given by where v is the tilt parameter, which is a funtion of x, with −1 ≤ v ≤ 1. The equation of state parameter p = (γ − 1)µ, 1 ≤ γ < 2 is chosen for the perfect fluid (unless otherwise stated). By convenience a function that depends only on x: is defined, which represents the energy density of the fluid measured by an observer associated with the homothetic symmetry. On the other hand, if the content of matter is that of a non-homogeneous scalar field, φ (t, x), with its self-interaction potential V (φ (t, x)), these must respect the homothecy of the conformally static symmetry associated with the line element (25), so they have to be of the form [15]: λ , where by convenience it is assumed λ > 0, such that for ψ > 0, U → 0 as λ → 0, which restrict the kind of scalar field potentials to be considered.
Using the metric (25), the lagrangian (20) becomes: The Lagrange multiplier (24a) is calculated as: The aether equation (24b) is reduced to: The trace of the intrinsic Ricci 3-curvature of the spatial 3surfaces orthogonal to u is given by Now the aether parameters are re-defined as [40]: corresponding to terms in the Lagrangian relative to expansion, shear scalar, acceleration and twist of the aether. To impose the condition of the aether (24b) is taken c a = 3c θ , (2c 1 + 3c 2 + c 3 − c 4 ) = 0. Therefore, the parameter space is reduced to a constant, c θ . The aether energy-moment tensor can be expressed by: c θ σ 2 are the effective energy density, isotropic pressure, energy flux, and anisotropic pressure of the aether, measured by an observer associated with the homothetic symmetry in the static universe M , ds 2 ; therefore, depending only on x.
The matter energy-momentum tensor is given by: Using the Einstein equations, the Jacobi identities and the contracted Bianchi identities, a system of ordinary differential equations for the frame vectors and the comutation functions, and an extra equation for the aether are obtained. The comoving gauge is chosen for aether; leaving as a degree of freedom a reparametrization of the spatial variables and the temporal variable. The final equations are: Propagation Equations: Equation for µ t : Auxiliary equation: Restriction: To simplify the notation the following parameters are introduced C 1 = 1 − 2c σ ,C 2 = 1 + 3c θ ,C 3 = 1 + c a , such that by choosing C 1 = 1,C 2 = 1,C 3 = 1, General Relativity is recovered. The condition c a = 3c θ , implies C 2 = C 3 , and the parameter C 1 does not appears explicitly in the equations, only on the definition of the Lagrange multiplier.
For an ideal gas with γ = 1 the matter energy-momentum tensor is simplified to: where v measures the inclination of the fluid relative to the 4-velocity of the aether. The equations reduce to with restrictions C 2 = 1 needs to be set to recover General Relativity . For this reason is natural to consider Due to the usual energy condition for the fluid, which is expressed as 1 < γ < 2, the second condition is satisfied. Together with the energy condition µ t ≥ 0, lead to By hypothesis C 2 > 0, θ 2 is the dominant quantity, and the terms on the left hand side of the above inequality are both non-negative. This suggests to considering θ -normalized equations.
The normalized equations of interest are: Propagation equations: Equation for Ω t : Auxiliary equation: Restriction: It is easily verified that the previous equations are invariant under the discrete transformation In particular, four specific models will be studied, these are: models with extreme tilt (58); presureless perfect fluid (72); the reduced system (77) in the invariant set A = v = 0; and the general system (75). The case v 2 = γ − 1 will be studied in section 5.1. The invariant sets v = ±1 will be analyzed in section 5.2. In section 5.3 are calculated invariant manifolds of P 6 using analytical tools. Section 5.4 is devoted to the study of the ideal gas (γ = 1). The general case v = 0 corresponding to tilted fluid will be studied in section 5.5. The equilibrium points with v = 0, A = 0 will be studied in section 5.6. In section 5.7 the sinks and sources for the model with perfect fluid are summarized. Finally, in section 5.8 results will be summarized and the relation with previous results in the literature will be discussed.
In Figure 1 the real parts of the eigenvalues of the equilibrium point M ± for C 2 = 1 are depicted, showing in general that it is a hyperbolic saddle. For γ = 6 5 , the eigenvalues of M − are {−2, 0, 0, −4}, so, it is non-hyperbolic.

Surface of non-extendibility of solutions
For γ > 1, v 2 = γ −1 represents a surface of non-extendibility of the solutions, and it is called sonic surface. The curve parametrized by Σ , is called sonic line. The solutions diverge in a finite time when the solutions approach the sonic surface v 2 = γ − 1.
The only way it can be passed through the sonic surface is when the numerator of the equation (44d) also vanishes, this is through the sonic line SL ± . In SL ± both the denominator and the numerator of (44d) are zero. This indicates the presence of a singularity of the system (44). As a difference with General Relativity, for 1 < γ < 2 and C 2 = γ 2 4(γ−1) 2 the system (75) admits the following equilibrium points: 4(γ−1) 3/2 , which lie on the sonic line. These points do not exist in General Relativity when 1 < γ < 2. When γ = 2,C 2 = 1 these points exist, and since γ = 2 the fluid behaves like stiff matter. Additionally, if γ = 2,C 2 = 1, these points correspond to models with extreme tilt (v = ε), SL 1 : On the sonic surface the inequality C 2 (A 2 + Σ 2 ) ≤ 1 must be satisfied, which corresponds to K ≥ 0, which imposes additional conditions on the parameters. To analyze locally the behavior of the solutions near this sonic line, the new "shock" variable ξ is introduced: that leads to the system The new variable is not monotonic since (γ − 1) − v 2 can change the sign, so the system is not suitable to do qualitative/asymptotic analysis of the system since it does not represent a dynamical system. However, the system is suitable for numerical integration in a neighborhood of the sonic line SL ± , and for the local stability analysis of SL ± at the perturbative level. SL ± can be parametrized by: Defining the equation of state parameter ω = γ − 1, we deduce that ω = v 2 0 . Defining the following linear perturbations: the evolution equations of the perturbations are given by: The eigenvalues of the equilibrium point (δ Σ , δ A , δ v ) = (0, 0, 0) are λ 1 = 0 and the roots λ 2 y λ 3 of the polynomial: A zero eigenvalue appears because it is a curve of equilibrium points. Also, the curve: The eigenvector associated with the zero eigenvalue: is parallel to the vector tangent to the curve at the point. Then the curve is normally-hyperbolic, so the stability can be studied considering only the signs of the nonzero eigenvalues.
For the stability analysis of the sonic line SL ± , the following invariant sets are identified: both invariant sets determine curves in the space of parameters v 0 , Σ 0 , which, due to the fact that they are invariant cannot be passed by orbits. For the discussion the following existence conditions will be imposed: In Figure 2 stability regions for are depicted in the space of parameters, for (a) Σ in a fixed point at the sonic curve. The unshaded region represents the region where A 0 < 0 or K 0 < 0, which is the non-physical region. The dotted red line corresponds to K = 0 and the thick blue line corresponds to A = 0. The region represented in gray color corresponds to the region where (δ Σ , δ A , δ v ) = (0, 0, 0) is a hyperbolic saddle. The region represented in black color corresponds to the region where (δ Σ , δ A , δ v ) = (0, 0, 0) is stable. Figure 2(b) reproduces the stability results shown in Figure 1 of [33] (with the exception of a strip in the parameter space, represented in gray, where the point equilibrium (δ Σ , δ A , δ v ) = (0, 0, 0) is a hyperbolic saddle, see upper right corner in Figure 2(b), whose analysis was omitted in [33]).
For SL + the existence condition is not verified, so they are outside the physical region and their analysis is omitted.

Invariant sets v = ±1
In this section, the following invariant sets are studied v = ±1, which corresponds to extreme tilt. Then, from the equations (47), (45), and (43) it follow: Afterwards, the following reduced 2-dimensional system is obtained: where ±1 denotes the sign of v.
The systems (58) are related through the simultaneous change of A → −A, and the sign " + " by the sign " − ". Therefore, without loss of generality, the positive sign " + " is studied, with A ≥ 0. Table 1 presents the qualitative analysis of the equilibrium points of the systems (58) corresponding to the cases of extreme tilt v = 1, −1, which are the following: , 2 is a hyperbolic source for C 2 > 0. Figure 3 shows some orbits of the system (58) with v = 1, −1, for different choices of parameters. The points P 5 , P 6 , P 7 are included, where P i denotes the symmetric points of the points P i . The dotted invariant line represents H − (resp. H + ) for C 2 = 1. The analytical results discussed above are confirmed. Figure 4 shows some orbits of the system (60) for different choices of parameters. Table 1 Qualitative analysis of the equilibrium points of the systems (58) corresponding to extreme tilt: v = 1, −1.

Invariant manifolds of P 6
Concerning the equilibrium point corresponds to its unstable manifold for C 2 < 1. This line is stable for C 2 = 1, and for C 2 > 1 that line belongs to the stable 2-dimensional manifold of P 6 . Introducing the change of variables: the following equations are obtained: where the following parameter is introduced: The eigenvalues of P 6 are − 2µ µ+1 , −2 . For µ > 0, i.e., C 2 > 1, the equilibrium point P 6 is a hyperbolic sink in the invariant set v = −1. For µ < 0, i.e., C 2 < 1, the equilibrium point P 6 has an unstable manifold tangent to the axis x. This unstable manifold of P 6 can be expressed locally by the graph: where h satisfies the initial value problem: The previous differential equation admits the first integral where C 1 is an integration constant. Imposing the condition h(0) = 0, it is obtained C 1 = 0. Solving the resulting equation for h are obtained two solutions: The last solution is discarded since h (x) = −1, which implies that the tangential condition h (0) = 0 is not fulfilled. Then, the unstable solution is given locally by the trivial solution: The dynamics on the invariant manifold is given by: whose solution passing by x(0) = x 0 , is: That is, the solutions generically approaches the origin as η → −∞, tending towards the past to P 6 . Figure 5 shows the 1-dimensional flow of (67) for −1 < µ < 0, where it is illustrated that the origin of the system (67) is stable. Then, applying the Unstable Manifold Theorem, it is confirmed that P 6 is a hyperbolic saddle, as it was previously commented in table 1.   (67), for µ = −0.7. All the flows are topologically equivalent for µ on the interval −1 < µ < 0.

Ideal gas γ = 1
For an ideal gas with γ = 1 (pressureless fluid; dust) the equations (38) and restrictions (39) become: where the expressions Ω t and r given by: were used. The restriction (47) is reduced to: This allows us to study the reduced 3-dimensional system: defined on the phase space: The equilibrium points of system (72) are the following: , 2 are hyperbolic saddles for C 2 > 0.
3/2 and λ 1 , λ 2 and λ 3 are the roots of the polynomial in λ : λ 3 . Figure 6 shows the real part λ i (which differs from the eigenvalues in an overall multiplicative factor k) corresponding to the equilibrium point P 8 for C 2 ≥ 3 4 and γ = 1, showing that in general it has saddle behavior.
These lines of equilibrium points, which do not exist in Einstein-aether theory (C 2 = 1), are associated with a change of causality of the homothetic vector field.
The stability criteria of the equilibrium points of the system (72) for pressureless fluid (γ = 1) and v = 0 are summarized in table 2.

General case
In the general case v = 0 is possible reduce the system's dimension when the restriction (47) in non-degenerated. The restrictions (47) and (45) can be globally solved for Ω t y K (assuming γv = 0 and γv 2 − v 2 + 1 = 0): Then, a 3-dimensional dynamical system is obtained: (75c) The dynamical system (75) admits some invariant sets. These are: v = ±1, corresponding to extreme tilt, and the invariant sets A = 0 and Σ = 0. The equilibrium points of the system (75) are the following.
(b) non-hyperbolic with three zero eigenvalues for C 2 = 1. (c) non-hyperbolic with a zero eigenvalue and two purely imaginary eigenvalues for 0 ≤ C 2 < 1.
This shows a hyperbolic saddle behavior for values of the parameters close to the values of General Relativity. For example, for a pressureless fluid(γ = 1) the λ i are approximately −2δ , δ + √ δ , δ − √ δ . For δ < 0 there are two complex imaginary eigenvalues with negative real parts and a positive real eigenvalue, while for δ > 0 there is a negative eigenvalue and the others have different signs. In the figure 7 the real part of λ i (which differ from the eigenvalues associated with the equilibrium point P 8 by the overall factor k) are represented graphically for C 2 ≥ 2+γ 4γ and γ ∈ [1, 2]. The figure illustrates that the equilibrium point is generally a hyperbolic saddle or is non-hyperbolic.
, and Noting that there is at least one change of sign in two eigenvalues, that is µ 1 µ 2 < 0, for 1 < γ < 2,C 2 > 0, it is concluded that it is a hyperbolic saddle.

Invariant set
The focus of this section is the stability analysis of the equilibrium points of the system (49) in the invariant set A = v = 0. In the following list the stability analysis is done with preserving the four eigenvalues. , 0, 0, 0 . The eigenvalues are (a) It is a hyperbolic sink for γ > 1, It is a saddle for 1 < γ < 2,C 2 > 0.
It is: (a) a hyperbolic source for 1 < γ < 2, 0 < C 2 < (2−3γ) 2 16(γ−1) 2 . (b) a hyperbolic saddle for: It is: (a) non-hyperbolic for: Figure 8 shows the real parts of the eigenvalues λ i for the equilibrium point , 0 . It represents static solutions for 1 ≤ γ ≤ 2 y C 2 ≥ 0. The figure shows that the equilibrium point is non-hyperbolic in the cases (a)-i,ii,iii previously described, or, it is a hyperbolic saddle. (47) is trivially satisfied. On the other hand, from equation (45) it follows

Reduced system
Imposing the energy condition Ω t ≥ 0 and choosing γ ∈ (1, 2], the following reduced dynamical system is obtained: defined in the phase space: The qualitative analysis of system (77) is given in Table 3.

Summary of sources and sinks for a perfect fluid
A summary of the equilibrium points classified as sinks or sources of the model with perfect fluid is presented.

Non-extensibility surface of solutions
SL − is stable for the regions described in Figure 2.

Discussion
In this section, timelike self-similar spherically symmetric metrics in Einstein-aether theory with perfect fluid as matter content were studied using the homothetic diagonal formulation, which gives the propagation equations (33) plus algebraic restrictions. The homothetic diagonal formalism has some disadvantages. Symmetry surfaces generally change causality. Then, in the homothetic diagonal formulation spacetime must be covered with two coordinate systems (two charts); one when homothetic Killing vector is timelike, and another when homothetic Killing vector is spacelike. These two regions have to be matched in the region where the Killing vector is null [33]. However, the formulation has more advantages than disadvantages. The main one is that it allows the field equations, which are a well-defined system of firstorder partial derivative equations (PDE), in two variables (from the 1 + 3 formalism), to be written as a system of ordinary differential equations using the symmetries that come from the Killing vectors. The resulting equations are very similar to those of the models with homogeneous hypersurfaces. In turn, it is possible to write these equations as a dynamical system, which makes it possible to study the model using the techniques of the qualitative theory of dynamical systems. This makes it possible to obtain a complete description in a phase space, which leads to a better understanding of the dynamics of the model. In this regard, the θ -normalized equations were presented.
. hyperbolic saddle (see text).  Table 3 summarizes the results of the qualitative analysis of the equilibrium points of the system (49) with v = A = 0. Line N 1 is included. The figure 9 shows some orbits in the phase space of the system (77) with v = A = 0 and 1 < γ ≤ 2.

θ -Normalized equations
The following θ -normalized variables are introduced: along with the radial coordinate: The parameter r is defined, analogously to the "Hubble gradient parameter" r, by In these variables, the dynamical system is given by: The restrictions are: These can be globally resolved for Ω t and K to get: Finally, the following reduced 5-dimensional system is obtained: If u = w = 0 and the limit λ → 0 is taken, system (75) is recovered. Given the computational difficulty of obtaining (analytically) the stability conditions for all the equilibrium points of the system (90), in the following sections some subcases of (90) of special interest will be studied: perfect fluid in the form of ideal gas (94) , the invariant set Σ = 0 (96), the extreme tilt case (97) and the invariant set A = v = 0 (98). An exhaustive analysis (analytical or numerical) of the stability conditions is provided for these particular cases. Relaxing the condition 1 ≤ γ ≤ 2, interestingly, a cosmological fluid in the form of an ideal gas with equation of state p m = (γ − 1)µ m , with γ = 2/3 describes a FLRW spacetime with non-zero curvature.
The equations are: with restrictions Afterwards, the restrictions are solved to find: Therefore, the following reduced 3-dimensional dynamical system is obtained : The equilibrium points of the system (94) are the following.
It is: (a) a hyperbolic source for: where the expressions for µ i depend on λ ,C 2 . They exist for C 2 ≥ λ 2 2 . In Figure 10, the real parts of the µ i 's are depicted, so, the equilibrium points are typically saddles or non-hyperbolic with three zero eigenvalues.

Invariant sets v = ±1
Assuming v = ε = ±1, the system is reduced to the following 4-dimensional dynamical system: The equilibrium points of the systems (97) for ε = ±1 are the following: librium points contains the equilibrium points P 1 , P 2 , P 3 , P 4 studied in section 4. It exists for λ > 0, 0 < C 2 ≤ 1 Σ 2 y Σ ∈ R. This line is a normally-hyperbolic invariant set. Indeed, the parametric curve can be expressed as: Its tangent vector evaluated at Σ 0 : is parallel to the eigenvector corresponding to the zero eigevalue: In this particular case, the stability can be studied considering only the signs of the real parts of the non-zero eigenvalues. In this way it is concluded that: (a) Q 15 is a hyperbolic source for: 2 , 0 where one eigenvalue is zero and the other three eigenvalues are the roots of the polynomial in µ: It is either a non-hyperbolic saddle or a center depending on the choice of parameters. Figure 11 graphically represents the real part of the µ i corresponding to the equilibrium point Q 17,18 , illustrating that the equilibrium points have saddle behavior or they are nonhyperbolic.
. It exists for: , for different choices of the parameter λ . This figure illustrates that the point has a general saddle behavior.
The real parts of f i 's are represented in figure 14 (left panel).
with ε = −1 exists for λ > 0,C 2 ≥ λ 2 2 with eigenvalues The real parts of the g i 's are represented in Figure 14 (right panel). Summarizing, according to Figure 14 these equilibrium points, Q 30 and Q 31 , are either sources or they are non-hyperbolic (with four pure imaginary eigenvalues).

Invariant set
At the invariant set A = v = 0, the equations are reduced to: with restriction For 1 < γ ≤ 2, we have the auxiliary equation:

Reduced system
For 1 < γ ≤ 2, the restriction (98e) can be globally solved for Ω t , leading to the reduced system: For u = w = 0 and λ → 0 the system (77) is recovered. Therefore, the equilibrium points, and their stability conditions studied in section 5.6.1 are retrieved. By definition K ≥ 0, and w ≥ 0 (if θ > 0). Given that the system (99) is invariant to the simultaneous change (u, λ ) → (−u, −λ ), it can be assumed λ > 0. In the following discussion the analysis is restricted to λ > 0, u ≥ 0 (the sign of u corresponds to the sign of Ψ , if θ > 0). The following lists contains the equilibrium points of the reduced system (99).    Fig. 13 Real parts of the eigenvalues of the equilibrium point (Σ , A, u, w) = − λ 2 2C 2 +λ 2 , 0, , for different choices λ .

Fig. 14
Real parts of f i (left panel) and g i (right panel) corresponding to the equilibrium points Q 30 : , ε = −1, respectively. Then, Q 30 and Q 31 , are either sources or they are non-hyperbolic (with four pure imaginary eigenvalues).
line is normally-hyperbolic, i. e., given the curve parametrization: the tangent vector at Σ 0 : is parallel to the eigenvector associated to the zero eigen- In this particular case, the stability of the curve of equilibrium point can be studied considering only the signs of the real part of the non-zero eigenvalues, concluding that: (a) L is a hyperbolic source for: i.
The following list contains new points that were not studied in section 5.6.1.
The real parts of the µ i 's are represented in the figure 15 for some values of C 2 , where it is shown that P 14 (λ ) is typically a hyperbolic saddle for the given values of the parameter C 2 (or is it nonhyperbolic). The following points are recovered: lim λ →0 P 13 (λ ) = P 13 and lim λ →0 P 14 (λ ) = P 14 .
It is assumed that λ > 0, such that for ψ > 0, U → 0 as λ → 0. The equations were normalized with the variable θ . Due to the computational complexity of the resulting problem, it was not possible to obtain and analytically treat all the equilibrium points of the system (90). Hence, only some special cases of physical interest were considered, being (94) corresponding to a perfect fluid in the form of an ideal gas, (96) corresponding to the solutions in the invariant set Σ = 0, (97) corresponding to the case of extreme inclination. Finally, the invariant set A = v = 0 sets of system (98) was studied. The hyperbolic points were completely classified according to their stability conditions.

Conclusions
In this paper the space of the solutions of the differential equations that result from considering perfect fluid and/ or scalar field as the matter content in the Einstein-aether theory was studied. Einstein-aether theory of gravity consists of General Relativity coupled to a vector field of unit time type, called the aether. In this effective theory, the Lorentz invariance is violated, but the locality and the covariance are preserved in the presence of the vector field.
In section 2 the 1 + 3 formalism was discussed. This formalism is useful to write the field equations as a system of partial differential equations in two variables for spherically symmetric metrics. Furthermore, using the homothetic diagonal formulation, it was possible to write the partial differential equations as ordinary differential equations using the fact that the metric adapts to homothetic symmetry. The resulting equations (with algebraic restrictions) are very similar to those of the models with homogeneous spatial hypersurfaces [32]. It was then possible to use the techniques of the qualitative theory of dynamical systems for the stability analysis of the solutions of the models. The analytical results were verified by numerical integration.
In the section 3 the Einstein-aether theory of gravity was presented, which contains the theory of General Relativity as a limit. Conformally static metrics were studied in Einstein-aether theory for models of physical interest, such as pressure-free perfect fluids, perfect fluids, and models with extreme tilt. The stability criteria of the equilibrium points of the dynamical systems were obtained and discussed, imposing restrictions on the parameter space. Phase portraits were also presented to illustrate the qualitative behavior of the solutions. The equilibrium points obtained by [33] are recovered as particular cases of the present model. In the notation Kernel sgn(v) sgn(Σ ) the kernel indicates the interpretation of the point: M,C represent the Minkowski spacetime; K represents a Kasner solution; T corresponds to static solutions; SL ± corresponds to a flat FLRW space and static orbits depending on the parameter γ. H is associated with a change of causality of the homothetic vector field. The following results were retrieved: SL ± : Sonic lines defined by A = − γε(γ(Σ +2)−2) 4(γ−1) 3/2 , v = ε √ γ − 1, were analyzed in section 5.1. As a difference with General Relativity, for 1 < γ < 2 and C 2 = γ 2 4(γ−1) 2 the system (75) admits the following equilibrium points: 4(γ−1) 3/2 , which lie on the sonic line. If γ = 2,C 2 = 1 these points exist, and since γ = 2 the fluid behaves like stiff matter. Additionally, if γ = 2,C 2 = 1, these points correspond to models with extreme tilt (v = ε), SL 1 : Σ = 1, A = −2, v = 1, and SL 2 : Σ = −1, A = 0, v = −1. SL ± corresponds to a flat FLRW space and static orbits depending on the parameter γ. T : (Σ , A, v) = −2 γ−1 3γ−2 , 0, 0 , (K, Ω t ) = γ 2 +4(γ−1) (3γ−2 , 4(γ−1) (3γ−2 , corresponds to P 13 for C 2 = 1. These results are of interest in Cosmology and Astrophysics. In section 6 conformally static metrics were studied in Einstein-aether theory for models with tilted perfect fluid and inhomogeneous scalar field with exponential potential, so the model contains the model studied in section 4 and thus contains the model studied in [32]. Particular cases of interest in Physics were studied, such as the perfect fluid in the form of an ideal gas, solutions with Σ = 0, models with extreme tilt and the invariant set A = v = 0. It was possible to study a more general model than the one studied in [32], and the results obtained by the authors were reproduced through the use of techniques from the qualitative theory of dynamical systems. A qualitative analysis of some invariant points was also made for models with timelike self-similar spherically symmetric metrics with a perfect fluid and a scalar field. New equilibrium points were obtained, and their stability conditions were found either numerically or analytically, by imposing restrictions on the parameter space.