Analytic solutions in Einstein-aether scalar field cosmology

In the context of Einstein-aether scalar field cosmology we solve the field equations and determine exact and analytic solutions. In particular, we consider a model proposed by Kanno and Soda where the aether and the scalar fields interact through the aether coefficient parameters, which are promoted to be functions of the scalar field. For this model, we write the field equations by using the minisuperspace approach and we determine the scalar field potentials which leads to Liouville--integrable systems. We solve the field equations for five families of scalar field potentials and, whether it is feasible, we write down the analytic solutions by using closed-form functions.

In [39] it was explored the impact of Lorentz violation on the inflationary scenario. More precisely, it is studied homogeneous but anisotropic solutions in the presence of a positive cosmological constant, with a Bianchi type I (Kasner-like) symmetry with three orthogonal principal directions of expansion, and with the aether tilted in one of the principal directions. In this model the inflationary stage is divided into two parts; the Lorentz violating stage and the standard slow-roll stage. In the first stage the universe expands as an exact de Sitter spacetime, although the inflaton field is rolling down the potential. Interestingly, exact Lorentz violating inflationary solutions can be found in the absence of an inflaton potential. To linear order in the anisotropy, the system relaxes exponentially to the isotropic, de Sitter solution. This approach was an special case of the perturbative treatments used in [42]. In [42], it was investigated large deviations from isotropy, maintaining homogeneity. It was found that, for generic values of the coupling constants, the aether and metric isotropizes if the initial aether hyperbolic boost angle and its time derivative in units of the cosmological constant are less than order O(1). For larger angles or larger angle derivatives, the behavior is strongly dependent on the values of the coupling constants. In general, there is a runaway behavior in which the anisotropy increases with time, and/or singularities occur. In [43] it was studied the Einsteinaether theory with an scalar inflaton coupled bilinearly to the expansion of the aether. There were determined the conditions for linearized stability, positive energy, and vanishing of preferred-frame post-Newtonian parameters, and examined whether all of these restrictions can be simultaneously satisfied. In a homogeneous and isotropic cosmology, the inflaton-aether expansion coupling leads to a driving force on the inflaton that is proportional to the Hubble parameter. This force affects the slow-roll dynamics, but still allows a graceful exit of inflation.
Einstein-aether theory have been applied also in various anisotropic and inhomogeneous models with many interesting results. In [49] were studied spherically symmetric cosmological models in Einstein-aether theory with a non-comoving perfect fluid source using a 1+3 frame formalism, in the context of inhomogeneous cosmological models. Adopting the comoving aether gauge it is derived the evolution equations in normalized variables to provide numerical computations and studying the local stability of the equilibrium points of the resulting dynamical system. Special emphasis was made on spatially homogeneous Kantowski-Sachs models, see also [50][51][52]. In [51] was studied the dynamics of spatially homogeneous (SH) Einstein-aether cosmological models with an scalar field with a self-interaction generalized harmonic potential, in which the scalar field is coupled to both the aether field expansion and shear scalars. The stability analysis indicated that there exists a range of values of the parameters where the late-time attractor corresponds to an accelerated expansion phase. For the analysis are considered spatially curvature and anisotropic perturbations. On the other hand, static anisotropic models for a mixture of a necessarily non-tilted perfect fluid with a barotropic equation of state (linear and/or polytropic equations of state) and a self-interacting scalar field were studied in [49,53,54]. In [55] it was presented the solution space of the field equations in the Einstein-aether theory for the case of a vacuum Bianchi Type V space-time. In this model the reduced equations not always admits a solution. Whenever a solution do exist, their physical interpretation was examined through the analysis of the behavior of Ricci and/or Kretschmann scalar, as well as with the identification of the effective energy momentum tensor in terms of a perfect fluid. There are cases in which no singularities appears and in other cases the effective fluid is isotropic. Friedmann-Lemaître-Robertson-Walker metric (FLRW) and a Locally Rotationally Symmetric (LRS) Bianchi Type III space-time were studied in [56]. It was examined whether the reduced equations do have a solution, and it was found that there are portions of the initial parameters space for which no solution is admitted by the reduced equations.
In [57] it is considered an Einstein-aether scalar field cosmological model where the aether and the scalar field are interacting through two different interactions proposed in the literature by Kanno and Soda [39] and by Donnelly et al. [43]. It was provided an extended dynamical systems analysis of the cosmological evolution. The reduced Lagrangians deduced from the full action are, in general, correctly describing the dynamics whenever solutions do exist. Furthermore, the cosmological evolution of the field equations in the context of Einstein-aether cosmology by including a scalar field in a spatially flat FLRW spacetime was studied in [58] by using dynamical system tools. The analysis was separated into two cases: a pressureless fluid source is included or it is absent. The limit of general relativity is fully recovered, while the dynamical system admits de Sitter solutions which can describe the past inflationary era and the future late-time attractor. Results for generic scalar field potentials were presented, while some numerical behaviors were given for specific potential forms.
The plan of the paper is as follows.
In Section 2, we present the cosmological model under consideration which is that of the Einstein-aether gravity with an scalar field coupled to the aether through an effective coupling B (φ) as defined by B (φ) = β 1 (φ)+3β 2 (φ)+β 3 (φ)−1 in terms of the aether parameters β 1 , . . . β 4 [39]. As far as for the physical space is concerned, we consider that it is described by the spatially flat FLRW metric. For the latter cosmological model we present the field equations and we give the minisuperspace description of the theory as well.
In Section 3, we determine exact solutions of the field equations of physical interest. Specifically we find the scalar field potential such that the scale factor of the FLRW metric describes either the de Sitter universe, or it describes an scaling solution. In addition, we study the stability of these solutions by calculating their first order perturbations around the exact solutions and analyzing their evolution.
The main results of our analysis are presented in Section 4. We assume the presence of a dust fluid in the cosmological model. We determine the functional forms of the scalar field potential such that the field equations are Liouville-integrable, with at least the existence of a second conservation law, quadratic in the momentum. We find five families of power-law potentials, for which we present the analytic solutions as functions in a close-form or in algebraic form by solving the Hamilton-Jacobi equations and reducing the dimensionality of the field equations. By using the results of Section 3, we infer the asymptotic behavior of the cosmological solutions, since we can relate the dominant terms of the scalar field potential with the exact solutions presented in Section 3. Recall that a system of polynomial differential equations is said to be Liouville-integrable, if it has first integrals given by elementary functions or integrals of elementary functions, that is, functions expressed in terms of combinations of exponential functions, trigonometric functions, logarithmic functions or polynomial functions (see, e.g., [59], in the context of Tolman-Oppenheimer-Volkoff approach for a relativistic star model with the isothermal equation of state p m = ρ m /n; which is Liouville-integrable if and only if n ∈ {−1, −3, −5, −6}).
In Appendix A, we present the five Liouville integrable scalar field potentials where the additional matter source in the cosmological fluid it is an ideal gas with equation of state p m = (γ − 1) ρ m . Note that when γ = 2 3 , our results describe the case of a non spatially flat FLRW spacetime.
Finally, in Section 5, we summarize the results and we draw our conclusions.

