Spatially Homogeneous Einstein-Aether Cosmological Models: Scalar Fields with a Generalized Harmonic Potential
B. Alhulaimi R. J. van den Hoogen A. A. Coley
gr-qc

[
    V.P. Gachok
    July 14, 2017
=================











§ INTRODUCTION
 



 §.§ Lorentz-violating Inflationary Cosmological Models


Lorentz invariance is a fundamental symmetry of General Relativity (GR) and standard particle physics and has been tested to a very high degree of accuracy.  However, there is growing evidence that some issues in quantum gravity (certain divergences, micro-causality) can be resolved by the removal of Lorentz invariance <cit.>. Furthermore, some approaches to quantum gravity may even desire a preferred rest frame in vacuum <cit.>.  On larger scales, explanations of Dark Energy and Dark Matter in the current cosmological paradigms based on GR might also be explained using an alternative theory of gravity <cit.>, such as those in which the Lorentz invariance requirement is relaxed. Indeed, in cosmology there is a natural frame associated with the cosmic microwave background, and therefore it is possible the Lorentz invariance assumption may be relaxed at late times.

The Einstein-Aether (AE) theory <cit.> is a proposed model of gravity which preserves general covariance and incorporates a violation of Lorentz invariance. The local space-time structure is determined by a dynamical time-like vector field, u^a (the aether), together with a metric tensor, g_ab. The field equations for this Einstein-Aether theory essentially consist of GR with a modified source due to the aether together with additional field equations describing the evolution of the aether vector field.  In standard Einstein-Aether theory, it is commonly assumed that the violation of Lorentz invariance is only within the gravity sector of the theory, while the matter sector continues to be coupled only to the metric, and hence remains Lorentz invariant.  However, it is natural to expect that Lorentz violations in the matter sector could also be permitted, albeit with the understanding that there are quite stringent constraints on such Lorentz violations <cit.>.

Assuming matter is determined by a scalar field, some researchers have investigated the potential changes that arise as a result of a violation of Lorentz invariance in the matter sector of the Einstein-Aether theory <cit.>. Barrow <cit.> investigated the effect of scalar field/aether field coupling in which the dependence of the scalar field in the potential was exponential in nature.  He found there are solutions with the possibility that the coupling parameter enables inflation in situations in which it would not otherwise occur. Sandin et al. used a similar ansatz for the scalar field potential in <cit.> and determined that there is a fundamental change in the future asymptotic behaviour when the coupling parameter becomes sufficiently large.  Donnelly and Jacobson <cit.> considered a chaotic inflationary scenario and determined that the coupling of the scalar field to the aether field can either slow down or speed up the evolution of inflation. Solomon and Barrow <cit.> completed a detailed analysis with no prescribed coupling between the scalar field and the aether field and found conditions on the potential that must be satisfied if one is to have stable slow-roll inflationary solutions.  Alhulaimi <cit.> generalized some of these results to include not only a coupling of the scalar field to the aether vector expansion, but also a coupling to the aether vector shear.  Where most others have coupled the scalar field to the aether field through the scalar field potential, Kanno and Soda <cit.> took a very different approach.   In their analysis, they assumed that the aether parameters c_i in equation (<ref>) are functions of the scalar field.  They found that it is possible to have inflation without a scalar field potential,  i.e., with a massless scalar field.

Inflation has become a well accepted, but not yet proven, mechanism which attempts to explain many cosmological issues <cit.>. A finite period of accelerated expansion (inflation) in the early universe is desirable to help address the isotropy, spatial homogeneity, horizon, and flatness problems <cit.>.   The standard inflationary model consists of a single massive scalar field that causes the universe to experience a period of exponential expansion early in its evolution.  A common inflationary scenario assumes a convex potential, such as a harmonic scalar field potential, in which inflation takes place during a period of “slow roll”, when the scalar field is decreasing very slowly in comparison to the expansion of the universe.  Not only is a period of accelerated expansion desirable at early times, but due to the Dark Energy phenomenon<cit.>, a period of accelerated expansion is also an attractive property to have at late times.

In this paper we shall investigate the dynamical evolution of a class of isotropic and anisotropic spatially homogeneous Einstein-Aether cosmological models containing a scalar field that is coupled to the aether field through the scalar field potential.  In particular, we explore the potential impact of Lorentz violation in the matter sector on the standard inflationary scenario <cit.>. More precisely, we study the inflationary scenario and investigate whether the inflationary solutions proposed <cit.> are stable when spatial curvature and anisotropy perturbations are considered. Further, we are also interested in the possibility of late time accelerated expansion in these models.




 §.§ Einstein-Aether Gravity


The action under consideration contains a Lagrangian describing Einstein-Aether gravity together with a Lagrangian for a matter field or fields (M)

    S=∫ d^4x√(-g)[1/8π Gℒ^AE+ℒ^M] .

The lagrangian ℒ^AE depends on the spacetime metric, g_ab, and the normalized aether vector field, u^a, and has the form  <cit.>:

    ℒ^AE=1/2R - K^ab_ab cd∇_au^c∇_bu^d  + λ(u^au_a+1)

where

    K^ab_abcd≡  c_1 g^ab g_cd + c_2δ_c^aδ_d^b+c_3δ_d^aδ_c^b+ c_4 u^a u^b g_cd.

We note that the parameters c_i defined here are the same as those used in <cit.> which are equal to half of the values of the c_i employed in <cit.> with an opposite sign for the c_4. In comparison to <cit.>, the c_i used here are 8π G times the values of c_i. The metric signature is assumed to be +2.

Let

    ℒ^U= - K^ab_ab cd∇_au^c∇_bu^d

then the variation of the action (<ref>) with respect to the aether vector field u^a yields

    -2λ u_a =  δℒ^U/δ u^a+8π G δℒ^M/δ u^a,

which when contracted with u^a, provides an explicit expression for the Lagrange multiplier

    2λ  =  δℒ^U/δ u^au^a+8π G δℒ^M/δ u^au^a.

Equation (<ref>) can then be used to eliminate the contribution of the Lagrange multiplier λ when calculating the effective energy momentum tensors.

Variation of the action (<ref>) with respect to g^ab, λ and a generalized matter field or fields Ψ, yields

    G_ab   =    T_ab^U + 8π G  T_ab^M,
    
    u^a u_a    =    -1, 
    δℒ^M/δΨ   =    0.

When the contributions from the Lagrange multiplier in equation (<ref>) are taken into account, expressions for the effective energy momentum tensors due to the aether vector field and the matter field become

    T_ab^U   =    -2δℒ^U/δ g^ab + g_abℒ^U+δℒ^U/δ u^cu^cu_au_b, 
    
    T_ab^M    =    -2δℒ^M/δ g^ab + g_abℒ^M+δℒ^M/δ u^cu^cu_au_b.


Given the form of the Lagrangian in equation (<ref>), the effective energy momentum tensor due to the Aether field <cit.> is

    T_ab^U   =    2∇_c(J_(a^acu^c_b) - J^c_c (au^c_b) -J_(ab)u^c ) 
            2c_1((∇_a u^c)(∇_b u_c) - (∇^c u_a)(∇_c u_b)) -2c_4u̇_au̇_b 
           -2( u^d∇_c J^c_cd+c_4u̇_cu̇^c   )u_au_b - g_ab(K^cd_cdef∇_c u^e ∇_d u^f),

where

    J^a_ab   =    -K^ac_acbd∇_c u^d,
    u̇^a    =    u^b∇_b u^a.





 §.§ Matter as a Scalar Field


Assuming that the matter component of the universe is a single scalar field having a potential that is assumed to be a function of the scalar field together with the expansion and shear scalars of the aether vector field, the matter Lagrangian becomes:

    L^M= -1/2 g^a b ∇_aϕ∇_bϕ - V(ϕ,θ,σ^2),

where θ=∇_a u^a is the expansion scalar and σ^2=1/2σ_abσ^ab is the shear scalar.  Again taking into account contributions from the Lagrange multiplier, equation (<ref>) yields the effective energy momentum tensor due to the scalar field

    T_ab^M   =   ∇_a ϕ∇_b ϕ - (1/2∇_aϕ∇^a ϕ + V )g_ab
    + θ V_θ g_ab+V̇_θ h_ab
    
           +(θ V_σ^2+V̇_σ^2)σ_ab+V_σ^2σ̇_ab-2σ^2V_σ^2u_au_b

where the terms V_θ and V_σ^2 are the partial derivatives of the scalar field potential with respect to θ and σ^2, respectively, and h_ab≡ g_ab+u_au_b.  If there is no coupling between the aether field and the scalar field via the potential, then V_θ=V_σ^2=0 and the energy momentum tensor reduces to the standard form for a minimally coupled scalar field.
In addition, the field equation (<ref>) yields the Klein-Gordon equation for the scalar field:

    ∇^a∇_aϕ-V_ϕ=0.










§ ISOTROPIC EINSTEIN-AETHER MODELS COUPLED TO A SCALAR FIELD
 



 §.§ The Spatially Homogeneous and Isotropic Model


We shall assume that the spacetime is spatially homogeneous and isotropic having spacetime coordinates [t,r,θ,ψ] and a metric of the form:

    ds^2=-dt^2+a(t)^2(1/1-kr^2dr^2+r^2dθ^2+r^2sin^2(θ)dψ^2),

where k takes on values {-1,0,1} for negative, zero and positive spatial curvature, respectively.  In a spatially homogeneous and isotropic cosmological model with comoving time, the aether vector field necessarily coincides with the rest frame defined by the Hubble expansion.  Specifically, this implies that in spatially homogeneous and isotropic models that the aether vector must be orthogonal to the three-dimensional spatial hyper-surfaces and takes the form u^a=(1,0,0,0).

With the above assumptions on the metric and the aether vector, the shear, the vorticity and the acceleration of the aether vector are zero and the covariant derivative

    ∇_b u_a=1/3θ (g_ab+u_au_b)

is simply determined by the expansion scalar

    θ = ∇_au^a = 3ȧ/a.


With the definition of T_ab^U in equation (<ref>), the effective energy density ρ^U and isotropic pressure, p^U,   due to the aether field are

    ρ^U    =     -1/3c_θθ^2 ,
    
    p^U    =    1/3c_θθ^2 +2/3c_θθ̇.

Where a new parameter c_θ= (c_1 +3c_2 +c_3), defined before in <cit.>, allows for some efficiencies in notation since the field equations are independent of any other linear combinations of the c_i.
The Einstein-aether field equations reduce to the following:

    0   =    -1/3(1+c_θ)θ^2+8π Gρ^M-3k/a^2, 
    
    0    =    -(1+c_θ)θ̇-1/3(1+c_θ) θ^2  -8π G/2(ρ^M+3p^M ),

where there still exists the freedom to choose some appropriate units.   Without loss of generality, new units can be chosen so that 8π G/1+c_θ=1 in which case the explicit dependence of the field equations on the aether parameter c_θ has been eliminated.



 §.§ The Scalar Field Potential

We shall consider a class of quadratic scalar field  potentials  of the form

    V(ϕ, θ)=1/2 m^2 ϕ^2 + μθϕ,

where the scalar field/Aether field coupling term μθϕ term can be interpreted as an incorporation of an external source, in this case the Aether, acting on the scalar field.  The effective energy density ρ^M and isotropic pressure p^M from equation (<ref>), are

    ρ^M    =    1/2ϕ̇^2+ 1/2 m^2 ϕ^2,
    
    p^M    =    1/2ϕ̇^2- 1/2 m^2 ϕ^2+μϕ̇.

The Klein-Gordon equation becomes

    0=ϕ̈+θϕ̇+m^2ϕ+μθ,

where we can more clearly see how μθ acts like an external source in the Klein-Gordon equation when compared to the usual non-coupled μ=0 version of the equation.





 §.§ The Dynamical System


The Einstein-aether field equations and the Klein-Gordon equation can be expressed as the following system of ordinary differential equations

    θ̇   =    -1/3θ^2 +m^2/2ϕ^2-ψ^2-3/2μψ,
    
    ϕ̇   =   ψ,
    
    ψ̇   =    -θψ-m^2ϕ-μθ.

with first integral

    θ^2/3 = m^2/2ϕ^2+1/2ψ^2-3k/a^2.


Equations (<ref>)-(<ref>), therefore, yield a three dimensional dynamical system for the variables (θ, ϕ,ψ)  depending on three parameters (k, m, μ) having a first integral given by equation (<ref>). Since the system of equations is invariant under the transformation (μ,ϕ,ψ) ↦ -(μ,ϕ,ψ), we can without loss of generality, assume that μ≥ 0.  Given that the phase space for the dynamical system defined in equations (<ref>)-(<ref>) with first integral (<ref>) is not  bounded, we employ dimensionless variables <cit.> which will transform the system into an autonomous system of differential equations on a bounded phase space.




 §.§ Qualitative Analysis




  §.§.§ Introducing Normalized Variables


Introducing a time variable  τ

    dτ/dt=√(1+θ^2)

and normalized variables

    D    ≡   θ/√(1+θ^2), 
    Φ   ≡   √(3/2)(m ϕ/√(1+θ^2)),
    Ψ   ≡   √(3/2)(ϕ̇/√(1+θ^2)),

the evolution equations (<ref>)-(<ref>) become

    D^'   =    (1-D^2)𝒳,  
    Φ^'   =    m Ψ√(1-D^2)-D Φ𝒳,  
    Ψ^'   =    -D Ψ- √(1-D^2)( m Φ+√(3/2)μ D ) -Ψ D 𝒳,

where the prime here indicates the differentiation with respect to τ and 𝒳 is given by the expression

    𝒳   =   θ̇/θ^2+1
       =    -1/3D^2-2/3Ψ^2+1/3Φ^2 -√(3/2)μΨ√(1-D^2)


The Friedmann equation (<ref>) becomes

    D^2-Φ^2-Ψ^2= -9 k/ a^2(1+θ^2).

Further, if k = -1,0, then it follows that

    0 ≤Φ^2+Ψ^2  ≤ D^2≤ 1.

That is, D, Φ, Ψ are bounded in the flat and negatively curved scenarios and the phase space is a compact set.  Hence forward, we shall restrict our analysis to k=0, -1 cases only.



  §.§.§ Invariant Sets and Monotonic Functions


The phase space can be subdivided into four disjoint invariant sets according to the curvature of the model and whether D=1 (θ→∞) or not.  A superscript “-” indicates that points in this set represent negatively curved models, while a superscript “0” indicates a flat model. The invariant sets are

    ^-    =   {(D,Φ,Ψ)|D<1,Φ^2+Ψ^2  <   D^2},
    ^0    =   {(D,Φ,Ψ)|D<1,Φ^2+Ψ^2  =   D^2},
    ^-    =   {(D,Φ,Ψ)|D=1,Φ^2+Ψ^2  < 1},
    ^0    =   {(D,Φ,Ψ)|D=1,Φ^2+Ψ^2  = 1},

the dimensions of which are 3, 2, 2, and 1, respectively.  We note the following closure properties of the sets

    ^-   =   ^-∪^0∪^-∪^0,
    ^-   =   ^-∪^0,
    ^0   =   ^0∪^0.

Further the invariant set ^- can be divided into three distinct pieces depending on whether Φ<0, Φ=0 or Φ>0.

If we define Λ_1 = D^2- Φ^2-Ψ^2 and Λ_2=D^2-1 then

    Λ_1^'/Λ_1   =    -2/3  D  ( 3𝒳+1 ),
    Λ_2^'/Λ_2   =    2D𝒳.

in which case the non-negative function W=(Λ_1)^2(Λ_2)^2 has the derivative

    W^' = -4/3WD.

Since W> 0 and W^'<0  in the set ^- we can conclude that there are no periodic orbits in this 3-dimensional invariant set. This implies that there are no equilibrium points in the set ^-, and any equilibrium points of the autonomous system of differential equations (<ref>)-(<ref>) will lie in the lower dimensional invariant sets ^0, ^-  or ^0.  We also note that in the invariant set ^-∩{Φ<0}, one can show that Φ^' <0, and similarly in the set ^-∩{Φ>0}, one can show that Φ^' >0.  This shows that there are no closed or periodic orbits in these sets.  The remaining portion ^-∩{Φ=0}, is 1 dimensional. No monotonic function has been found in the set ^0 and consequently the most interesting qualitative behaviour for this autonomous system of differential equations occurs in ^0.



  §.§.§ Equilibrium Points

The equilibrium points and a non-isolated line of equilibria for the system (<ref>)-(<ref>) with the value of 𝒳  and their stability are summarized in Table (<ref>).




  §.§.§ Stability of Equilibrium point P_0


Evaluating the linearization matrix of the system (<ref>)-(<ref>) at P_0 gives us the following eigenvalues

    λ_1   =0,
       λ_2,3   = ±√(6)/2√(μ^2-μ_c^2).

Note that, if μ > μ_c then P_0 is a saddle. But, if  μ < μ_c then  all the eigenvalues  have  zero real part which implies that the local qualitative behaviour at P_0 is not determined by its linearization.  However a perturbative solution near P_0  can be found, and fortunately an analysis of the first order solution is sufficient to determine the local stability of P_0 when μ < μ_c.


We first introduce new scaled variables (d,ϕ,ψ) such that

    D=ϵ(d-μ/μ_cϕ),  Φ=ϵϕ,  Ψ=ϵψ,

where ϵ is assumed to be small, to determine a leading order approximation to the solution of the equations near P_0.  We note that the ϕ and ψ variables that are employed in this subsection are not the original variables used to describe the scalar field and its derivative.  Using our new dependent variables (<ref>), and expanding (<ref>)-(<ref>) as a power series in ϵ we derive the following

    d^'   =   ϵ/3(-d^2-2ψ^2+ ϕ^2+2μ/μ_cdϕ-μ^2/μ_c^2ϕ^2)+O(ϵ^2), 
    ϕ^'   =   √(6)/2μ_cψ +O(ϵ^2), 
    ψ^'   =   √(6)/2(- μ d -μ_cϕ +μ^2/μ_cϕ)+ ϵ(-d ψ+μ/μ_cϕψ)+O(ϵ^2),

where we kept only terms up to linear order in ϵ.  To proceed with the construction of a perturbative solution, we employ the method of multiple scales <cit.>.


In the method of multiple scales with two time scales, the original fast time τ and a second slow time η=ϵτ, each dependent variable is expressed as

    x ≡    x(τ, η)=    x_0 (τ, η)+ϵ    x_1 (τ, η) + O(ϵ^2)

and using the chain rule, derivatives become expanded as

    x^'=    x_0τ+   ϵ (x_0η      + x_1τ )  + O(ϵ^2).

The (') indicates the ordinary derivative of the variable with respect to time τ while the subscripts τ and η denote partial derivatives. A valid perturbative solution is obtained by ensuring that the solution remains bounded at all orders of ϵ.

Using equation (<ref>) and (<ref>) for variables (d,ϕ,ψ) and substituting into  (<ref>) and matching powers of ϵ yields the following system of partial differential equations for the zeroth order [O(ϵ^0)] terms

    d_0τ   =     0, 
    
    ϕ_0τ   =   √(6)/2μ_cψ_0, 
    
    ψ_0τ   =   √(6)/2(-μ d_0 - μ_cϕ_0+μ^2/μ_cϕ_0),

and the following system of partial differential equations for the first order [O(ϵ^1)] terms

    d_1τ   =   1/3(-d_0^2-2ψ_0^2 + ϕ_0^2 +2μ/μ_cd_0ϕ_0-μ^2/μ_c^2ϕ_0^2)-d_0 η, 
    
    ϕ_1τ   =   √(6)/2μ_cψ_1-ϕ_0 η,
    
    ψ_1τ   =   √(6)/2(-μ d_1-μ_cϕ_1 +μ^2/μ_cϕ_1)
    +(-d_0ψ_0+μ/μ_cϕ_0ψ_0)-ψ_0 η.

Solving the partial differential equations for the Zeroth order terms yields

    d_0(τ,η)    =    B(η), 
    ϕ_0(τ,η)    =    A(η)cos(λτ-Λ(η))-μμ_c/μ_c^2-μ^2B(η),
    ψ_0(τ,η)    =    -√(6)λ/3μ_cA(η)sin(λτ-Λ(η)),

where λ=√(6)/2√(μ_c^2-μ^2) and A(η), B(η) and Λ(η) are as yet undetermined functions of the slow time η.  Solving the partial differential equations for the first order terms, and restricting ourselves to only bounded solutions, determines a set of ordinary differential equations for the unknown functions A(η), B(η) and Λ(η),

    A_η   =   -1/2μ_c^2/μ_c^2-μ^2AB, 
    
           B_η   =   -1/3(μ_c^2/μ_c^2-μ^2B^2 + μ_c^2-μ^2/2μ_c^2A^2),
    Λ_η   =    0.


Therefore, in terms of the original variables the first term of the perturbative solution is

    D(τ)    =   ϵ(- μ/μ_cA(η)cos(λτ-Λ(η))+μ_c^2/μ_c^2-μ^2B(η) ), 
    Φ(τ)    =   ϵ(  A(η)cos(λτ-Λ(η))-μμ_c/μ_c^2-μ^2B(η)),
    Ψ(τ)    =   ϵ(-√(μ_c^2-μ^2)/μ_cA(η)sin(λτ-Λ(η))),

where the functions A(η), B(η) and Λ(η) satisfy the differential equations (<ref>), and due to (<ref>) are bounded by

    B(η)^2≥(μ_c^2-μ^2)^2/μ_c^4A(η)^2,

where we note that if B(η)→ 0 then we also have A(η)→ 0.

We are interested in determining the asymptotic behaviour as τ→∞.  We observe that the phase shift Λ(η) is a constant and has no effect on the future dynamics.  The fast time τ essentially describes the oscillations of the scalar field, which to first order in ϵ has a period of T=2π/λ.  We note that the period of these oscillations T ∼ 1/√(μ_c^2-μ^2), gets longer as the strength of the coupling parameter μ is increased towards μ_c.

We also observe that the amplitude of the oscillations A(η), and the vertical shift B(η) are functions of the slow time η and consequently the amplitude and vertical shift drift slowly in comparison to the oscillatory changes. For initial values of B(η)>0 we see that both A(η),B(η)→ 0 as η→∞.  That is, the amplitude of the oscillations and the vertical shift both slowly decrease to zero, indicating that the point P_0 is stable when μ<μ_c.




  §.§.§ Stability of Equilibrium points in ^-∪^0


Unfortunately, while we have an autonomous system of differential equations defined on a compact set, the system is not differentiable at any points in the invariant set D^-∪ D^0. We note that any equilibrium points in the invariant set D^-∪ D^0  represent asymptotic states in which θ→∞.  In order to determine the local qualitative behaviour at these equilibrium points, we replace variable D with

    T=1/√(1+θ^2)=√(1-D^2).

The evolution equations (<ref>)-(<ref>) become

    T^'   =    -T√(1-T^2)𝒳, 
    Ψ^'   =    -√(1-T^2)Ψ- T ( m Φ+√(3/2)μ√(1-T^2)) -Ψ√(1-T^2)𝒳,
    Φ^'   =    m Ψ T-√(1-T^2)Φ𝒳,

with

    𝒳 = -1/3(1-T^2)-2/3Ψ^2+1/3Φ^2 -√(3/2)μΨ T.

The value D=1 for the equilibrium points P_1,2,3,4,5 is simply replaced with T=0.  With this transformation we are able to locally determine the qualitative behaviour of each the equilibrium points in the invariant set D^-∪ D^0.

The eigenvalues of the points P_1 and P_2 are 1,4/3, 1 which implies these points are generally sources.  The eigenvalues of the point P_5 is 1/3,1/3, -2/3 which implies this point is generally a saddle.  Further, the eigen-directions that span the T=0 invariant set, are associated with one positive and one negative eigenvalue.  Therefore this equilibrium point is a saddle in the T=0 set.

The eigenvalues of the point P_3 and P_4 are 0,-1, -2/3 which implies that we cannot determine the general behaviour of this point without resorting to additional analysis.  However, the eigen-directions that span the T=0 invariant set, are associated with the two negative eigenvalues.  Therefore these equilibrium points are sinks in the T=0 set.

To complete the analysis of the qualitative behaviour near P_3 and P_4 we calculate the center manifold <cit.>.  In this case the center manifold is a one dimensional curve that must lie in the ^0 invariant set.  The center manifold for P_3 can be parameterized as

    T    =    T 
    Φ   =    1 - (1/2 + 3/4(μ+μ_c)^2)T^2 + O(T^4)
    Ψ   =    -√(6)/2(μ+μ_c)T+3/8√(6)μ_c(μ+μ_c)^2T^3+O(T^4)

The leading order term of the dynamical system restricted to the center manifold reduces to

    T'= 3/2μ_c(μ+μ_c)T^3.

Since T'>0 for T>0, P_3 is unstable along its center manifold. It is therefore a saddle in the full three dimensional phase space.


The center manifold for P_4 can be parameterized as

    T    =    T 
    Φ   =    -1 + (1/2 + 3/4(μ-μ_c)^2)T^2 + O(T^4)
    Ψ   =    -√(6)/2(μ-μ_c)T-3/8√(6)μ_c(μ-μ_c)^2T^3+O(T^4)

The leading order term of the dynamical system restricted to the center manifold reduces to

    T'= -3/2μ_c(μ-μ_c)T^3.

If μ<μ_c then T'>0 for T>0 and P_4 is unstable along its center manifold.  However, if μ>μ_c then T'<0 for T>0 and P_4 is stable along its center manifold. Therefore P_4 is a saddle when μ<μ_c and a sink when μ>μ_c in the full three dimensional phase space.



  §.§.§ The Bifurcation Value


When μ=μ_c, there is a line of equilibria given by (D,Φ,Ψ)=(s, -s,0) where 0<s<1 having endpoints P_0 and P_4.  The eigenvalues of points on this line of equilibria are 0, -s, -2/3s which implies that points on this line are attractors when it exists. We observe that as the parameter μ increases towards its bifurcation value, μ_c, the stability of the point P_0 is transferred to the point P_4 via this line of equilibria.



  §.§.§ Heteroclinic Sequences


Very often one is not only interested in the past and future behaviour of a system of differential equations, but one is also interested in the intermediate behaviour of the system.  One technique to analyze the intermediate behaviour is to describe the heteroclinic sequences that are possible <cit.>.  We note that for every heteroclinic sequence there exists a set of orbits that are arbitrarily close to that sequence.  Figures (<ref>) describes the possible heteroclinic sequences that are possible. Again we see how the intermediate and future behaviour changes as the parameter μ changes.






 §.§ Inflation and Accelerated Expansion

As an indicator of the accelerated expansion we introduce the deceleration parameter

    q ≡ - a ä/ȧ^2=- (3θ̇/θ^2+1).

It follows that the deceleration parameter can also be expressed in terms of the normalized bounded variables as follows;

    q= -1/D^2(-2 Ψ^2 +Φ^2-3√(3/2)μ√(1-D^2)Ψ).

The sign of the deceleration parameter indicates the nature of the expansionary evolution. If q > 0, then the cosmological expansion is decelerating, while negative values of q indicate an accelerating or inflationary dynamics. See Table (<ref>) for a summary of the values of q for each equilibrium point.




 §.§ Numerical Analysis


It is constructive to illustrate a few numerical solutions for the three different regimes of future asymptotic behaviour, μ<μ_c [see Figure (<ref>)], μ=μ_c [see Figure (<ref>)], and μ>μ_c [see Figure (<ref>)]. In each case the integrations are done in the full 3-dimensional phase space. A few initial conditions are selected to show different past and future asymptotic behaviours.











 §.§ Observations


In the spatially homogeneous and isotropic case analyzed here we observe that the past dynamics are independent of the strength of the coupling parameter μ. We find in the zero curvature and the negative curvature models that there are two possible asymptotic behaviours to the past, one in which (D,Φ,Ψ)→ (1,0,1), P_1, and one in which (D,Φ,Ψ)→ (1,0,-1), P_2.  These past attractors represent a massless scalar field FRW model <cit.>.


We also observe that the future asymptotic state depends on the strength of the coupling parameter μ.  For weak coupling of the Aether field to the scalar field, i.e, , μ<μ_c=√(2/3)m, the dynamics are similar to that found when there is no coupling of the Aether field to the scalar field, i.e., when μ=0.  If μ<μ_c then P_0 is the stable attractor: orbits oscillate and slowly decay in amplitude towards this final non-inflationary asymptotic state.

For strong coupling of the Aether field to the scalar field, μ>μ_c, the dynamics are very different.  When μ>μ_c the scalar field does not come to rest at the minimum of the potential: the strength of the Aether interaction forces a different final equilibrium state.  If μ>μ_c then in both the zero curvature and the negative curvature models we find that the stable attractor in these models changes from P_0 to P_4.  In this scenario, we find that the square of the scalar field and the square of the expansion scalar scale together as

    m^2/2ϕ^2 ∼1/3(1+θ^2)

and consequently grow without bound.  We also observe that all orbits (excepting for the exceptional orbits) experience some period of accelerating expansion as they evolve to their final asymptotic state, P_4.  In the next section we will add anisotropy to these models to determine if the addition of anisotropy changes these observations.









§ ANISOTROPIC EINSTEIN-AETHER MODELS COUPLED TO A SCALAR FIELD
 



 §.§ A Class of Diagonal Bianchi Models


In order to investigate Einstein-Aether cosmological models with a scalar field having a potential with interaction terms that contain both the expansion and shear of the aether vector requires a broader class of space-time geometries which include anisotropy.  For our purposes, the one parameter family of spatially homogeneous and anisotropic diagonal Bianchi type VI_h space-times provides an arena to determine the effect of adding a shear interaction term into the scalar field potential in the Einstein-Aether models studied in the previous section.   The metric is assumed to have the form

    ds^2=-dt^2+a(t)^2dx^2+b(t)^2e^2(h-1)xdy^2+c(t)^2e^2xdz^2.

There are three special classes are worth mentioning: if the parameter h=2 then the metric is Bianchi type V which has the negatively curved isotropic models (see previous section) as a subcase, if h=1 then it is Bianchi type III, and if h=0 then it is Bianchi type VI_0.

In a spatially homogeneous and isotropic cosmological model with comoving time, the aether vector field necessarily coincides with the rest frame defined by the Hubble expansion. Therefore in deviations from spatially homogeneity and isotropy, one could assume that the preferred frame for the aether approximately coincides with the cosmological rest frame defined by the Hubble expansion. However, if one relaxes only the isotropy assumption, then the aether vector would be tilted away from the hyper-surface normal of the spatially homogeneous hyper-surfaces which could allow for a richer set of dynamical behaviours. However, it has been argued in <cit.> that during cosmological evolution, the two frames will come into alignment. Given these arguments, we will assume that the aether vector is aligned with the hyper-surface unit normal to the surfaces of homogeneity and is of the form u^a=(1,0,0,0).

With the above assumptions on the metric and the aether vector, the vorticity and the acceleration of the aether vector are zero, and the covariant derivative

    ∇_b u_a=σ_ab+1/3θ(g_ab+u_a u_b),

is simply determined by the expansion scalar

    θ = ∇_au^a = ȧ/a+ḃ/b+ċ/c,

and the shear tensor

    σ_ab=u_(a;b)-1/3θ (g_ab+u_a u_b).

The shear tensor has the form σ^a_ab=[0,σ_1,σ_2,-(σ_1+σ_2)] where

    σ_1   =   1/3(2ȧ/a-ḃ/b-ċ/c),
    σ_2   =   1/3(-ȧ/a+2ḃ/b-ċ/c).

We note that the shear scalar is

    σ^2=1/2σ_abσ^ab=σ_1^2+σ_2^2+σ_1σ_2.


With the definition of T_ab^U in equation (<ref>), the effective energy density ρ^U, isotropic pressure, p^U, energy flux q_a^aU, and anisotropic stress π_ab^abU   due to the aether field are

    ρ^U    =     -1/3c_θθ^2 -2c_σσ^2,
    
    p^U    =    1/3c_θθ^2 +2/3c_θθ̇-2c_σσ^2,
    
    q_a^aU   =   0,
    π_ab^abU    =    2c_σ (σ̇^a_ab+θσ^a_ab),

Where new parameters c_θ= (c_1 +3c_2 +c_3) and c_σ=c_1+c_3, defined before in <cit.>, allow for some efficiencies in notation since the field equations are independent of any other linear combinations of the c_i.

The Einstein-aether field equations reduce to the following set of linearly independent equations

    0   =    -1/3(1+c_θ)θ^2+(1-2c_σ)σ^2+h^2-h+1/a^2+8π Gρ^M, 
    
    0    =    -(1+c_θ)θ̇-1/3(1+c_θ) θ^2 -2(1-2c_σ)σ^2 -8π G/2(ρ^M+3p^M ),
    
    0   =    -(h+1)σ_1+(h-2)σ_2-8π G/a^2 q_2^aM
    
    0    =    (1-2c_σ)σ̇_1+(1-2c_σ)θσ_1 +(h-2)^2/3a^2 -8π Gπ^M_1, 
    
    0    =    (1-2c_σ)σ̇_2+(1-2c_σ)θσ_2 +(h-2)(h+1)/3a^2 -8π Gπ^M_2,

where there still exists the freedom to choose some appropriate units.  Without loss of generality, one can choose new units so that 8π G/1+c_θ=1, and a new parameter C=1-2c_σ/1+c_θ, in which case the explicit dependence of the field equations on the aether parameter c_θ has again been eliminated.  Assuming that GR and Einstein Aether theory have equivalent PPN parameters and that there is stable positive energy modes but no vacuum Čerenkov radiation imposes some constraints on the values of c_θ and C.  See Appendix <ref> for details.



 §.§ The Scalar Field Potential

We shall consider a class of quadratic scalar field potentials  of the form

    V(ϕ, θ)=1/2 m^2 ϕ^2 + μθϕ+νσϕ.

In this case, from equation (<ref>), the effective energy density ρ^M, isotropic pressure p^M, energy flux q_a^aM, and anisotropic stress π_ab^abM due to the scalar field are

    ρ^M    =    1/2ϕ̇^2+ 1/2 m^2 ϕ^2,
    
    p^M    =    1/2ϕ̇^2- 1/2 m^2 ϕ^2+μϕ̇-νϕσ,
    
    q_a^aM    =     0 
    π_ab^abM    =    (νϕ/2σ) (σ̇^a_ab+θσ^a_ab)
     +d/dt(νϕ/2σ)σ^a_ab.

where we define π^M_i such that π_ai^ibM=[0,π^M_1,π^M_2,-(π^M_1+π^M_2)].
The final field equation comes from the Klein-Gordon equation which becomes

    0=ϕ̈+θϕ̇+m^2ϕ+μθ+νσ.





 §.§ The Dynamical System


Given that there is no energy flux either from the existence of the aether field or from the non-minimal coupling of the aether to the scalar field potential, we are able to use equation (<ref>) to rewrite all equations in terms of the anisotropy scalar σ:

    σ_1 = h-2/√(3)√(h^2-h+1)σ,     σ_2 = h+1/√(3)√(h^2-h+1)σ.

Hence, the final form of the Einstein-aether field equations and the Klein-Gordon equation yield the following dynamical system

    θ̇   =    -1/3θ^2 -2Cσ^2-ψ^2+m^2/2ϕ^2-3/2μψ+3/2νϕσ,
    
    σ̇   =   -σθ+ν/2C(ψ+θϕ)
              +(h-2)/√(3)√(h^2-h+1)1/C(1/3θ^2-Cσ^2-1/2ψ^2-m^2/2ϕ^2),
    
    ϕ̇   =   ψ
    
    ψ̇   =    -θψ-m^2ϕ-μθ-νσ.

with first integral

    θ^2/3 = m^2/2ϕ^2+1/2ψ^2+Cσ^2+h^2-h+1/a^2,


Equations (<ref>)-(<ref>), therefore, yield a four dimensional dynamical system for the variables (θ,σ, ϕ,ψ)  depending on five parameters (h, C, m, μ, ν) having a first integral given by equation (<ref>). We restrict our analysis to the diagonal Bianchi type V models (h=2).  We note that even in this very simple anisotropic case when we compare the evolution equation for the shear with would would happen if the scalar field potential did not have an interaction term, we notice that the negatively curved FRW models are no longer an invariant subset of the Bianchi type V system (<ref>)-(<ref>).  Since the system of equations when h=2 is invariant under the transformations (μ,σ,ϕ,ψ) ↦ -(μ,σ,ϕ,ψ), and (ν,σ) ↦ -(ν,σ), we can without loss of generality, assume that both μ≥ 0 and ν≥0.  Given that the phase space for the dynamical system defined in equations (<ref>)-(<ref>) with first integral (<ref>) is not  bounded, we employ dimensionless variables <cit.> which will transform the system into an autonomous system of differential equations on a bounded phase space.



 §.§ Qualitative Analysis




  §.§.§ Introducing Normalized Variables


Introducing a time variable τ

    d τ/dt≡√(1+θ^2),

and  normalized variables

    D          ≡   θ/√(1+θ^2)
    Σ   ≡   √(3)σ/√(1+θ^2),
    Φ   ≡   √(3/2)(m ϕ/√(1+θ^2)),
    Ψ   ≡   √(3/2)(ϕ̇/√(1+θ^2)).

then, the Bianchi type V evolution equations in equations (<ref>)-(<ref>) when h=2  become

    D^'   =    (1-D^2)𝒳,
    Φ^'   =    m √(1-D^2)Ψ-D Φ𝒳
    Ψ^'   =    -D Ψ(1+𝒳)- √(1-D^2)[  m Φ+√(3/2)μ D+ ν/√(2)Σ],
    Σ^'   =    - Σ D(1+  𝒳) + ν/√(2)  C ( Ψ√(1-D^2) +D Φ/m),

where  the prime here indicates the differentiation with to respect to τ.   𝒳  is given by the  expression

    𝒳   =   θ̇/1+θ^2
       =   -1/3D^2-2/3Ψ^2+1/3Φ^2 -2/3 C Σ^2 -√(3/2)μ√(1-D^2)Ψ + ν/√(2) mΣΦ.

The Friedmann equation (<ref>) when h=2 becomes

    D^2-Φ^2-    C Σ^2 -Ψ^2=9/a^2 ( 1+ θ^2).

where if C≥0, then it follows that

    0 ≤Φ^2+Ψ^2 + C Σ^2  ≤ D^2≤ 1.

That is, D, Φ, Ψ, Σ are bounded if C>0 and the phase space is a bounded set.  Hence forward, we shall restrict our analysis to C≥0 case only.



  §.§.§ Invariant Sets and Monotonic Functions

The phase space can be subdivided into four disjoint invariant sets according to the curvature of the model and whether D=1 (θ→∞) or not.  A superscript “-” indicates that points in this set represent negatively curved models, while a superscript “0” indicates a flat model. The invariant sets are

    ^-    =   {(D,Φ,Ψ,Σ)|D<1,   Φ^2+Ψ^2+CΣ^2  < D^2},
    ^0    =   {(D,Φ,Ψ,Σ)|D<1,   Φ^2+Ψ^2+CΣ^2  = D^2},
    ^-       =   {(D,Φ,Ψ,Σ)|D=1,   Φ^2+Ψ^2+CΣ^2  < 1},
    ^0       =   {(D,Φ,Ψ,Σ)|D=1,   Φ^2+Ψ^2+CΣ^2  = 1}.

The dimensions of which are 4, 3, 3, and 2, respectively.  The invariant set ^- represents the Bianchi type V models while the invariant set ^0 represents the Bianchi type I models.  Again, the shear interaction term in the scalar field potential plays a significant role, since Σ=0 will be an invariant set only if ν=0.

Similar to the isotropic case, if we define Λ_3 = D^2- Φ^2-Ψ^2-CΣ^2 and Λ_4=D^2-1 then

    Λ_3^'/Λ_3   =    -2/3  D  ( 3𝒳+1 ),
    Λ_4^'/Λ_4   =    2D𝒳.

in which case the non-negative function W=(Λ_3)^2(Λ_4)^2 has the derivative

    W^' = -4/3WD.

Since W> 0 and W^'<0  in the set ^- we can conclude that there are no periodic orbits in this 4-dimensional invariant set. This also implies that there are no equilibrium points in the set ^-, and any equilibrium points of the autonomous system of differential equations (<ref>)-(<ref>) will lie in the lower dimensional invariant sets ^0, ^-  or ^0.





  §.§.§ Equilibrium Points

The equilibrium points for  the system (<ref>)-(<ref>) with the value of 𝒳  and their stability are summarized in Table (<ref>). Comparing the equilibrium points in Table (<ref>) with those found in Table (<ref>) we see that the points P_0 and P_5 represent the same equilibrium states in both tables. The points P_1 and P_2 represents the same equilibrium state in both tables, however, in Table (<ref>) they are actually two special (isotropic) points on a non-isolated circle of generally non-isotropic equilibria given by C^*.  The points P_3 and P_4 in Table (<ref>) reduce to P_3 and P_4 in Table (<ref>) when ν=0.  Further, the dynamical behaviour and stability is analogous to the stability of P_3 and P_4 in Table (<ref>). Similar to the isotropic case, there is a line of equilibria L_04 when μ=μ_cν that connects P_0 and P_4.  We note that the local stability of the equilibrium points depends on the bifurcation value μ_cν^2=μ_c^2 +ν^2/3C. Recall μ_c=√(2/3)m is the bifurcation value found in the isotropic case studied earlier, and so since C and ν are both positive, the bifurcation value for the anisotropic case, is always a bit larger than in the isotropic case. Further, by choosing smaller and smaller values of C, one can increase the value of bifurcation value μ_cν.







  §.§.§ Stability of Equilibrium Point P_0


Evaluating the linearization matrix of the system (<ref>)-(<ref>) at P_0 gives us the following eigenvalues

    λ_1,2   =0,
       λ_3,4   = ±√(6)/2√(μ^2-μ_c ν ^2).

Note that, if μ >μ_c ν then P_0 is a saddle. But, if  μ < μ_c ν  then  all the eigenvalues  have  zero real part which implies that the local qualitative behaviour at P_0 is not determined by its linearization.  However a perturbative solution near P_0  can be found, and fortunately an analysis of the first order solution is sufficient to determine the local stability of P_0 when μ  < μ_c ν .


We first introduce new scaled variables (d,ϕ,ψ,σ) such that

    D=ϵ(d-μ/μ_cϕ),  Φ=ϵϕ,  Ψ=ϵψ,   Σ= ϵ(σ +ν/√(3) C μ_cϕ)

where ϵ is assumed to be small, to determine a leading order approximation to the solution of the equations near P_0.  We note that the ϕ,ψ and σ variables that are employed in this subsection are not the original variables used to describe the scalar field and its derivative.  Using our new dependent variables (<ref>), and expanding (<ref>)-(<ref>) as a power series in ϵ we derive the following

    d^'   =   ϵ/3(-d^2-2ψ^2+ ϕ^2+2μ/μ_cdϕ-μ^2/μ_c^2ϕ^2 +ν^2 / 3 C μ_c^2ϕ^2 -2 C σ^2 -√(3)ν/ 3 μ_cϕσ)+O(ϵ^2), 
    ϕ^'   =   √(6)/2μ_cψ +O(ϵ^2), 
    ψ^'   =   √(6)/2(- μ d -μ_cϕ +μ^2/μ_cϕ - ν/√(3)σ -ν^2/ 3 C μ_cϕ)+ ϵ(-d ψ+μ/μ_cϕψ)+O(ϵ^2),
    σ^'    =   ϵ( -σ d +μ/μ_cσϕ).

where we kept only terms up to linear order in ϵ.  To proceed with the construction of a perturbative solution, we employ the method of multiple scales <cit.>.

Using equation (<ref>) and (<ref>) for variables (d,ϕ,ψ,σ) and substituting into  (<ref>) and matching powers of ϵ yields the following system of partial differential equations for the zeroth order [O(ϵ^0)] terms

    d_0τ   =     0, 
    
    ϕ_0τ   =   √(6)/2μ_cψ_0, 
    
    ψ_0τ   =   √(6)/2(- μ d_0 -μ_cϕ_0 +μ^2/μ_cϕ_0 - ν/√(3)σ_0 -ν^2/ 3 C μ_cϕ_0),
    σ_0τ   =     0

and the following system of partial differential equations for the first order [O(ϵ^1)] terms

    d_1τ   =   1/3(-d_0^2-2ψ_0^2+ ϕ_0^2+2μ/μ_cd_0ϕ_0-μ^2/μ_c^2ϕ_0^2 +ν^2 / 3 C μ_c^2ϕ_0^2 -2 C σ_0^2 -√(3)ν/ 3 μ_cϕ_0σ_0)-d_0 η, 
    
    ϕ_1τ   =   √(6)/2μ_cψ_1-ϕ_0 η,
    
    ψ_1τ   =   √(6)/2(- μ d_1 -μ_cϕ_1 +μ^2/μ_cϕ_1 - ν/√(3)σ_1 -ν^2/ 3 C μ_cϕ_1)
    +(-d_0ψ_0+μ/μ_cϕ_0ψ_0)-ψ_0 η، 
    σ_1τ   =   (-d_0σ_0 +μ/μ_cσ_0ϕ_0) -σ_0 η.

Solving the partial differential equations for the Zeroth order terms yields



    d_0(τ,η)    =    B(η), 
    ϕ_0(τ,η)    =    A(η)cos(λτ-Λ(η))-μ_c/(μ_c ν^2 -μ^2)(μ B(η)+ν/√(3) S(η)),
    ψ_0(τ,η)    =    -√(6)λ/3μ_cA(η)sin(λτ-Λ(η)),
    σ_0(τ,η)     =    S(η).

where λ=√(6)/2√(μ_c ν^2-μ^2) and A(η), B(η), S(η) and Λ(η) are as yet undetermined functions of the slow time η.  Solving the partial differential equations for the first order terms, and restricting ourselves to only bounded solutions, determines a set of ordinary differential equations for the unknown functions A(η), B(η), S(η) and Λ(η).  If we replace B(η) with the linear combination 
    B̅(η)=μ_cν^2B(η)+μν/√(3)S(η)
 then the resulting set of differential equations become

    A_η   =   -1/21/μ_c ν^2-μ^2 A B̅, 
    B̅_η   =   -1/3(1/μ_cν^2-μ^2B̅^2 + μ_cν^2(μ_cν^2-μ^2)/2μ_c^2A^2 + 2μ_c C S^2 ),
    
    S_η   =   -            1/μ_c ν^2-μ^2 S B̅, 
    Λ_η   =    0.


Therefore, in terms of the original variables the first term of the perturbative solution is

    D(τ)    =   ϵ[ B(η) - μ/μ_cA(η)cos(λτ-Λ(η))
       +μ/(μ_cν^2-μ^2)(μ B(η)+ν/√(3) S(η)) ], 
    
    Φ(τ)    =   ϵ[ A(η)cos(λτ-Λ(η))
        -μ_c/(μ_cν^2-μ^2)(μ B(η)+ν/√(3) S(η)) ],
    
    Ψ(τ)    =   ϵ[-√(μ_c ν^2-μ^2)/μ_cA(η)sin(λτ-Λ(η))],
    
    Σ( τ)    =   ϵ[ S(η)+ν/√(3) C μ_c   A(η)cos(λτ-Λ(η))
        -ν/√(3)C(μ_cν^2-μ^2)(μ B(η)+ν/√(3) S(η))],

where the functions A(η), B(η), S(η), and Λ(η) satisfy the differential equations (<ref>), and due to (<ref>) are bounded by

    B̅(η)^2 ≥(μ_cν^2-μ^2/μ_c^2)  A(η)^2  + μ_c^2/μ_cν^2CS(η)^2,

where we note that if B̅(η)→ 0 then we also have A(η)→ 0 and S(η)→ 0 which then also implies that B(η)→ 0.

We are interested in determining the asymptotic behaviour as τ→∞.  We observe that the phase shift Λ(η) is a constant and has no effect on the future dynamics.  The fast time τ essentially describes the oscillations of the scalar field, which to first order in ϵ has a period of T=2π/λ.  We note that the period of these oscillations T ∼ 1/√(μ_c ν^2-μ^2), gets longer as the strength of the coupling parameter μ is increased.

We also observe that the amplitude of the oscillations A(η), and the vertical shift B(η) and the shear term S(η) are functions of the slow time η and consequently the amplitude, vertical shift, and shear term drift slowly in comparison to the oscillatory changes. For initial values of B(η)>0 we see that A(η),B(η),S(η) → 0 as η→∞.  That is, the amplitude of the oscillations, the vertical shift, and the shear term all slowly decrease to zero, indicating that the point P_0 is stable when μ< μ_c ν.




  §.§.§ Stability of Equilibrium Points in ^-∪^0

Unfortunately, while  we have an autonomous system of differential equations defined on a compact set, the system is not differentiable  at any points in the invariant set ^-∪^0. In order to determine the local behaviour at these  points, we replace variable D with

    T=1/√(1+ θ^2)= √(1-D^2).

The evolution  equations (<ref>)-(<ref>) become

    T^'   =    -T √(1-T^2)𝒳,  
    Φ^'   =    m TΨ-√(1-T^2)Φ𝒳 ,   
    Ψ^'   =    -√(1-T^2)Ψ(1+𝒳)- T (  m Φ+√(3/2)μ√(1-T^2)+ ν/√(2)Σ), 
    Σ^'   =    - Σ√(1-T^2)(1+  𝒳) + ν/√(2)  C( Ψ T +√(1-T^2)Φ/m),

with

    𝒳=-1/3(1-T^2)-2/3Ψ^2+1/3Φ^2 -2/3  C Σ^2 -√(3/2)μ  TΨ +√(3/2)ν/mΣΦ.

The value D=1 for the equilibrium points P_1 P_2, P_3, P_4, P_5 and C^* is simply replaced with T=0. With this transformation we are able to locally determine the qualitative behaviour of each the equilibrium points.

The eigenvalues of the linearization at the equilibrium points P_1, P_2, and C^* are 0, 1,1, 4/3 which implies these points are unstable and are sources. The zero eigenvalue indicates the non-isolated nature of the circle of equilibria C^*.

The eigenvalues of the linearization at the equilibrium point P_5  is 1/3, -2/3, 1/3, -2/3 which implies that  P_5 is  generally saddle (unstable). Further, the eigen-directions that span the T=0  invariant set, are associated with one positive and two negative eigenvalues. Therefore this equilibrium point is a saddle within the  T=0 invariant set.

The eigenvalues of the linearization at the equilibrium points P_3 and P_4 are 0, -1,-1, -2/3 which implies that we cannot  determine the general behaviour of these points without resorting to additional analysis. However, the eigen-directions that span the T=0 invariant set, are associated with the three negative eigenvalues. Therefore, these equilibrium points are sinks in the T=0 set. One method to complete the determination of the general behaviour near P_3 and P_4 is to calculate the center manifold <cit.>.

The center manifold for P_3 can be parameterized as

    T    =    T 
    Φ   =   μ_c/μ_cν-μ/μ_cν(1+3/4(μ+μ_cν)^2)T^2+O(T^4),
    
    Ψ   =     -√(6)/2(μ+μ_cν)T+√(6)/2μ_cν(1/2+3/4(μ+μ_cν)^2)T^3+O(T^4),
    
    Σ    =    ν/√(3)C1/μ_cν-ν/√(3)C1/μ_cν(1+3/4(μ+μ_cν)^2)T^2+ O(T^4)

The leading order term of the dynamical system restricted to the center manifold reduces to

    T'= 3/2μ_cν(μ+μ_cν) T^3.

Since T'>0 for T>0, P_3 is unstable along its center manifold. Therefore P_3 is a saddle within the full four dimensional phase space.


The center manifold for P_4 can be parameterized as

    T    =    T 
    Φ   =    -μ_c/μ_cν+μ/μ_cν(1+3/4(μ-μ_cν)^2)T^2+O(T^4),
    
    Ψ   =     -√(6)/2(μ-μ_cν)T-√(6)/2μ_cν(1/2+3/4(μ-μ_cν)^2)T^3+O(T^4),
    
    Σ    =     -ν/√(3)C1/μ_cν+ν/√(3)C1/μ_cν(1+3/4(μ-μ_cν)^2)T^2+ O(T^4)

The leading order term of the dynamical system restricted to the center manifold reduces to

    T'= -3/2μ_cν(μ-μ_cν) T^3.

If μ > μ_c ν then T'<0 and P_4 is stable along its center manifold. It is therefore a sink when μ > μ_c ν, in the full four dimensional phase space and a saddle otherwise.




  §.§.§ The Bifurcation Value


If μ=μ_cν then there is a line of equilibria given by (D,Φ,Ψ,Σ)=(s,-μ_c/μ_cνs,0,-ν/√(3)C1/μ_cνs) where 0<s<1 having endpoints P_0 and P_4.  The eigenvalues of points on this line of equilibria are 0, -s,-s, -2/3s which implies that points on this line are attractors when it exists. We observe that as the parameter μ increases towards its bifurcation value, μ_cν, the stability of the point P_0 is transferred to the point P_4 via this line of equilibria.



  §.§.§ Heteroclinic Sequences


Very often one is not only interested in the past and future behaviour of a system of differential equations, but one is also interested in the intermediate behaviour of the system.  One technique to analyze the intermediate behaviour is to describe the heteroclinic sequences that are possible <cit.>.  We note that for every heteroclinic sequence there exists a set of orbits that are arbitrarily close to that sequence.  Figure (<ref>) describe the possible heteroclinic sequences.  Again we see how the intermediate behaviour changes as the parameter μ changes.






 §.§ Inflation and Accelerated Expansion

It follows that the deceleration parameter (<ref>) can also be expressed in terms of the normalized bounded variables in this case as follows;

    q= -1/D^2(-2 C Σ^2-2 Ψ^2 +Φ^2-3√(3/2)μ√(1-D^2)Ψ+3 ν/√(2)ΦΣ).

The sign of the deceleration parameter indicates the nature of the expansionary evolution. If q > 0, then the cosmological expansion is decelerating, while negative values of q indicate an accelerating or inflationary dynamics. See Table (<ref>) for a summary of the sign of q for each equilibrium point.





 §.§ Numerical Analysis


It is constructive to illustrate a few numerical solutions for the three different regimes of future asymptotic behaviour, μ<μ_cν [see Figure (<ref>)], μ=μ_cν [see Figure (<ref>)], and μ>μ_cν [see Figure (<ref>)]. In each case the integrations are done in the full 4-dimensional phase space. The initial conditions are selected to show different past and future asymptotic behaviours and are the same as those used in the isotropic case, in that here we initially set Σ(0)=0.  We also do not show any phase portraits in this case as they are not as illustrative in higher dimensions as in the isotropic case.











 §.§ Observations


In the spatially homogeneous and isotropic case analyzed here we observe that the past dynamics are independent of the strength of the coupling parameters μ and ν.  We find in the zero curvature (Bianchi type I) and the negative curvature models (Bianchi type V) that the past asymptotic state is one which the anisotropy is non-trivial.  The past solution for both the zero and negative curvature models is the  Jacobs’ Bianchi type I non-vacuum massless scalar field solution <cit.>.

We also observe that the future asymptotic state depends on the strength of the coupling parameter μ and ν.  For weak coupling of the Aether field to the scalar field, i.e., μ < μ_c ν, the dynamics are similar to but not the same as that found when there is no coupling of the Aether field to the scalar field, i.e., when μ=0 and ν=0.  Having ν>0 drives the system towards intermediate states, P_3 and P_4, that are anisotropic in nature.  If μ < μ_c ν then P_0 is the stable attractor: orbits oscillate and slowly decay in amplitude towards this final isotropic non-inflationary asymptotic state.

Similar to the isotropic case analyzed in Section <ref>, for strong coupling of the Aether field to the scalar field, μ > μ_c ν, the dynamics are very different.  When μ>μ_cν the scalar field does not come to rest at the minimum of the potential: the strength of the Aether interaction forces a different final equilibrium state. If μ > μ_c ν then we find that the stable equilibrium point in these models changes from the isotropic point P_0 to the anisotropic point P_4 if ν>0 which is isotropic if ν=0.  In this case, we find that the square of the scalar field, the square of the shear scalar, and the square of the expansion scalar scale together as

    m^2/2ϕ^2    ∼   1/3μ_2^2/μ_cν^2(1+θ^2)
    σ^2    ∼   1/9ν^2/C^2μ_cν^2(1+θ^2)

and consequently grow without bound.  We also observe that all orbits (excepting for the exceptional orbits) experience some period of accelerating expansion as they evolve to their final asymptotic state which is consistent with the isotropic case. The fundamental difference in the anisotropic case when μ > μ_c ν is that the future asymptotic state need not be isotropic.

Bianchi type V models in the standard inflationary scenario in GR (C=1,μ=0,ν=0) isotropize as a rule. In GR, having interaction terms in the scalar field potential, equation (<ref>), in which (C=1, μ>μ_cν,ν>0), changes this rule to one in which the future asymptotic behaviour has accelerated expansion but is not isotropic. Similarly, in the Einstein Aether theory in the standard inflationary scenario (0<C<1,μ=0,ν=0), one observes once again that the models will isotropize to the future.  However, just as in GR,  if (0<C<1,μ>μ_cν,ν>0) then the models also have a future asymptotic behaviour which has accelerated expansion but is not isotropic.









§ DISCUSSION
 



 §.§ Slow Roll Inflation


In the standard slow roll inflationary scenario in GR, inflation occurs at intermediate times during a period of slow roll (in which ϕ̈≪θϕ̇ and ϕ̇^2 ≪θ^2) and where the anisotropy, if present, is insignificant when compared to the expansion.  However, the existence of a non-trivial coupling of the Aether field to the scalar field changes this scenario, and in particular if ν>0, then there can be a significant departure from the standard scenario.

To find the slow roll inflationary attractor with scalar field/aether field coupling, we start off with all the same assumptions as above except one. The anisotropy is not assumed to be insignificant, but changes in the anisotropy are assumed to be small during slow roll inflation.  With these assumptions the slow roll solutions take on the form

    θ   =   3/2μ_cν|ϕ|
    σ   =   ν/2cϕ
    ϕ̇   =    -μ_cν((ϕ)+μ/μ_cν)

We see that if μ<μ_cν, then the slow roll solution is stable as ϕ and ϕ̇ have different signs (same for σ and σ̇). The number of e-foldings N that can take place are:

    N    =   1/3∫_ϕ_i^ϕ_f θ/ϕ̇ dϕ
       =   1/4(1+(ϕ)μ/μ_cv)(ϕ_i^i2-ϕ_f^i2)

We note that if (ϕ)<0, then N can be made large for fixed initial and final endpoints by choosing μ≲μ_cν.
It appears that the existence of aether field/scalar field coupling terms of the nature studied here do not change the possibility of a period of slow-roll inflation at intermediate times.  We do note, however, that if ν>0 then the inflation is anisotropic in nature.  Indeed, slow roll inflation is possible even when the mass of the scalar field is zero by simply choosing the coupling parameters to satisfy μ < μ_cν|_m=0=ν/√(3C).  In some sense, the slow roll inflationary expansion in this scenario is driven by the coupling to the shear.




 §.§ Final Comments


We have investigated cosmological models in the Einstein-Aether theory in which scalar field matter is coupled to the aether through the scalar field potential. We have been especially interested in possible accelerated expansion and inflationary behaviour in a class of spatially homogeneous cosmological models. In particular, we have studied scalar field models in which the scalar-field potential depends on the time-like aether vector field through its expansion and shear.  We have observed that in the isotropic case, by choosing appropriate units, the dynamics are independent of the aether parameters.  The existence of the aether is to essentially re-normalize the gravitational constant G_c=G/(1+c_θ) in these cosmological settings. Further, in both the isotropic and anisotropic models studied here, we find that the past asymptotic state does not depend on the value of the aether parameters, c_i or on the value of the scalar field coupling parameters μ and ν.

However, in the anisotropic model the intermediate states and in some cases the final states do depend on a single combination of Aether parameters c_θ and c_σ through the parameter C.  Indeed, these intermediate and final states increasingly become more anisotropic as a result of decreasing the Aether parameter C from its maximum value C=1.  We also note that the future asymptotic states in both the isotropic and anisotropic models depend on the value of the scalar field/aether field coupling parameters μ and ν in the scalar field potential.  For sufficiently small values of the parameter μ, the both the isotropic and anisotropic models experience a period of slow-roll inflation at intermediate times, even when the scalar field is massless. For sufficiently large values of the parameter μ, the future asymptotic state changes to one which has accelerated expansion at late times.  Indeed, it is possible to have an accelerated expansion at late times even when the mass of the scalar field is zero, provided μ>0 in the isotropic case, and μ>ν/√(3C) in the anisotropic case.  In both of these isotropic and anisotropic models, the accelerated expansion at late times is a direct result of the scalar field/aether field coupling.  Further, it must be noted that in the anisotropic case, having a non-zero coupling parameter ν causes the future asymptotic state to be anisotropic.

The scalar field/aether field coupling parameters μ and ν in the scalar field potential modify the slow roll inflationary dynamics for a sufficiently small μ <μ_cν, which adds a driving force which can slow down or speed up (depending on the sign of scalar field initially) the slow roll inflation <cit.>.  In the anisotropic case there are further refinements to the slow roll regime (which can occur for a sufficiently small non-zero parameter ν in the potential).  Additionally, in the anisotropic case, the shear coupling causes the slow roll inflationary solution to be anisotropic in nature.

Recall, that a period of accelerated expansion is desirable at early intermediate times for inflationary purposes, but a period of accelerated expansion is also an attractive feature to have at late times to describe the effects of Dark Energy.  Here, in all cases (isotropic or anisotropic) or (zero curvature or negative curvature),  if μ is sufficiently small then there is a period of slow roll inflation at intermediate times, and if μ is sufficiently large, there will be accelerated expansion at late times.  Further, if ν>0, these statements are true even when the scalar field is massless.







This project is supported in part by the Atlantic Association for Research in the Mathematical Sciences through a Collaborative Research Grant.  We thank Theodore Kolokolnikov for his guidance on some technical points. BA would also like to thank the Government of Saudi Arabia for financial support. RvdH thanks the Department of Mathematics and Statistics at Dalhousie University for their kind hospitality.  AAC is supported by the Natural Sciences and Engineering Research Council of Canada.




10

Mattingly:2005re
D. Mattingly, Modern tests of Lorentz invariance,
  http://dx.doi.org/10.12942/lrr-2005-5Living Rev. Rel. 8
  (2005) 5, [https://arxiv.org/abs/gr-qc/0502097gr-qc/0502097].

Jacobson:2000xp
T. Jacobson and D. Mattingly, Gravity with a dynamical preferred
  frame, http://dx.doi.org/10.1103/PhysRevD.64.024028Phys.
  Rev. D64 (2001) 024028,
  [https://arxiv.org/abs/gr-qc/0007031gr-qc/0007031].

Liberati:2013xla
S. Liberati, Tests of Lorentz invariance: a 2013 update,
  http://dx.doi.org/10.1088/0264-9381/30/13/133001Class. Quant.
  Grav. 30 (2013) 133001, [https://arxiv.org/abs/1304.57951304.5795].

Jacobson:2008aj
T. Jacobson, Einstein-aether gravity: A Status report, PoS
  QG-PH (2007) 020, [https://arxiv.org/abs/0801.15470801.1547].

Zlosnik:2006zu
T. G. Zlosnik, P. G. Ferreira and G. D. Starkman, Modifying gravity with
  the Aether: An alternative to Dark Matter,
  http://dx.doi.org/10.1103/PhysRevD.75.044017Phys. Rev. D75 (2007) 044017, [https://arxiv.org/abs/astro-ph/0607411astro-ph/0607411].

Clifton:2011jh
T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, Modified Gravity
  and Cosmology,
  http://dx.doi.org/10.1016/j.physrep.2012.01.001Phys. Rept.
  513 (2012) 1–189, [https://arxiv.org/abs/1106.24761106.2476].

Jacobson:2004ts
T. Jacobson and D. Mattingly, Einstein-Aether waves,
  http://dx.doi.org/10.1103/PhysRevD.70.024003Phys. Rev. D70 (2004) 024003, [https://arxiv.org/abs/gr-qc/0402005gr-qc/0402005].

Jacobson:2010mxa
T. Jacobson, Extended Horava gravity and Einstein-aether theory,
  http://dx.doi.org/10.1103/PhysRevD.82.129901,
  10.1103/PhysRevD.81.101502Phys. Rev. D81 (2010) 101502,
  [https://arxiv.org/abs/1001.48231001.4823].

Jacobson:2010mxb
T. Jacobson, Erratum: Extended hořřava
  gravity and einstein-aether theory [phys. rev. d 81, 101502 (2010)],
  http://dx.doi.org/10.1103/PhysRevD.82.129901Phys. Rev. D 82 (Dec, 2010) 129901.

Garfinkle:2007bk
D. Garfinkle, C. Eling and T. Jacobson, Numerical simulations of
  gravitational collapse in Einstein-aether theory,
  http://dx.doi.org/10.1103/PhysRevD.76.024003Phys. Rev. D76 (2007) 024003, [https://arxiv.org/abs/gr-qc/0703093gr-qc/0703093].

Garfinkle:2011iw
D. Garfinkle and T. Jacobson, A positive energy theorem for
  Einstein-aether and Hořava gravity,
  http://dx.doi.org/10.1103/PhysRevLett.107.191102Phys. Rev.
  Lett. 107 (2011) 191102, [https://arxiv.org/abs/1108.18351108.1835].

Barrow:2012qy
J. D. Barrow, Some Inflationary Einstein-Aether Cosmologies,
  http://dx.doi.org/10.1103/PhysRevD.85.047503Phys. Rev. D85 (2012) 047503, [https://arxiv.org/abs/1201.28821201.2882].

Sandin:2012gq
P. Sandin, B. Alhulaimi and A. Coley, Stability of Einstein-Aether
  Cosmological Models,
  http://dx.doi.org/10.1103/PhysRevD.87.044031Phys. Rev. D87 (2013) 044031, [https://arxiv.org/abs/1211.44021211.4402].

Donnelly:2010cr
W. Donnelly and T. Jacobson, Coupling the inflaton to an expanding
  aether, http://dx.doi.org/10.1103/PhysRevD.82.064032Phys.
  Rev. D82 (2010) 064032, [https://arxiv.org/abs/1007.25941007.2594].

Solomon:2013iza
A. R. Solomon and J. D. Barrow, Inflationary Instabilities of
  Einstein-Aether Cosmology,
  http://dx.doi.org/10.1103/PhysRevD.89.024001Phys. Rev. D89 (2014) 024001, [https://arxiv.org/abs/1309.47781309.4778].

Solomon:2015hja
A. R. Solomon, Cosmology Beyond Einstein.
PhD thesis, Cambridge U., Cham, 2015.
https://arxiv.org/abs/1508.068591508.06859.
10.1007/978-3-319-46621-7.

Alhulaimi:2017
B. Alhulaimi, Einstein-aether Cosmological Scalar Field Models.
PhD thesis, Dalhousie University, Halifax, Nova Scotia, February,
  2017.

Kanno:2006ty
S. Kanno and J. Soda, Lorentz Violating Inflation,
  http://dx.doi.org/10.1103/PhysRevD.74.063505Phys. Rev. D74 (2006) 063505, [https://arxiv.org/abs/hep-th/0604192hep-th/0604192].

Olive:1989nu
K. A. Olive, Inflation,
  http://dx.doi.org/10.1016/0370-1573(90)90144-QPhys. Rept. 190 (1990) 307–403.

Linde:1987
A. D. Linde, Inflation and Quantum Cosmology,  in 300 Years of
  Gravity, pp. 604–630, 1987.

Riess:1998cb
Supernova Search Team collaboration, A. G. Riess et al.,
  Observational evidence from supernovae for an accelerating universe
  and a cosmological constant,
  http://dx.doi.org/10.1086/300499Astron. J. 116 (1998)
  1009–1038, [https://arxiv.org/abs/astro-ph/9805201astro-ph/9805201].

Perlmutter:1998np
Supernova Cosmology Project collaboration, S. Perlmutter et al.,
  Measurements of Omega and Lambda from 42 high redshift supernovae,
  http://dx.doi.org/10.1086/307221Astrophys. J. 517 (1999)
  565–586, [https://arxiv.org/abs/astro-ph/9812133astro-ph/9812133].

Coley:2015qqa
A. A. Coley, G. Leon, P. Sandin and J. Latta, Spherically symmetric
  Einstein-aether perfect fluid models,
  http://dx.doi.org/10.1088/1475-7516/2015/12/010JCAP 12
  (2015) 010, [https://arxiv.org/abs/1508.002761508.00276].

Latta:2016jix
J. Latta, G. Leon and A. Paliathanasis, Kantowski-Sachs Einstein-æther
  perfect fluid models,
  http://dx.doi.org/10.1088/1475-7516/2016/11/051JCAP 1611
  (2016) 051, [https://arxiv.org/abs/1606.085861606.08586].

Coley:2003mj
A. A. Coley, Dynamical systems and cosmology, vol. 291.
Kluwer, Dordrecht, Netherlands, 2003,
  http://dx.doi.org/10.1007/978-94-017-0327-710.1007/978-94-017-0327-7.

wainwright_ellis2005
J. Wainwright and G. Ellis, Dynamical Systems in Cosmology.
Cambridge University Press, 2005.

hinch1991
E. Hinch, Perturbation Methods.
Cambridge Texts in Applied Mathematics. Cambridge University Press,
  1991.

kevorkian2013
J. Kevorkian and J. Cole, Perturbation Methods in Applied Mathematics.
Applied Mathematical Sciences. Springer New York, 2013.

nayfeh2000
A. Nayfeh, Perturbation methods.
Physics textbook. John Wiley & Sons, 2000.

wiggins2003
S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and
  Chaos.
Texts in Applied Mathematics. Springer New York, 2003.

Alhulaimi:2013sha
B. Alhulaimi, A. Coley and P. Sandin, Anisotropic Einstein-aether
  cosmological models, http://dx.doi.org/10.1063/1.4802246J.
  Math. Phys. 54 (2013) 042503.

Carruthers:2010ii
I. Carruthers and T. Jacobson, Cosmic alignment of the aether,
  http://dx.doi.org/10.1103/PhysRevD.83.024034Phys. Rev. D83 (2011) 024034, [https://arxiv.org/abs/1011.64661011.6466].

Foster:2005dk
B. Z. Foster and T. Jacobson, Post-Newtonian parameters and constraints
  on Einstein-aether theory,
  http://dx.doi.org/10.1103/PhysRevD.73.064015Phys. Rev. D73 (2006) 064015, [https://arxiv.org/abs/gr-qc/0509083gr-qc/0509083].











§ APPENDIX:  CONSTRAINTS ON THE EINSTEIN AETHER PARAMETERS C_I
 

According to <cit.> the PPN parameters for Einstein's General Relativity and Einstein Aether theory are identical if

    c_2    =   -2c_1^2-c_1c_3+c_3^2/3c_1, 
    
    c_4    =    -c_3^2/c_1.

If one also assumes the squared speeds of massless modes relative to the Aether rest frame must be super-luminal to avoid vacuum Čerenkov radiation <cit.>, then two additional constraints must be satisfied

    0    ≤    c_1+c_3 ≤1/2, 
    
    0    ≤    c_1-c_3 ≤c_1+c_3/3[1-2(c_1+c_3)],

which when combined with (<ref>) are sufficient to show both positive energy modes and linear stability.  Recall that the c_i employed here are one-half of the values used in <cit.>. In terms of c_σ and c_θ used in this paper, assuming that Einstein Aether and GR are equivalent for weak fields, i.e., equation (<ref>), then these constraints (<ref>) become

    0    ≤    c_σ≤1/2
    c_σ/3c_σ-2   ≤    c_θ≤ 0

and in terms of c_σ and C

    0    ≤    c_σ≤1/2, 
    
    1-2c_σ   ≤    C ≤ 1-3/2c_σ.




Donnelly and Jacobson included the coupling parameter μ in their analysis of the parameter constraints and showed that these constraints are relaxed when μ>0 <cit.> and are automatically satisfied when the PPN parameters match those of GR.