EINSTEIN-AETHER SCALAR FIELD COSMOLOGY
We consider the Einstein-aether scalar field theory with Action Integral [39]: where S Aether describes the terms of the aether field u µ as follows Coefficients β 1 , β 2 , β 3 and β 4 define the coupling between the aether field and the gravitational field. In Einsteinaether theory, the coefficients are constants, though in this model, coefficients β 1 , β 2 , β 3 and β 4 define a coupling between the aether field u µ and the scalar field φ (x µ ), by promoting themselves to be functions of φ. Additionally, function λ is a Lagrange multiplier which ensures the unitarity, u µ u µ + 1 = 0, of the aether field u µ . In large scales the universe it is assumed to be isotropic and homogeneous described by the spatially flat FLRW metric, with line element where a (t) is the scale factor, N (t) is the lapse function while the Hubble function is defined as H (t) = 1 Nȧ a , where a dot denotes total derivative with respect the variable t.
For the aether field u µ = 1 N δ µ t , and the line element (3), the Action Integral (1) is simplified as follows [39]: where L N, a,ȧ, φ,φ is the point-like Lagrangian while function B (φ) is defined as B (φ) = β 1 (φ) + 3β 2 (φ) + β 3 (φ) − 1, and we have assumed that the scalar field φ inherits the symmetries of the spacetime such that φ = φ (t). Variation with respect the variables a and φ of the Action Integral (4) gives the second-order field equations Equation (7) is the modified Klein-Gordon equation for the scalar field φ, while equation (6) is the modified second Friedmann equation. Moreover, variation of (4) with respect to the variable N produces the modified first Friedmann equation, that is, the constraint equation, The field equations (6), (7) and (8) can be written as follows andφ where k eff = 1 B(φ) , and ρ eff and p eff describe the energy density and the pressure of the effective fluid, defined as We observe that there are similarities with the Scalar-tensor theories as mentioned in [39], indeed k eff is not a constant but changes in time, however the two theories are different. The effective k eff is the only contribution of the aether field in the first Friedmann equation, because the effective energy density ρ eff is that of the scalar field, i.e. ρ eff = ρ φ . On the other, hand from second Friedmann equation we see that the 2B ,φ Hφ modifies the effective pressure from that of the scalar field, that is, p eff = p φ + 2B ,φ Hφ. Consequently, the parameter for the effective equation of state it is defined as As far as equation (11) is concerned, this readṡ which looks like the particle creation, bulk viscosity or varying vacuum theories [60][61][62][63][64][65][66]. Positive values of B φ indicate particle annihilation while negative values of B φ indicate particle creation.
In the presence of an additional fluid source, such that of an ideal gas p m = (γ − 1) ρ m which we assume that it is not interacting with the scalar field or with the aether field, the effective energy density and pressure terms are modified as with the additional conservation equationρ from which we infer ρ m = ρ m0 a −3γ , ρ m0 is an integration constant. In the following Section, we assume that the additional matter source is that of a dust fluid, that is, γ = 1 and ρ m = ρ m0 a −3 .

EXACT SOLUTIONS
In this section we present some exact solutions of the field equations. In particular we determine the functional forms of the potential V (φ) and the function φ (t) by incorporating the requirement that the de Sitter solution a (t) = a 0 e H C t and the scaling solution a (t) = a 0 t p , are special solutions of the field equations.
Recall that we have assumed that N (t) = 1. In addition we assume that there is not any contribution in the cosmological fluid by the ideal gas, i.e. ρ m0 = 0. Because there are only two independent equations and there are three unknown functions, namely, φ (t) , V (t) and B (t), we proceed further by defining the exact form of B (φ (t)). In particular we select B (φ (t)) = 6B 0 φ 2 .

De Sitter solution
The exponential scale factor a (t) = a 0 e H0t solves the field equations (6)- (8) if and only if [39] that is In order to determine the stability of the de Sitter solution we substitute B0H0t + εδφ (t) in the field equations and we linearize around ε → 0. We end with the linearized system with closed-form solution Therefore we conclude that the expanding de Sitter universe is stable when 0 < B 0 < 1 24 , while when H 0 < 0 the exact solution is stable when B 0 > 1 24 .

Scaling solution
In a similar way, we find that the scaling solution a (t) = a 0 t p √ B 0 satisfies the field equations (6)-(8) when [39] that is We take linear perturbations around the exact solution as before and for the perturbations we find δa t R , δφ = t S , where from which we infer that the scaling solution is attractor when 0 . We proceed with the presentation of the analytic solutions.

ANALYTIC SOLUTIONS
For the lapse function N = 1, and when dust fluid is included in the model, the point-like Lagrangian (5) is written as which describes the motion of a particle in a two-dimensional space, where now the constraint equation (9) correspond to the Hamiltonian conservation law for Lagrangian (29) with value the ρ m0 . The equation of motions depend on two unknown functions, the B (φ) and the V (φ). Function V (φ) is a potential term, while function B (φ) defines the geometry of the two-dimensional space where the motion of the point-like particle occurs.
The authors of [39] considered the function B (φ) in the particular form B (φ) = 6B 0 φ 2 , and, in our work, this specific function will be selected as well. The reason is that B (φ) = B 0 φ 2 simplifies the dynamics such that the minisuperspace defined by the kinetic part of Lagrangian (29), will be a maximally symmetric two-dimensional space with zero curvature, that is, a two-dimensional flat space. Therefore, the field equations describes a typical dynamical system of Classical Mechanics.
As we commented before, we follow [39] and we set V (φ) = V 0 φ 2 . For that specific functional forms of B (φ) and V (φ), the field equations are written We define the canonical variables such that equations (30), (31), (32) are written in simpler expressions as follows x from which we derive the analytic solution in closed-form functions with constraint ρ m0 = −V 0 x 0 y 0 cosh (x 1 − y 1 ), and ω 2 = V 0 1 8B0 − 2 . Because V 0 is positive, we conclude that when B 0 > 1 16 we have a bounced universe, while when 0 < B < 1 16 the scale factor a (t) is described by hyperbolic functions.
Consider now the special case where x 1 = y 1 = 0, then we find the scale factor From the latter scale factor we see that the Λ-cosmology is not recovered since which is the solution of General Relativity of an ideal gas with a cosmological constant term. Recall, that the field φ contributes in the cosmological fluid while function B (φ) affect the total fluid source. When B 0 = 1 16 , we introduce the new canonical variables where the field equations take the form When the scalar field is massless, i.e. V (φ) = 0, then field equations are reduced, and the generic solution can be easily constructed by equations (35), (36) for V 0 = 0 and the transformation rule (33).
These are not the only functional forms of the potential V (φ) for which we can write the analytic solutions of the field equations. Some power-law functions V (φ) and their analytic solutions are presented in what it follows. Specifically, the potentials for which we shall present the analytic solution of the field equations are Potentials (43)-(47) have a common feature, they are Liouville-integrable, for which the field equations admit an additional conservation law for each potential, more specifically a quadratic conservation law, different for each potential. The V (φ) = V 0 φ 2 is also a super-integrable potential. Another superintegrable potentials we observe arē Integrable cosmological models in modified theories of gravity have been widely studied in the literature, and have been drawn the attention of cosmologists and of the mathematicians ever. The main reason is that analytic solutions can be used as toy models in order to understand the main properties of a proposed cosmological theory [67][68][69][70][71][72][73][74][75][76][77].
In order to describe Nature, we need a large number of free parameters or boundary conditions, which makes numerical treatment worthless. Because in general the equations which describe a specific theory are nonlinear such that numerical solutions may be sensitive on small changes of the initial conditions. Consequently, we refer to analytic techniques in order to understand the generic properties of a propose physical theory. Hence, the knowledge for the existence and the determination of analytical or exact solutions for a given dynamical system is important for the detailed study and understanding of the given physical theory.
According to the results of the previous Section, we observe that according to which term of the potential dominates, the behavior of the scale factor will be near to that of the scaling solution or to that of an exponential solution. For instance consider the potential V B (φ). For large values of φ, it follows V B (φ) φ 2 from which we infer that solution approaches the de Sitter universe, on the other hand, for small values of φ, V B (φ) φ − 1 2 √ B 0 , from which we infer that the scale factor behaves like that of the scaling solution. As far as potentials V D (φ) , V E (φ), are concerned, we remark that they have two terms, where only scaling solutions are described.
The method that we apply in order to determine the analytic solutions is based on canonical coordinates, as described in the example V (φ) = V 0 φ 2 . In the following, we present the analytic solutions.
In the special case where y 2 = 0 exact solution is simplified For the latter two solutions, namely (48), (49) and (51), (52) √ B0 < 0, for large time the dominated term is e ωt , which means that the scale factor for large values of t, approaches that of the de Sitter universe a (t) = a 0 e H0t .
For the latter exact solution (51), (52), in Fig. 1 we present the qualitative behaviour of the Hubble function H (z) and of the parameter for the equation of state for the effective fluid w (z) in terms of the redshift 1 + z = 1 a . Now, if we assume that V 0 = 0, then for the super-integrable potentialV A (φ) we find the exact solution with constraint equation ρ m0 − 8B0 1−16B0 y 1 x 1 = 0. For solution (53), (54) for large values of time the scale factor has a power-law behaviour a (t) = a 0 t p , where p = p (B 0 ).
which was the one studied before. Hence, we continue by presenting the analytic solution for potential V B (φ).
If x 2 = 0, the closed-form solution is y (t) = y 1 e ωt + y 2 e −ωt + 4 from which we infer, similarly as for the potential V A (φ), that the scale factor for large values of t, it is approximated by that of the de Sitter universe. When V 0 = 0, the analytic solution is found to be with constraint condition ρ m0 − 8B0 1−16B0 y 1 x 1 = 0, while when x 1 = 0 the generic analytic solution is The latter solutions are physically accepted if and only if ρ m0 = 0, that is, there is not any contribution by the dust fluid in the cosmological fluid. Consequently from (33) we infer that the scale factor for the solutions with V 0 = 0 have a power-law expression.
For B = 1 16 , we work with the variables {u, v} which are defined by expression (40). Hence, the field equations are reduced to the following systemü − 8V 1 e 3v = 0, from where it follows the generic solutions where D + (t) is the Dawson function defined as D + (t) = e −t 2 t 0 e r 2 dr, while from the constraint equation it follows From solution (68), (69) which gives which means that for large values of t and for positive v 1 , the solution behaves like that of the de Sitter universe. For the scale factor (70) we present in Fig. 2 the qualitative behavior of the Hubble function H (z), as well as the behavior of the equation of state parameter w (z) of the effective fluid in terms of the redshift 1 + z = 1 a .

Potential VD (φ)
For potential V D (φ), in the canonical coordinates {x, y} the constraint equation, i.e. the Hamiltonian of the dynamical system is written as where {p x , p y } = {ẏ,ẋ} and V 0 ,V 1 ,ρ m0 = 1 8B0 − 2 (V 0 , V 1 , ρ m0 ). The dynamical system admits the additional conservation law The Action which follows as a solution of the Hamilton Jacobi equation is calculated as such that the analytic solution of the field equations is given by the following system of two first-order ordinary differential equationsẋ where in the special case whereρ m0 = 0,V 0 = 0 it becomeṡ or equivalently, or A special solution of the latter equation is the power-law expression y t p , p = 2 , which leads to a power law scale factor. We remark that when B 0 = 1 16 for the potential The procedure that we follow to write the analytic solution for potential V E (φ) is based on the derivation of the Action by solving the Hamilton-Jacobi equation, as we did for potential V D (φ).
In the canonical coordinates {x, y} the constraint equation reads where V 0 ,V 1 ,ρ m0 = 1 8B0 − 2 (V 0 , V 1 , ρ m0 ). The quadratic conservation law admitted by the field equations is Consequently, the Action is calculated as where the reduced equations areẋ There are similarities of the latter solution with that of potential V D (φ), but for different value of the constant B 0 , specifically by replacing mathematically B 0 → i 4 B 0 . For B 0 = 1 16 in the canonical coordinates {u, v} the constraint equation becomes can be produced by a point-like Lagrangian. There are similarities with scalar-tensor theories although they are quite different theories. We have focused on the construction of scalar field potentials to see whether the field equations are Liouvilleintegrable, that means that the field equations can be solved in quadratures. Consequently, we investigated the functional forms of the scalar field potential where the field equations admit conservation laws quadratic in the momenta. By using the second conservation law we were able to write the analytic solution of the field equations for that specific scalar field potentials and whenever it was feasible, we have expressed the scale factor and the scalar field in closed-form functions.
We have found five families of scalar field potentials which are Liouville-integrable, and that admits conservation laws quadratic in the momenta, and they are in the form V A (φ) = V 0 φ p + V 1 φ r , where p, r are constants. For each dominant term of the potential the analytic solution for the scale factor behaves like a power-law function or like an exponential function which describes the de Sitter universe when the dominant power has the value two. We remark that we have selected a specific interaction function between the aether and the scalar fields. The interaction form that we selected has also geometric origins since for that function in the minisuperspace description, the dynamical variables of the field equations evolve in a two-dimensional space of maximally symmetry; in particular in two-dimensional flat space of Lorentzian signature. That is also a condition that we have assumed, in order the field equations to admit conservation laws quadratic in the momentum.
For some of the close-form solutions that we have found, we have studied the qualitative behavior of the Hubble factor, and we have presented the evolution of the equation of state parameter of the effective fluid in terms of redshift. From which we found that the effective fluid can cross the phantom divide line and behaves like a quintom field [3][4][5][6][7][8][9][10][11][12][13][14][15][16] or as a phantom field [17][18][19][20][21]. However the final attractor for that solutions is that of the de Sitter universe. More analysis should be done in that models in order to specify their physical viability, specifically if they can be contrasted against of cosmological observations. However such an analysis extends the scopes of this work and will be published elsewhere.
The transformation of the canonical variables in this model now is defined as As this point it is important to mention that the above results include the case where the Einstein-aether scalar field theory it is on a FLRW spacetime with nonzero spatially curvature K, that it is true when γ = 2 3 and ρ m0 = K.