Dynamical Analysis of an Integrable Cubic Galileon Cosmological Model
Alex Giacomini, Sameerah Jamal, Genly Leon, Andronikos Paliathanasis, Joel Saavedra
aa r X i v : . [ g r- q c ] M a y Dynamical Analysis of an Integrable Cubic Galileon Cosmological Model
Alex Giacomini, ∗ Sameerah Jamal, † Genly Leon, ‡ Andronikos Paliathanasis,
1, 4, § and Joel Saavedra ¶ Instituto de Ciencias F´ısicas y Matem´aticas, Universidad Austral de Chile, Valdivia, Chile School of Mathematics and Centre for Differential Equations,Continuum Mechanics and Applications, University of the Witwatersrand, Johannesburg, South Africa Instituto de F´ısica, Pontificia Universidad Cat´olica de Valpara´ıso, Casilla 4950, Valpara´ıso, Chile Institute of Systems Science, Durban University of Technology,PO Box 1334, Durban 4000, Republic of South Africa
Recently a cubic Galileon cosmological model was derived by the assumption that the field equa-tions are invariant under the action of point transformations. The cubic Galileon model admitsa second conservation law which means that the field equations form an integrable system. Theanalysis of the critical points for this integrable model is the main subject of this work. To performthe analysis, we work on dimensionless variables different from that of the Hubble normalization.New critical points are derived while the gravitational effects which follow from the cubic term arestudied.
PACS numbers: 98.80.-k, 95.35.+d, 95.36.+xKeywords: Cosmology; Galileon; Dynamical system
1. INTRODUCTION
A theory which has drawn the attention of the scientific society in the last few years is the Galileon gravity [1, 2].It belongs to the modified theories of gravity in which a noncanonical scalar field is introduced and the field equationsare invariant under the Galilean transformation. The Action integral of the Galileon gravity belongs to the Horndeskitheories [3], which means that the gravitational theory is of second-order [4]. The vast applications of study forGalileons in the literature covers all the areas of gravitation physics from neutron stars, black holes, acceleration ofthe universe, for instance see [4–21] and references therein.In this work we are interested in the cosmological scenario and specifically in the so-called Galileon cosmology [22–24]. In cosmology, the Galileon field has been applied in order to explain various phases of the evolution of the universe[25–29]. Specifically the new terms in the gravitational Action integral can force the dynamics in a way such that themodel fits the observations. The mechanics can explain the inflation era [30–36] as also the late time-acceleration ofthe universe [37–42]. Last but not least the growth index of matter perturbations have been constrained in [43].As we mentioned in the previous paragraph, Galileon gravity belongs to the Hondenski theories and specificallythere is an infinite number of different models which can be constructed from a general Lagrangian. A simple modelis the cubic Galileon model [11–17] where the Action Integral is that of a canonical scalar field plus a new term whichhas a cubic derivative on the Galileon field. The theory can be seen as a first extension of the scalar field cosmology.Due to this cubic term, the nonlinearity and the complexity of the field equations is increased dramatically. Recentlyin [44] a cubic Galileon model was derived which admits an additional conservation law and where the field equationsformed an integrable dynamical system.Integrability is an important issue in all the areas of physics and mathematical sciences. The reason is that whilea dynamical system can be studied numerically, it is unknown if an actual solution which describes the “orbits”exist. The integrable cubic Galileon model admits special solutions which describes an ideal gas universe, that is,power-law scale factors. While this is similar to the solution of the canonical scalar field, we found that the power ofthe power-law solution is not strongly constrained by the Galileon field and that is because of the cubic term. Onthe other hand a special property of that model is that when the potential in the Action integral dominates, then thecubic term disappears which mean that the theory approach is that of a canonical scalar field. However as we shallsee from our analysis, the existence of the conservation law provides new dynamics which have not been investigatedpreviously. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected]
The scope of this work is to study the evolution of the integrable cubic Galileon cosmological model. For that weperform an analysis of the critical points. In particular, in Section 2 we briefly discuss the cubic Galileon cosmologyand we review the integrable case that was derived before in [44]. Section 3 includes the main material of our analysiswhere we rewrite the field equations in dimensionless variables. We define variables different from that of the H -normalizaton where we find that the dynamical system is not bounded. Because of the latter property the criticalpoints at the finite and the infinite regions are studied. At the finite region we find various critical points which candescribe the expansion history of our universe as also the matter dominated era. Appendices A, B, C, D, E and Finclude important mathematical material which justify our analysis. One important property of the integrable modelthat we study is that there is a limit in which the terms in the field equations (which follow from the cubic termof the Galileon Lagrangian) vanish and the model is then reduced to that of a canonical scalar field cosmologicalmodel. Hence, in order to study the effects of the cubic term in Section 4 we perform an asymptotic expansion of thesolution when the cubic term dominates the universe. In Section 5 we extend our analysis to the case where an extramatter term is included in the gravitational Action integral. Finally, we discuss our results and draw our conclusionsin Section 6.
2. CUBIC GALILEON COSMOLOGY
The cubic Galileon model is defined by the following Action integral S = Z d x √− g (cid:18) R − ∂ µ φ∂ µ φ − V ( φ ) − g ( φ ) ∂ µ φ∂ µ φ (cid:3) φ (cid:19) (1)which has various cosmological and gravitational applications.In the cosmological scenario of a homogeneous and isotropic universe with zero spatially curvature the line elementof the spacetime is that of the FLRW metric ds = − dt + a ( t ) (cid:0) dx + dy + dz (cid:1) , (2)where a ( t ) is the scale factor of the universe.Indeed for this line element, variation with respect to the metric tensor in (1) provides the gravitational fieldequations 3 H = ˙ φ (cid:16) − g ( φ ) H ˙ φ + g ′ ( φ ) ˙ φ (cid:17) + V ( φ ) , (3)and 2 ˙ H + ˙ φ (cid:16) g ′ ( φ ) ˙ φ − g ( φ ) H ˙ φ + g ( φ ) ¨ φ (cid:17) = 0 , (4)while variation with respect to the field φ , provides the (modified) “Klein-Gordon” equation¨ φ (cid:16) φ g ′ ( φ ) − H g ( φ ) ˙ φ + 1 (cid:17) + ˙ φ (cid:18)
12 ˙ φ g ′′ ( φ ) − g ( φ ) ˙ H − H g ( φ ) (cid:19) + 3 H ˙ φ + V ′ ( φ ) = 0 , (5)which describes the evolution of the field, and H = ˙ aa . Recall that we have assumed that the Galileon field inheritsthe symmetries of the spacetime, that is if K µ is an isometry of (2), i.e. [ K, g µν ] = 0, then φ inherits the symmetriesof the spacetime if and only if [ K, φ ] = 0. Therefore φ is only a function of the “ t ” parameter, i.e. φ ( x µ ) = φ ( t ).An equivalent way to write the field equations is by defining fluid components such as energy density and pressurewhich corresponds to the Galileon field. Indeed if we consider the energy density ρ G = ˙ φ (cid:16) − g ( φ ) H ˙ φ + g ′ ( φ ) ˙ φ (cid:17) + V ( φ ) , (6)and the pressure term p G = ˙ φ (cid:16) g ¨ φ + g ,φ ˙ φ (cid:17) − V ( φ ) , (7)the field equations take the form G µν = T ( G ) µν , where T ( G ) µν is the energy momentum tensor corresponding to theGalileon field and T ( G ) µν = ρ G u µ u ν + p G ( g µν + u µ u ν ) (8)in which u µ = δ µt is the normalized comoving observer ( u µ u µ = − T ( G ) µν ; ν = 0, that is,˙ ρ DE + 3 H ( ρ DE + p DE ) = 0 . (9)Last but not least the dark energy equation-of-state parameter is defined as follows. w DE ≡ p DE ρ DE = 13 H " ˙ φ (cid:16) g ¨ φ + g ,φ ˙ φ (cid:17) − V ( φ ) . (10)It can be shown that with a proper election of g ( φ ), the equation-of-state parameter w DE can realize the quintessencescenario, the phantom one and cross the phantom divide during the evolution, which is one of the advantages ofGalileon cosmology. In general the specific functions of g ( φ ) and V ( φ ) are unknown and for different functions therewill be a different evolution. Recently in [44] two unknown functions were specified by the requirement that the gravitational field equationsform an integrable dynamical system. Specifically the following functions were found V ( φ ) = V e − λφ and g ( φ ) = g e λφ . (11)There exists a symmetry vector which provides, with the use of Noether’s second theorem, the following conservationlaw for the field equations I = − (cid:18) a ˙ a − λ a ˙ φ (cid:19) + g e λφ a ˙ φ − λ g a e λφ ˙ a ˙ φ . (12)Because of the nonlinearity of the field equations the general solution cannot be written in a closed-form. Howeverfrom the symmetry vector invariant curves have been defined and by using the zero-th order invariants some powerlaw (singular solutions) have been derived. In particular the following solutions were obtained: a ( t ) = a t p , φ ( t ) = 2 λ ln ( φ t ) , g = λ (cid:0) − λ p (cid:1) p − φ , V = φ (cid:18) λ + p (3 p − (cid:19) (13)and a , ( t ) = a t , φ , ( t ) = ± √
63 ln( φ t ) , V = 0 , λ , = ±√ . (14)These solutions are special solutions since they exist for specific initial conditions. In order to study the generalevolution of the system we perform an analysis in the phase-space.A phase-space analysis for this cosmological model has been performed previously in [16], however the integrable casewith g ( φ ) and V ( φ ) given by the expressions (11) were excluded from [16]. Moreover there is a special observationin the integrable case in the sense that V ( φ ) g ( φ ) = const . The latter means that when V ( φ ) dominates, g ( φ )becomes very small and the cubic Galileon model reduces to that of a canonical scalar field which can mimic also thecosmological constant when V ( φ ) ≫ ˙ φ .In the following we write the field equations in new dimensionless variables and we perform our analysis.
3. EVOLUTION OF THE DYNAMICAL SYSTEM
From equation (3) one immediately sees that the Hubble function H ( t ) can cross the value H ( t ) = 0, from negativeto positive values, or vice-verca. This means that the standard H -normalizaton is not useful and new variables haveto be defined. We introduce the new variables x = ˙ φ p H + 1) , y = V e − λφ H + 1) , z = H √ H + 1 , (15)and the parameter α = g V so we obtain the three-dimensional dynamical system x ′ = f ( x ) f ( x ) , y ′ = f ( x ) f ( x ) , z ′ = f ( x ) f ( x ) (16)where x = ( x, y, z ) and functions f − f are defined by the following expressions f ( x, y, z ) = y (cid:16) x z + 3 xz (cid:0) − y + z − (cid:1) + √ λy (cid:17) + 6 α x (cid:16) λ x z − √ λx − √ λxz + 6 z (cid:17) + αx y (cid:16) λx z − √ x (cid:0) λ + 12 z + 3 (cid:1) − λxz (cid:0) y − z (cid:1) + 3 √ (cid:0) yz + y − z + z (cid:1)(cid:17) , (17a) f ( x, y, z ) = 2 y (cid:16) α λ x z − √ α λx (cid:0) z + 1 (cid:1) + 18 αx z (cid:0) λy + 2 αz (cid:1) − √ αx y (cid:0) λ + 3 z (cid:1) +3 x yz (cid:0) y − αλ (cid:0) y − (cid:0) z + 1 (cid:1)(cid:1)(cid:1) + √ xy (cid:0) αz (cid:0) y − z (cid:1) − λy (cid:1) + 3 y z (cid:0) z − y (cid:1)(cid:17) , (17b) f ( x, y, z ) = 3 (cid:0) z − (cid:1) (cid:16) y (cid:16) x + 2 αx (cid:16) √ z − λx (cid:17) + z (cid:17) + 4 α x (cid:16) λ x − √ λxz + 3 z (cid:17) +2 αxy (cid:16) λx − √ x z + 2 λxz − √ z (cid:17) − y (cid:17) , (17c) f ( x, y, z ) = 2 (cid:16) αx (cid:16) αx + 2 λxy − √ yz (cid:17) + y (cid:17) , (17d)and prime denotes the new derivative dfdτ ≡ f ′ ≡ ˙ f √ H +1 .Interestingly, the system (16) admits the first integral2 αx (cid:16) λx − √ z (cid:17) + y ( x − z ) + y = 0 (18)which is the Friedmann’s first equation and constrains the evolution of the solution. Thus, the dynamics are restrictedto a surface given by (18). For a fixed value of y , the first and last equations in (16) are invariant under the discretesymmetry ( x, z, τ ) → − ( x, z, τ ). Thus, the fixed points related by this discrete symmetry have the opposite dynamicalbehavior. By definition, y ≥ x = ˙ φ √ H , y = √ V ( φ ) √ H , z = g ( φ ) H ˙ φ. Since we have chosen here V ( φ ) and g ( φ ) such that g ( φ ) V ( φ ) = α we have the relations x = xz , y = √ yz , z = √ αxzy , and the extra relation αx − q y z = 0 . This implies that the fixed points A ± , B ± , C and D investigatedin detail in [16], do not exist in our scenario since the values of their coordinates ( x , y , z ) do not satisfy the aboveextra relation.Furthermore, some cosmological parameters with great physical significance are the effective equation of stateparameter w tot ≡ p tot ρ tot = w DE (because we set ρ m = 0) and the ‘deceleration parameter” q ≡ − − ˙ HH = 12 + 32 w tot . (19)Some conditions for the cosmological viability of the most general scalar-tensor theories has to be satisfied by extendedGalileon dark energy models; say the model must be free of ghosts and Laplacian instabilities [45–47]. In the specialcase of the action (1) (in units where κ ≡ πG = 1), we require for the avoidance of Laplacian instabilities associatedwith the scalar field propagation speed that [46] c S ≡ w H − w − w w + 9 w ≥
0; (20) w ≡ g ˙ φ + 2 H, w ≡ φ (cid:20)
12 + g ,φ ˙ φ − Hg ˙ φ (cid:21) − H . (21)Meanwhile, for the absence of ghosts it is required that Q S ≡ (4 w + 9 w )3 w > . (22)Finally, we have from the Eqs. (10), (21) and (22) that the phantom phase can be free of instabilities and thuscosmologically viable, as it was already shown for Galileon cosmology [46]. The fixed points/ fixed lines at the finite region of the system (16), and a summary of their stability conditions arepresented in table I. In table II we display several cosmological parameters for the fixed points at the finite region ofthe system (16). The discussion about the physical interpretation of these points and the points at infinity is left forsection 3.4TABLE I: Summary of the stability conditions of the fixed points at the finite region of the system (16). Where P ( x ) = 6(1 + αλ ) x + 3 √ α + λ ) x + 2(3 + λ ) x + √ λ . We used the acronyms “A.S.” for Asymptotically Stableand “A.U.” for Asymptotically Unstable. Label: Coordinates ( x, y, z ) Existence Eigenvalues Stability P : (0 , ,
0) Always Undetermined Numerical analysis P : (cid:16) − √ λ , , − (cid:17) Always 0 , − , − λ < −√ < λ < √ −√ < λ < λ > √ P : (cid:16) √ λ , , (cid:17) Always 0 , , λ < −√ < λ < √ −√ < λ < λ > √ P : (0 , , λ = 0 − , − , P : (0 , , − λ = 0 3 , , P ( x c ) : (cid:16) x c , q λx c + x c + 1 , − (cid:17) Where x c = 0 is a real root of P ( x )such that q λx c + x c + 1 ≥ P ( x c ) : (cid:16) − x c , q λx c + x c + 1 , (cid:17) Where x c = 0 is a real root of P ( x )such that q λx c + x c + 1 ≥ P : (0 , , −
1) Always − , − , − (cid:0) αλ − (cid:1) Sink for αλ > Saddle otherwise P : (0 , ,
1) Always 3 , , (cid:0) αλ − (cid:1) Source for αλ > Saddle otherwise P : (cid:18) − √ λ , , − (cid:19) α = 0 , λ = 0 (cid:0) − λ , (cid:0) − λ (cid:1) , − λ − (cid:1) sink for 0 < λ < saddle otherwise P : (cid:18) √ λ , , (cid:19) α = 0 , λ = 0 (cid:0) λ , − (cid:0) − λ (cid:1) , λ + 3 (cid:1) source for 0 < λ < saddle otherwise P : (cid:18) − √ λ , λ , − (cid:19) λ + 3 αλ − , , − saddle P : (cid:18) √ λ , λ , (cid:19) λ + 3 αλ − , , − saddle P ( z c ) : (cid:0) , z c , z c (cid:1) λ = 0 (0 , − z c , − z c ) A.S. for 0 < z c ≤ − ≤ z c < P ( z c ) : (cid:0) βz c , √ αβz c , z c (cid:1) λ = 0 , β = √ (cid:0) √ α − √ α − (cid:1) , α > q (0 , − z c , − z c ) A.S. for 0 < z c ≤ − ≤ z c < P ( z c ) : (cid:0) βz c , √ αβz c , z c (cid:1) λ = 0 , β = √ (cid:0) √ α + √ α − (cid:1) , α > q (0 , − z c , − z c ) A.S. for 0 < z c ≤ − ≤ z c < TABLE II: Cosmological parameters for the fixed points at the finite region of the system (16).
Label (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) Physical interpretation P (cid:16) − , , , ( a λ + b ( − λ )) τ b (3 b − a ) λ + O ( τ − ) , ( a λ + b ( − λ )) τ b (3 b − a ) λ + + O ( τ − ) (cid:17) a ( t ) ≈ ( b − b t ) b / (cid:0) a + O (cid:0) t − ) (cid:1)(cid:1) . a = − λ ( λ − α ( λ − )) ±√ q λ ( α ( αλ − λ +3 ) + λ ) λ ( α ( λ − ) − λ ) , b = , or The Galileon mimics radiation for b = .a = , b = , or The Galileon mimics matter for b = . b = a , a (cid:0) a λ (cid:0) α (cid:0) λ − (cid:1) − λ (cid:1) + 1 (cid:1) = 0 Powerlaw-solution for b = a = 0 . The Galileon mimics dust for4 a λ + b (cid:0) − λ (cid:1) = 0. c s < P , (cid:0) − , , λ , , (cid:1) The Galileon mimics dust. c s < P , (cid:0) − , − , , − , − (cid:1) de Sitter solution. c s < , Q S < P , ( x c ) see section 3.4 Accelerated solution for x c < , λ < − x c +2 √ x c , or x c > , λ > − x c +2 √ x c P , (cid:0) − , − , , , (cid:1) The Galileon mimics dust. c s < , Q S < P , (cid:18) − ( λ − ) , λ − , λ , , (cid:0) λ + 2 (cid:1)(cid:19) The Galileon mimics dust. c s ≥ , Q S > < λ < P , (cid:18) λ − λ +1) ( λ − ) λ − λ − , − ( λ − λ − )( λ − ) , , − λ , (cid:19) The Galileon mimics dust c s ≥ , Q S > ≤ λ < (cid:0) √ (cid:1) ≈ . P ( z c ) (cid:0) − , − , , − , − (cid:1) de Sitter solution.˙ φ ≈ , a ( t ) ≈ a e tzc √ − zc , z c = ± q V V .c s < , Q S < P , ( z c ) see section 3.4 de Sitter solution.∆ φ ≈ √ βz c t √ − z c , a ( t ) ≈ a e tzc √ − zc , z c = ± √ V √ √ αβ + V . We proceed with the determination of critical points at the infinity.
Because the phase space of the system is unbounded we introduce the Poincar`e compactification and a new timederivative f ′ X = x p x + y , Y = y p x + y , f ′ → (1 − X − Y ) f ′ . (23)The dynamics on the “cylinder at infinity” can be obtained by setting X = cos θ, Y = sin θ ; the dynamics in thecoordinates ( θ, z ) is governed by the equations θ ′ = h ( θ, z ) := λ z sin( θ ) cos ( θ ) , (24a) z ′ = h ( θ, z ) := λ (cid:0) z − (cid:1) cos ( θ ) . (24b)We linearize around a given fixed point on the “cylinder at infinity” by introducing X = cos θ − ε , Y = sin θ − ε ,with ε ≪ , ε ≪
1. Notice that 1 − X − Y ≃ ε cos θ + 2 ε sin θ , so to examine the stability of the fixed points atthe cylinder, and from the interior of it, we have to estimate how r = ε cos θ + ε sin θ evolves, not only the stabilityin the plane ( θ, z ). We obtain for ε ≪ , ε ≪ r ′ = r (cid:2) λ z cos ( θ )(cos(2 θ ) − (cid:3) . (25)TABLE III: Summary of the stability conditions of the fixed points at infinity of the system (16). Label: Coordinates ( θ, z ) Coordinates (
X, Y, z ) r ′ /r ( λ , λ ) Stability Q : (cid:0) π , z c (cid:1) (0 , , z c ) 0 (0 ,
0) Nonhyperbolic Q : (0 , −
1) (1 , , −
1) 2 λ (cid:0) − λ , − λ (cid:1) Saddle Q : ( π, −
1) ( − , , −
1) 2 λ (cid:0) − λ , − λ (cid:1) Saddle Q : (0 ,
1) (1 , , − λ (cid:0) λ , λ (cid:1) Saddle Q : ( π,
1) ( − , , − λ (cid:0) λ , λ (cid:1) Saddle
The stability condition of a fixed point along r is then r ′ /r <
0. The full stability of the above fixed points issummarized in the table III.
Let us complete our analysis by performing some numerical simulations. Specifically we choose the constants of themodel to satisfy the conditions g = − λ (cid:0) λ p − (cid:1) p − φ , V = φ (cid:18) λ + p (3 p − (cid:19) , hence, α = − λ ( λ + p (3 p − )( λ p − ) p − . Moreover we impose the condition p > , which guarantees the stability of the perturbation of the scaling solution[44]. Because we have chosen α ≥
0, this leads to the “allowed” region on the parameter space defined byi) λ < −√ , p ≥ q λ − λ + , orii) λ = −√ , p > , oriii) −√ < λ < , p ≥ λ , oriv) 0 < λ < √ , < p ≤ λ , orv) λ > √ , < p ≤ q λ − λ + ,as displayed in figure 1. FIG. 1: “Allowed” region on the parameter space.FIG. 2: Poincar`e projection of the system (16) on the invariant set z = −
1. The green dot corresponds to the points P ( x c ). In the special case λ = 6 , p = we have 1 + αλ = 0, thus, the polynomial P ( x ) is quadratic and there areonly two roots of P ( x ) = 0. The blue contour is defined by f ( x, y, −
1) = 0. As shown in the figures this line issingular and attracts some orbits. The brown solid line corresponds to intersection of the invariant surface2 αx (cid:0) λx − √ z (cid:1) + y ( x − z ) + y = 0 and the invariant set z = −
1. In the top figures, P attracts some orbits, butothers are attracted by one the of green points associated to P ( x c ). P is not the attractor of the whole phase spacesince αλ < . In the bottom figures P is the attractor, not just in this invariant set but in the whole phase spacesince αλ > . The red thick dashed line denotes the local center manifold of P .In figure 2, it is presented a Poincar`e projection of the system (16) on the invariant set z = −
1. The green dotscorrespond to the points P ( x c ) (that we solved numerically). In the special case λ = 6 , p = we have 1 + αλ = 0,thus, the polynomial P ( x ) is quadratic and there are only two roots of P ( x ) = 0. The blue contour is definedby f ( x, y, −
1) = 0. As shown in the figures this line is singular, and attracts some orbits. The brown solid linecorresponds to intersection of the invariant surface 2 αx (cid:0) λx − √ z (cid:1) + y ( x − z ) + y = 0 and the invariant set z = −
1. In the top figures, P attracts some orbits, but others are attracted by one the of green points associated to P ( x c ). P is not the attractor of the whole phase space since αλ < . In the bottom figures P is the attractor, notjust in this invariant set but in the whole phase space since αλ > . In this section we discuss the stability conditions, cosmological properties and physical meaning of the (lines of)fixed points in both finite and infinite regions. • P : (0 , ,
0) always exists. To analyze the stability we resort to numerical examination.In the Appendix B we proved, using normal forms calculations, that the fixed point P corresponds to thecosmological solution (B15). Moreover, in order to improve the range and accuracy we calculate the diagonalFIG. 3: Poincar`e projection of the system (16) on the invariant surface (18) in the coordinates( X, Y, z ) = (cid:18) x √ x + y , y √ x + y , z (cid:19) . The red continuous line corresponds to the exact solution (28). The origin(represented by a blue dot) attracts this line. The vector field (16) is projected onto the surface.first order Pade approximants [1 / ˙ φ ( t ) , [1 / H ( t ) , [1 / φ ( t ) , around t = ∞ . This yields the following approximate expressions˙ φ ( t ) ≈ λ (2 a − b ) (cid:18) a − b − a t + 2 a + 3 b t − b + 3 b ln tt (cid:19) , (26a) H ( t ) ≈ b (cid:18) − tt − b − b t (cid:19) , (26b) φ ( t ) ≈ ln (cid:16) V αλ (2 a − b ) +2 λ (2 a − b ) +3 b (9 b − (cid:17) λ + 2 ln tλ , (26c)while the scale factor is calculated to be a ( t ) ≈ a e b t (cid:18) t (cid:19) − b t ( b − b t ) b / = ( b − b t ) b / (cid:0) a + O (cid:0) t − ) (cid:1)(cid:1) . (27)Due to the new conservation law (12), the allowed values of the constants a , a , b , b are:0i) a = − λ (cid:0) λ − α (cid:0) λ − (cid:1)(cid:1) ± √ p λ ( α (6 αλ − λ + 3) + λ )6 λ ( α (2 λ − − λ ) , b = 29 ,a (6 − b ) + 36 a − , a −
23 + I a = 0 , orii) a = 13 , b = 49 , αλ − αλ − λ = 0 , I a + 4 a − b = 0 , oriii) a (cid:0) a λ (cid:0) α (cid:0) λ − (cid:1) − λ (cid:1) + 1 (cid:1) = 0 , b = 43 a , b = 43 a , I = 0 . It is interesting to note that the power-law solution (13) satisfies the condition( x ( τ ) , y ( τ ) , z ( τ )) = q λ p p + t , V φ ( p + t ) , p p p + t . (28)Additionally, after the substitution of the functional forms of ( x ( τ ) , y ( τ ) , z ( τ )) in (28), and the substitutionof g = − λ ( λ p − ) p − φ , V = φ (cid:0) λ + p (3 p − (cid:1) it follows that the restriction (28) is satisfied for all the values of t .There exists a relation between τ and t given by τ = p p + t − p ln p (cid:16)p p + t + p (cid:17) t , (29)such that t → ∞ implies τ → ∞ . Thus, as τ → ∞ , this power-law solution approaches the P as t → ∞ . Forlarge τ we can invert to have t = 14 (cid:16)p p + 4( p ln( p ) + τ ) + 2 p ln( p ) + 2 τ (cid:17) = τ + p ln( p ) + p τ − p ln( p )2 τ + O (cid:0) τ − (cid:1) . (30)Thus, we can take as approximation for large τ : x = q λτ − q p ln( p ) λτ + O (cid:0) τ − (cid:1) , y = V φ τ + O (cid:0) τ − (cid:1) , z = pτ − p ln( p ) τ + O (cid:0) τ − (cid:1) . Furthermore λφ = (2 ln( φ ) + 2 ln τ ) + 2 p ln( p ) τ − p (cid:0) ln ( p ) − (cid:1) τ + O (cid:0) τ − (cid:1) , ˙ φ = 2 λτ − p ln( p )) λτ + O (cid:0) τ − (cid:1) . These features are represented in Figure 3. There it is shown the Poincar`e projection of the system (16) onthe invariant surface (18) in the variables (
X, Y, z ) = (cid:18) x √ x + y , y √ x + y , z (cid:19) . The red continuous linecorresponds to the exact solution (28). The origin (represented by a blue dot) attracts this line. The vectorfield (16) is projected onto the surface.The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:16) − , , , ( a λ + b ( − λ )) τ b (3 b − a ) λ + O ( τ − ) , ( a λ + b ( − λ )) τ b (3 b − a ) λ + + O ( τ − ) (cid:17) . The scale factor satisfy a ( t ) ≈ ( b − b t ) b / (cid:0) a + O (cid:0) t − (cid:1)(cid:1) . The Galileon mimics radiation for b = . It mimics matter for b = . It is a Powerlaw-solution for b = a = 0 . Finally, the Galileon mimicsdust for 4 a λ + b (cid:0) − λ (cid:1) = 0. Furthermore c s <
0. This point has not been obtained previously in [16] andin [17], since in these works the authors used H -normalization, which fails obviously when H = 0.1 • P : (cid:16) − √ λ , , − (cid:17) always exists. The eigenvalues are 0 , − , −
6, so, the points are nonhyperbolic. Using theCenter Manifold Theory we have proven that P is stable for λ < −√ < λ < √
6, and saddle for −√ < λ < λ > √
6. (see appendix C). This point represents kinetic dominated solutions with H → −∞ .The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:0) − , , λ , , (cid:1) . The Galileon mimics dust. Furthermore c s < • P : (cid:16) √ λ , , (cid:17) always exists. The eigenvalues are 0 , ,
6, so, the points are nonhyperbolic. This point has theopposite dynamical behavior of P . Using the Center Manifold Theory (in a similar way as in the appendix Cwe can prove that P is unstable for λ < −√ < λ < √
6, and saddle for −√ < λ < λ > √ H → + ∞ . The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:0) − , , λ , , (cid:1) . The Galileon mimics dust. Furthermore c s <
0. So, the Laplacian instabilities associated withthe scalar field propagation speed cannot be avoided for this solution [46]. • The fixed point P : (0 , ,
1) exists if λ = 0. As before, it corresponds to the special case V = V , and thecoupling function becomes constant too, g = g ; it is a de Sitter solution, but now H → + ∞ . The eigenvaluesare − , − ,
0, so it is nonhyperbolic. Using the Center Manifold Theory we obtain that it is asymptoticallystable (for details see Appendix D). Furthermore, perturbations from the equilibrium grow or decay algebraicallyin time, not exponentially as in the usual linear stability analysis. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:0) − , − , , − , − (cid:1) . As it can be seen, c s < , Q S <
0. Thus, this solution suffers from Laplacian instabilitiesand the presence of ghosts. • The fixed point P : (0 , , −
1) exists if λ = 0. The eigenvalues are 3 , ,
0, so it is nonhyperbolic. To analyze theirstability we resort to numerical examination or use the Center Manifold Theory. This fixed point correspondsto a de Sitter solution driven by a cosmological constant, since V = V , and the coupling function becomesconstant too, g = g (although, H → −∞ , so, it is not cosmological viable). This point has the oppositedynamical behavior of P , thus it is asymptotically unstable. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:0) − , − , , − , − (cid:1) . The associated cosmological solution corresponds to a de Sitter solution. However, c s < , Q S <
0. Thus, this solution suffers from Laplacian instabilities and the presence of ghosts. • For each choice ǫ = ±
1, there are 1, 2 or 3 isolated fixed points of the form P , ( x c ) : (cid:16) ǫx c , q λx c + x c + 1 , − ǫ (cid:17) ,where x c are the nonzero real roots of the polynomial P ( x ) = 6(1 + αλ ) x + 3 √ α + λ ) x + 2(3 + λ ) x + √ λ ,satisfying q λx c + x c + 1 ≥
0. To analyze their stability we resort to numerical examination.The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P , ( x c ) are c s ∗ = − P ( x c ) / ( Q ( x c ) Q ( x c )), P ( x c ) = 144 α x c + 24 α x c y c (cid:0) λx c + 2 √ (cid:1) + 6 α x c y c (cid:0) (cid:0) λ − (cid:1) x c + 6 √ λx c + 18 y c + 11 (cid:1) +2 αx c y c (cid:0) λx c + λx c (21 y c −
8) + 12 √ x c + √ y c − (cid:1) + y c (cid:0) − x c + 3 y c − (cid:1) , Q ( x c ) = 3 (cid:0) α x c + 2 αx c y c (cid:0) λx c + √ (cid:1) + y c (cid:1) , and Q ( x c ) = 24 α x c − αx c y c (cid:0) √ − λx c (cid:1) + (cid:0) x c − (cid:1) y c , where y c = q λx c + x c + 1 ≥ Q ∗ S = ( ( α +4 αλ +1 ) x c +4 √ ( α ( λ − ) +2 λ ) x c +(8 λ (3 α + λ )+15) x c +2 √ λ − α ) x c − ( λ +1 ) x c − √ λx c − )( − √ αx c + √ λx c +3 x c +3 ) , Ω ∗ DE = x c ( x c ( αλ +1) x c +3 √ α + λ ) x c +2 λ +9 ) +2 √ λ ) +3 √ λx c +3 x c +3 , ω ∗ DE = − x c (cid:0) x c + √ λ (cid:1) −
1, and q ∗ = (cid:0) − x c (cid:0) x c + √ λ (cid:1) − (cid:1) .The fixed points P , ( x c ) represents accelerated solution for x c < , λ < − x c +2 √ x c , or x c > , λ > − x c +2 √ x c .Since the above analytical expressions for c s and Q S are quite complicated, we have resorted to numericalinvestigation. To represent the regions of physical interest, we proceed in the following way. Recall that thepolynomial P ( x ) = 6(1 + αλ ) x + 3 √ α + λ ) x + 2(3 + λ ) x + √ λ has a discrete number of roots x c (1, 2,or 3 depending on the parameters α and λ ). Since this polynomial is linear in α , we have for each value of x c ( x c = 0 , x c = − √ λ ) the relation α := α ( x c , λ ) = − ( x c + x c ) + √ λ ( x c +1 ) +2 λ x c x c ( λx c + √ ) , α >
0. So, we can represent theregions of physical interests on the parameter space ( x c , λ ), rather than in the plane ( α, λ ) (to avoid solving agenerically third order polynomial in x c using Cardano’s formulas, with the subsequent numerical errors issue).The above procedure leads to the conditions c s ≥ , Q S > c s < , Q S ≤
0, for P , ( x c ) as displayed inFigure 4. In order to cover all the possible values of the parameters we have made the representation in the(compact) variables (cid:18) x c √ x c + λ , λ √ x c + λ (cid:19) .2FIG. 4: Parameter space that leads to the conditions c s ≥ , Q S > c s < , Q S ≤
0, for P , ( x c ). • The fixed point P : (0 , , −
1) always exists. The eigenvalues are − , − , − (cid:0) αλ − (cid:1) . Thus, it is a sink for αλ > or saddle otherwise. This solution is dominated by the Hubble scalar with H → −∞ . The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:0) − , − , , , (cid:1) . The Galileon mimics dust. However, the fluid satisfies theconditions c s < , Q S <
0. Thus, this solution suffers from Laplacian instabilities and the presence of ghosts. • The fixed point P : (0 , ,
1) always exists. The eigenvalues are 3 , , (cid:0) αλ − (cid:1) . Thus, it is a source for αλ > or a saddle otherwise. This solution is dominated by the Hubble scalar with H → + ∞ . The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:0) − , − , , , (cid:1) . The Galileon mimics dust. However, the fluid satisfies theconditions c s < , Q S <
0. Thus, this solution suffers from Laplacian instabilities and the presence of ghosts. • The fixed point P : (cid:18) − √ λ , , − (cid:19) exists for α = 0 , λ = 0. The eigenvalues are (cid:0) − λ , (cid:0) − λ (cid:1) , − λ − (cid:1) .Thus, it is a sink for 0 < λ < and a saddle otherwise. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:16) − λ − , λ − , λ , , (cid:0) λ + 2 (cid:1)(cid:17) . The Galileon mimics dust. Furthermore, c s ≥ , Q S > < λ <
2. In this region of the parameter space, the cosmological solution is free of Laplacian instabilities and ghosts. • The fixed point P : (cid:18) √ λ , , (cid:19) exists for α = 0 , λ = 0. The eigenvalues are (cid:0) λ , − (cid:0) − λ (cid:1) , λ + 3 (cid:1) .It is a source for 0 < λ < or a saddle otherwise. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:16) − λ − , λ − , λ , , (cid:0) λ + 2 (cid:1)(cid:17) . The Galileon mimics dust. Furthermore, c s ≥ , Q S > < λ <
2. In this region of the parameter space, the cosmological solution is free of Laplacian instabili-ties and ghosts. • The fixed point P : (cid:18) − √ λ , λ , − (cid:19) exists for λ + 3 αλ − , , − . It is asaddle. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:18) λ − λ +1) ( λ − ) λ − λ − , − ( λ − λ − ) ( λ − , , − λ , (cid:19) . TheGalileon mimics dust. Furthermore, c s ≥ , Q S > ≤ λ < (cid:0) √ (cid:1) ≈ . • The fixed point P : (cid:18) √ λ , λ , (cid:19) exists for λ + 3 αλ − , , − . It is a saddle.The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P are (cid:18) λ − λ +1) ( λ − ) λ − λ − , − ( λ − λ − ) ( λ − , , − λ , (cid:19) . The Galileonmimics dust. Furthermore, c s ≥ , Q S > ≤ λ < (cid:0) √ (cid:1) ≈ . • The line of fixed points P ( z c ) : (cid:0) , z c , z c (cid:1) exists for λ = 0. The eigenvalues are (0 , − z c , − z c ). Thus, it isnonhyperbolic. Using the Center Manifold Theory we obtain that it is asymptotically stable for 0 < z c ≤ − ≤ z c < φ, H are approximately constant, such that˙ φ ≈ , H ( t ) ≈ z c √ − z c , a ( t ) ≈ a e tzc √ − zc , z c = ± r V V , where the potential is constant V ( φ ) = V since λ = 0. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P ( z c ) are (cid:0) − , − , , − , − (cid:1) . However, c s < , Q S <
0. Thus, this solution suffers from Laplacian instabilities and thepresence of ghosts. • The line of fixed points P ( z c ) : (cid:0) βz c , √ αβz c , z c (cid:1) exists for λ = 0 , β = √ (cid:0) √ α − √ α − (cid:1) .The eigenvalues are (0 , − z c , − z c ). Thus, it is nonhyperbolic. The stability has to be ex-amined numerically or using Center Manifold calculations. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P ( z c ) are c s = ( α ( √ α − − α ( − α ( α ( α ( √ α − − α ) +7 ) − √ α − ))) +6 ) z c ( α ( √ α − − α ) +2 )( ( ( α ( √ α − − α ) +1 ) α +7 ) z c − ( √ α − − α ) √ − z c ) , Q S = ( √ α − − α ) √ − z c +3 ( α ( α −√ α − α − ) − ) z c ( α ( √ α − − α ) +1 ) z c , Ω DE = 1, ω DE = ( √ α − − α )( αz c −√ − z c ) z c , and q = − α/ (1 + α ) to cover all real values of α ).FIG. 5: Parameter space where P ( z c ) is free from Lapacian instabilities and it is ghost -free. • The line of fixed points P ( z c ) : (cid:0) βz c , √ αβz c , z c (cid:1) exists for λ = 0 , β = √ (cid:0) √ α + √ α − (cid:1) .The eigenvalues are (0 , − z c , − z c ). Thus, it is nonhyperbolic. The stability has to be ex-amined numerically or using Center Manifold calculations. The values of (cid:0) c s , Q S , Ω DE , ω DE , q (cid:1) for P ( z c ) are c s = ( − α ( α ( α ( α ( α ( √ α − α ) − ) − √ α − ) +51 ) +19 √ α − )) z c ( α ( √ α − α ) − )( ( α ( α ( √ α − α ) − ) − ) z c − ( √ α − α ) √ − z c ) , Q S = ( α ( α ( √ α − α ) − ) − ) z c − ( √ α − α ) √ − z c ( − α ( √ α − α )) z c , Ω DE = 1, ω DE = ( √ α − α )( √ − z c − αz c ) z c , and q = − P ( z c ) always satisfy c s Q S <
0. Thus, these solutions suffers either from Laplacianinstabilities or from the presence of ghosts.4For the lines of fixed points P ( z c ) and P ( z c ) the cosmological solutions satisfy˙ φ ( t ) ≈ √ βz c √ − z c , H ( t ) ≈ z c √ − z c , φ ( t ) − φ ≈ √ βz c t √ − z c , a ( t ) ≈ a e tzc √ − zc , z c = ± √ V p √ αβ + V , where β = √ α ∓√ α − √ respectively, and the potential is constant V ( φ ) = V since λ = 0.Finally, the fixed points/ fixed lines at infinity are • (cid:0) π , z c (cid:1) with eigenvalues (0 , • (0 , −
1) with eigenvalues (cid:0) − λ , − λ (cid:1) . This point attracts nearby orbits lying on the cylinder for all λ = 0.However, if we take into account the stability along r , the point is generically a saddle. • ( π, −
1) with eigenvalues (cid:0) − λ , − λ (cid:1) . This point attracts nearby orbits lying on the cylinder for all λ = 0.However, if we take into account the stability along r , the point is generically a saddle. • (0 ,
1) with eigenvalues (cid:0) λ , λ (cid:1) . This point repels nearby orbits lying on the cylinder for all λ = 0. However,if we take into account the stability along r , the point is generically a saddle. • ( π,
1) with eigenvalues (cid:0) λ , λ (cid:1) . This point repels nearby orbits lying on the cylinder for all λ = 0. However,if we take into account the stability along r , the point is generically a saddle.
4. ASYMPTOTIC EXPANSIONS IN THE REGIME WHERE THE CUBIC DERIVATIVE TERMDOMINATES
An important question about the present model is: in what phase of the cosmological evolution are the extra cubicinteraction terms expected to play an important role? Another interesting question is: how does the cosmologicalpredictions of our model differ from those of standard FRW cosmology driven by a scalar field with exponentialpotential? To answer these questions we proceed as follows.To show in what phase of the cosmological history it is expected that the cubic derivatives play a role, let us consider g ≫
1, and take the limit g to infinity in the equations for ˙ H, ¨ φ . We obtain the approximations¨ φ + 12 λ ˙ φ = 0 , ˙ H − λH ˙ φ + 3 H + 16 λ ˙ φ = 0 (31)that admits the solution H s ( t ) = 6 H − λ ˙ φ (cid:16) ( t − t ) (cid:16) H − λ ˙ φ (cid:17) + 2 (cid:17) + λ ˙ φ λ ˙ φ ( t − t ) + 6 , φ s ( t ) = φ + 2 λ ln λ ˙ φ t − t ) ! , (32)where we have chosen the initial conditions H ( t ) = H , φ ( t ) = φ , ˙ φ ( t ) = ˙ φ .Furthermore, we have a s ( t ) = a (cid:18) λ ˙ φ ( t − t ) (cid:19) (cid:18) t − t ) (cid:16) H − λ ˙ φ (cid:17)(cid:19) . (33)Taking initial conditions such that ˙ φ = 0 or 6 H − λ ˙ φ = 0, we obtain that a ( t ) ∝ t , that is, it corresponds toa radiation dominated universe. On the other hand, if we choose initial conditions such that 3 H − λ ˙ φ = 0, weobtain that a ( t ) ∝ t , that is, it corresponds to a universe dominated by dark matter. For other values of the initialconditions we can model the transition from a radiation dominated universe to a matter dominated one. However,from the leading terms in the Friedmann equations as g → ∞ we obtain − a ˙ φ e λφ (cid:16) λ ˙ φ − H (cid:17) = 0 . (34)Since we are looking for universes with a >
0, it immediately follows that Galileon modifications are particularlyrelevant for the radiation-dominated (early time universe) as shown in [16]. On the other hand, for g →
0, the modelis well-suited for describing the late-time universe, and we recover the standard quintessence results found elsewhere,5e.g., in the seminal work [49]. Now let us use the above solution as a seed solution to construct an asymptoticexpansion for large g when the potential is turned on.Let us define H = H s + εh + O ( ε ) φ = φ s + ε Φ + O ( ε ) , ε = g − . (35)and let us assume ˙ φ = 0 (without losing generality we can set φ = 0). Then we have at zeroth-order the equation4 ˙ φ (cid:16) λ ˙ φ − H (cid:17)(cid:16) λ ˙ φ ( t − t ) + 2 (cid:17) (cid:16) − H ( t − t ) + λ ˙ φ ( t − t ) − (cid:17) = 0 . (36)Since we have assumed ˙ φ = 0, the solution is H = λ ˙ φ . This implies that a ( t ) ∝ t . That is, it correspondsto a universe dominated by a dust fluid/dark matter. Substituting back the value for H in the Klein-Gordon, theRaychaudhuri and Friedmann equations, respectively, we have the following first-order equations:˙ φ (cid:16) − λ + 3 λ ˙ φ (cid:16) ( t − t )Φ ′′ ( t ) (cid:16) λ ˙ φ ( t − t ) + 4 (cid:17) + 2Φ ′ ( t ) (cid:16) λ ˙ φ ( t − t ) + 2 (cid:17)(cid:17) + 12Φ ′′ ( t ) + 6 (cid:17) − V = 0 , (37a)9 h ′ ( t ) (cid:16) λ ˙ φ ( t − t ) + 2 (cid:17) + λ (cid:16) − λ + 3 λ ˙ φ Φ ′ ( t ) (cid:16) λ ˙ φ ( t − t ) + 2 (cid:17) + 6 (cid:17) = 0 , (37b)2 (cid:16) ˙ φ (cid:16) − φ (cid:16) λ ˙ φ ( t − t ) + 2 (cid:17) ( λ Φ ′ ( t ) − h ( t )) + λ − (cid:17) − V (cid:17) = 0 . (37c)The solution of (37) that satisfies Φ( t ) = Φ , h ( t ) = h is given by:Φ = Φ + 4 (cid:16) V − h ˙ φ (cid:17) λ ˙ φ (cid:16) λ ˙ φ ( t − t ) + 1 (cid:17) − (cid:16) V − h ˙ φ (cid:17) λ ˙ φ + (cid:16)(cid:0) λ − (cid:1) ˙ φ + 12 V (cid:17) ln (cid:16) λ ˙ φ ( t − t ) + 1 (cid:17) λ ˙ φ , (38a) h = 3 h ˙ φ − V φ (cid:16) λ ˙ φ ( t − t ) + 1 (cid:17) + 2 V φ (cid:16) λ ˙ φ ( t − t ) + 1 (cid:17) , (38b)where we observe that lim t →∞ Φ( t ) = ∞ unless V ˙ φ = − λ , and lim t →∞ h ( t ) = 0 . Furthermore the relative errors aredefined to be E ( H s ) = | H − H s | H , E ( φ s ) = | φ − φ s | φ , E ( ˙ H s ) = | ˙ H − ˙ H s | ˙ H , E ( ˙ φ s ) = | ˙ φ − ˙ φ s | ˙ φ , (39)and satisfies the conditions lim t →∞ E (Φ s ( t )) = lim t →∞ E ( ˙Φ s ( t )) = ǫ (cid:16)(cid:0) λ − (cid:1) ˙ φ + 12 V (cid:17) λ ˙ φ , (40a)lim t →∞ E ( H s ( t )) = lim t →∞ E ( H s ( t )) = 2 V ǫλ ˙ φ . Thus, the relative errors can be small enough for V ≪ ˙ φ , λ = ±√ V ≪ ˙ φ , ˙ φ ≫ V and large values of ˙ φ ). We see that for V ≫ ˙ φ and finite ˙ φ , the relative errors are large,and the approximation fails. In the former case the potential makes the ideal gas solution unstable, as expected fora universe dominated by a dust fluid/dark matter. That is an important result since the cubic term in the Galileonfield provides a matter epoch.
5. THE GALILEON MODEL WITH MATTER
Until now we have considered the case with vacuum (i.e. ρ m = 0), and we have found at some fixed points c s < Q S < c s ≥ Q S >
0. In this section we study if the matter stabilizes the Galileon field, in the sense thatit helps to restore the above conditions. To begin with, the definition of c s changes to c S ≡ w H − w − w − ρ m w +9 w ; when ρ m > w tot = w DE .By using the variables x = ˙ φ p H + 1) , y = V e − λφ H + 1) , z = H √ H + 1 , (41)and the parameter α = g V , we obtain a dynamical system the same as before (16): x ′ = f ( x ) f ( x ) , y ′ = f ( x ) f ( x ) , z ′ = f ( x ) f ( x ) , (42) x = ( x, y, z ) and functions f − f are defined by (17). We have an auxiliary evolution equation for˜Ω m = ρ m H + 1) , ρ m = ρ m, a − . (43)The system can be written in the compact form x ′ = 2 qz (cid:0) αx z + √ y (cid:1) − √ αλx + 12 αx z − √ x y + √ y (cid:0) y − z (cid:1) αx , (44a) y ′ = y (cid:16) q + 1) z − √ λx (cid:17) , (44b) z ′ = (cid:0) z − (cid:1) , (44c)˜Ω ′ m = ˜Ω m z (cid:0) q + 1) z − (cid:1) (44d)where q , the deceleration parameter, can be expressed in terms of the phase-space variables ( x, y, z ). Observe thatthe evolution equation for matter decouples from the other evolution equations. The subtle difference, including dustmatter, is that the restriction (18) now becomes y + 2 x α ( −√ z + xλ ) + y ( x − z + ˜Ω m ) = 0 . (45)Since y ≥ , ˜Ω m ≥
0, the motion of the particle in the phase space is not just restricted to the given surface y + 2 x α ( −√ z + xλ ) + y ( x − z ) = 0, as it is the case in (16) and (17). Instead, the phase space is now the set n ( x, y, z ) : y + 2 x α ( −√ z + xλ ) + y ( x − z ) ≤ , y ≥ o . (46)There are recovered the fixed points for the vacuum case, that is, with ˜Ω m = 0; but some existence conditions changewhich is the main difference between the two cases. The existence conditions for P , P changes to α/λ ≥
0; theexistence conditions of P , P are changed by λ + 3 αλ − ≥
0. The stability conditions changes accordingly. Theeigenvalues and stability conditions for P , P remains unaffected when matter is included. The eigenvalues for P , changes to ± (cid:18) − , − √ λ (6 α ( αλ − − λ )+244 √ α +1 λ , + √ λ (6 α ( αλ − − λ )+244 √ α +1 λ (cid:19) , but the point will still be a saddle.The results shown in table II remains unaltered. We conclude then, that the matter background has no imprintin the Galileon field concerning the avoidance of Laplacian instabilities and concerning the absence of ghosts. Theconditions c s ≥ , Q S >
6. CONCLUSIONS
In this paper we have studied cubic Galileon cosmology with an extra conservation law from the dynamical systemsperspective. The novelty of this model is that was derived by the application of Noether’s theorem in the gravitationalLagrangian. Firstly, we have noticed that the fixed points A ± , B ± , C and D investigated in detail in [16], do not existin our scenario. So, our analysis has its own right, and it is complementary to all the analysis done before. This new7scenario admits power-law solutions. We have found a new asymptotic solution given by the fixed point solution P ,which satisfies H ≈ b (cid:18) − tt − b − b t (cid:19) , φ ≈ ln (cid:16) V αλ (2 a − b ) +2 λ (2 a − b ) +3 b (9 b − (cid:17) λ + 2 ln tλ , while the scale factor is calculated to be a ( t ) ≈ a e b t t b t ( b − b t ) b / = ( b − b t ) b / (cid:0) a + O (cid:0) t − ) (cid:1)(cid:1) .This point has not been obtained previously in [16] and in [17], since in these works the authors used H -normalization, which fails obviously when H = 0. This solution attracts the exact power-law solution previouslydescribed in [44] for the proper choice of free parameters. Moreover, we have obtained several solutions ( P - P , P , P , P ( x c ), P ( z c )), that violates the conditions c s ≥ , Q S > • The solutions P , ( x c ) can satisfy the conditions c s ≥ , Q S >
0, denoted by the yellow region with dot-dashedboundary in the Figure 4. Although, they can violate these conditions, for example, it is possible to have c s < , Q S ≤ • The fine-tuned solutions P and P , that exists in the unmodified (quintessence) scenario can be sink, respec-tively, source for 0 ≤ λ ≤ , and they satisfy c s ≥ , Q S > < λ < • The solutions P and P that exists for λ + 3 αλ − c s ≥ , Q S > ≤ λ < (cid:0) √ (cid:1) ≈ . • The line of fixed points P ( z c ) : (cid:0) βz c , √ αβz c , z c (cid:1) exists for λ = 0 , β = √ (cid:0) √ α − √ α − (cid:1) , and thecosmological solution is given by˙ φ ( t ) ≈ √ βz c √ − z c , H ( t ) ≈ z c √ − z c , φ ( t ) − φ ≈ √ βz c t √ − z c , a ( t ) ≈ a e tzc √ − zc , z c = ± √ V p √ αβ + V . This solution is free of Laplacian instabilities and ghosts-free in the region displayed in figure 5.Furthermore we have investigated the fixed points at infinity, and all of them are saddle points, so they have norelevance for the early or late time universe.The Galileon modifications are particularly relevant for the radiation-dominated (early time universe) as shown in[16] by investigating the regime where the coupling parameter satisfy g ≫
1. On the other hand, for g → g when the potential is turned on.All the previous results were found for vacuum Galileon ( ρ m = 0). We have seen that for some points c s <
0, andat some fixed points Q S < c s ≥ Q S >
0. Moreover, we have investigated if the matter stabilizesthe Galileon field, in the sense that it helps to restore the above conditions. When we include background matterin the form of dust, we recovered the fixed points for ˜Ω m = 0 but some existence conditions change, especially for P , P and P , P . The stability conditions changes accordingly. The eigenvalues and stability conditions for P , P remains unaffected when matter is included. The eigenvalues for P , changes, but the point will still bea saddle. The results shown in table II remains unaltered. We conclude then, that the matter background has noimprint in the Galileon field concerning the avoidance of Laplacian instabilities and concerning the absence of ghosts.The conditions c s ≥ , Q S > z > H >
0) for accommodating alate-time accelerated expansion phase. For simplicity, we will restrict the analysis to the invariant set of the flow ofthe system (16) given by z = 1. In this invariant set the total equation of state parameter is w tot = − y . Now, for thespecial parameters choice α = λ (3 − λ ) , the fixed points P ( x c ) merges to one fixed point, denoted by dS , in the plane( x, w tot ) with coordinates ( x, w tot ) = (cid:16)q λ, − (cid:17) that mimics a de Sitter solution. It is an attractor for 0 < λ < q .For 0 . < λ < . c s > , Q S ≥ w tot , x ), representing the typical behavior of our model for α = λ (3 − λ ) anddifferent choices of λ .We obtain also a transient phase given by the fixed point D := ( x, w tot ) = (cid:18) √ λ , (cid:19) , that mimics a dust-like solution.Furthermore, we have a scaling solution S := ( x, w tot ) = (cid:18) √ λ , − λ (cid:19) , which can be an attractor for λ < − √ or q < λ < √
3, or a saddle otherwise. The scaling solution is accelerating for − √ < λ < −√ < λ < √
3. Thephysical conditions c s > , Q S ≥ . < λ < . D := ( x, w tot ) = (cid:16) √ λ , (cid:17) which is unstable. In the Figure 6 it is portrayed in the plane w tot vs. x , the typical behavior of our model for α = λ (3 − λ ) and different choices of λ . On the left panel the attractor of the de Sitter point, dS , such that thevalue w tot = − D and the scaling solution S that crosses the phantom divide line (represent by a red dashed line),that eventually enter the w = − S is a stable spiral, and the de Sitter point is a transient one. There are trajectoriesconnecting the dust point D and the de Sitter solution dS , crossing the phantom divide. As in the previous case,some orbits remains on the phantom region. As we have briefly shown, this model has quite interesting cosmologicalfeatures, resembling the usual late-time dynamics usually found in conventional quintessence models for a non-zerovalue of α .Different subjects of study are still open for the integrable model in which new critical points exist. Hence, therequirement of the existence of a conservation law has a physical interpretation in the evolution of the model. Howeverthe exact nature of the physical interpretation for this model is still unknown, such an analysis is still in progress. Acknowledgments
This work was financial supported by FONDECYT grants, 1150246 (AG) and 3160121 (AP). AP thanks theDurban University of Technology and the University of the Witwatersrand for the hospitality provided while thiswork was performed. SJ would like to acknowledge the financial support from the National Research Foundation ofSouth Africa with grant number 99279. AP and SJ also acknowledge support from the DST-NRF Centre of Excellencein Mathematical and Statistical Sciences (CoE-MaSS). Opinions expressed and conclusions arrived at are those of theauthors and are not necessarily to be attributed to the CoE-MaSS. GL acknowledges DI-VRIEA for financial supportthrough Proyectos VRIEA Investigador Joven 2016 and Investigador Joven 2017.9
Appendix A: Nonhyperbolic fixed point at infinity
For analyzing the stability of Q we proceed as follows. Dividing by the positive quantity λ cos ( θ ), we have thatthe flow of the system (25) is equivalent to the flow of dθdT = ˆ h ( θ, z ) := z sin( θ ) cos( θ ) , (A1a) dzdT = ˆ h ( θ, z ) := (cid:0) z − (cid:1) , (A1b) drdT = ˆ h ( r, θ, z ) := r [ z (cos(2 θ ) − . (A1c)To improve the range and accuracy we calculate the diagonal first order Pade approximants[1 / ˆ h ( θ ) , [1 / ˆ h ( θ ) , [1 / ˆ h ( θ ) , around θ = π/
2, which yields the following approximate expressions dθdT ≈
12 ( π − θ ) z, dzdT ≈ − z , drdT ≈ − rz. (A2)Imposing the initial conditions z (0) = − δz, θ (0) = π δθ, r (0) = δr, δz ≪ , δθ ≪ , δr ≪ , we obtain the solution z = δze − T − , θ = δθe T + π , r = δre T . (A3)Observe that the small perturbations δθ, and δr grow exponentially, and z → − T → + ∞ . At the cylinderat infinity ( δr = 0) the orbits near Q are attracted by the invariant set z = −
1, but the angle departs fromthe equilibrium value θ = π as time goes forward. This means that Q is the local past attractor at the circle z = 1 , X + Y = 1 , Y ≥
0, as it is shown in figure 2. If we move to the interior of the phase space δr >
0, the orbitsare repelled by the cylinder at infinity at an exponential rate, although z → − T → + ∞ . Hence, Q is genericallya saddle. Appendix B: Fixed points lying on singular surfaces
The system (16) blows-up when f ( x, y, z ) = 0, and specially at the fixed points P , P , P . For this reason we hadexamined their stability numerically or by taking limits, since they lead to indeterminacy of the form 0 / P we introduce the new variables u = z − √ xλ , v = z , w = y , such that thesystem can be written in the canonical form u ′ v ′ w ′ = uvw + λu (cid:0) α (cid:0) λ − (cid:1) − λ (cid:1) + uv (cid:0) λ (cid:0) λ − α (cid:0) λ − (cid:1)(cid:1) − (cid:1) + v (cid:0) α (cid:0) λ − (cid:1) λ − λ + 3 (cid:1) (cid:0) − αλ (2 u − v ) − λ (2 u − v ) − v (cid:1) λ w (2 u − v ) + O ( | ( u, v, w ) T | ) . (B1)Given H , the vector space of the homogeneous polynomials of second order in u = ( u , u , u ), let us consider thelinear operator L ( ) J associated to J = , that assigns to h ( u ) ∈ H , the Lie bracket of the vector fields Ju h ( u ): L ( ) J : H → H h → L J h ( u ) = Dh ( u ) Ju − Jh ( u ) . (B2)where H the real vector space of vector fields whose components are homogeneous polynomials of degree 2. Thecanonical basis for the real vector space of 3-dimensional vector fields whose components are homogeneous polynomialsof degree 2 is given by H = span u , u u , u u , u , u u , u , u , u u , u u , u , u u , u u , u u , u u , u , u u , u . (B3)By computing the action of L ( ) J on each basis element on H we have L ( ) J (cid:0) H (cid:1) = span u u , u u , u , u , u u , u u , u , u u , u , u u , u u , u . (B4)Thus, the second order terms that are linear combinations of the six vectors in (B4) can be eliminated [48]. Todetermine the nature of the second order terms that cannot be eliminated we must compute the complementary spaceof (B4) which is G = span u , u u , u , u , u u , u . Following the above reasoning, we propose a transformation uvw = σ σ σ + h (1)2 ( σ , σ , σ ) h (2)2 ( σ , σ , σ ) h (3)2 ( σ , σ , σ ) (B5)where the h ( i )2 , i = 1 , , σ , σ , σ ). Up to this point h ( i )2 , i = 1 , , O (2) terms as much as possible.We choose h (1)2 ( σ , σ , σ ) h (2)2 ( σ , σ , σ ) h (3)2 ( σ , σ , σ ) = λ σ (2 αλ + 1) − λ σ σ (2 αλ + 1) + σ (cid:0) αλ + 2 λ + 3 (cid:1) . (B6)Hence, the normal form of the system (B1) is σ ′ = λσ (cid:0) α (cid:0) λ − (cid:1) − λ (cid:1) + σ σ (cid:18) λ (cid:0) λ − α (cid:0) λ − (cid:1)(cid:1) − (cid:19) + 98 σ (cid:0) α (cid:0) λ − (cid:1) λ − λ + 3 (cid:1) + O (3) , (B7a) σ ′ = σ + O (3) , (B7b) σ ′ = σ σ ( αλ + 1) − σ σ (cid:0) λ ( αλ + 1) + 3 (cid:1) + O (3) (B7c)1We propose the following anzats: σ = a τ + a τ + O ( τ − ) , σ = b τ + b τ + O ( τ − ) , σ = c τ + c τ + O ( τ − ) (B8)as τ → ∞ . Substituting back and comparing the coefficient of equal powers of τ we find the solutions: a = − p b (8 λ ( − αλ + 6 α + λ ) + 9 b (2 λ ( α (2 λ (3 α + λ ) − − λ ) + 3) −
12) + 4 − b (cid:0) αλ − αλ − λ + 3 (cid:1) + 24 λ ( α (2 λ − − λ ) ,c = 0 , c = − b , (B9a)or a = − − p b (8 λ ( − αλ + 6 α + λ ) + 9 b (2 λ ( α (2 λ (3 α + λ ) − − λ ) + 3) −
12) + 4 − b (cid:0) αλ − αλ − λ + 3 (cid:1) + 24 λ ( α (2 λ − − λ ) ,c = 0 , c = − b . (B9b)This approximated solution has the correct number of arbitrary constants: a , b , b . Then, we obtain x = (3 b − a ) λ √ τ + (3 b − a ) λ √ τ + O ( τ − ) , (B10a) y = 4(2 a − b ) αλ + 2(2 a − b ) λ + 3 b (9 b − τ + O ( τ − ) , (B10b) z = 3 b τ + 3 b τ + O ( τ − ) . (B10c)From here it is clear that the origin attracts the orbits as τ → ∞ .On the other hand, we get˙ φ ( t ( τ )) = (3 b − a ) λτ + (3 b − a ) λτ + O ( τ − ) , (B11a) H ( t ( τ )) = 3 b τ + 3 b τ + O ( τ − ) , (B11b) φ ( t ( τ )) = ln (cid:16) V a − b ) αλ +2(2 a − b ) λ +3 b (9 b − (cid:17) + 2 ln τλ + O ( τ − ) . (B11c)Furthermore dtdτ = 1 √ H + 1 = 1 − b τ + 3 (cid:0) b − b (cid:1) τ + O ( τ − ) (B12)which implies ∆ t = t ( τ ) − t = τ − b ln τ + 3 (cid:0) b − b (cid:1) τ + O ( τ − ) . (B13)Inverting the above expression we have τ = 3 (cid:0) b − b (cid:1) t + 9 b ln ∆ t t + 34 b ln ∆ t + ∆ t + O (∆ t − ) . (B14)Finally, ˙ φ ( t ) = − b λ (3 b − a ) ln t t + λ (3 b − a ) t + λ (3 b − a ) t + O ( t − ) , (B15a) H ( t ) = − b ln t t + 3 b t + 3 b t + O ( t − ) , (B15b) φ ( t ) = ln (cid:16) V αλ (2 a − b ) +2 λ (2 a − b ) +3 b (9 b − (cid:17) λ + 2 ln tλ + O ( t − ) , (B15c) a ( t ) = a t b (cid:0) b ln t + 3 b − b (cid:1) t + O (cid:0) t − (cid:1)! . (B15d)2For simplicity we set t = 0 , ∆ t = t .On the other hand, the new conservation law g a e λφ ˙ φ (cid:16) H − λ ˙ φ (cid:17) λ + 2 a (cid:16) λH − ˙ φ (cid:17) λ + I = 0 , has to be satisfied. Substituting the expression (B15), we obtain I a + 18 t b − (cid:0) a − b ) (cid:0) b − b (cid:1) + 3 b (9 b − a − b ) ln t + 32 a − b (cid:1) + (4 a − b ) t b − + O (cid:0) t − (cid:1) = 0 . where a = a ( b , λ, α ). Thus, for eliminating the term ∝ t b − , we have the following choices: • We set the exponent to zero, i.e., we set b = . Then a is reduced to the two values: a ± ( λ, α ) = − λ (cid:0) λ − α (cid:0) λ − (cid:1)(cid:1) ± √ p λ ( α (6 αλ − λ + 3) + λ )6 λ ( α (2 λ − − λ ) . Substituting b = , and specifying a = a ± ( λ, α ), the constraint becomes I a + a (cid:0) − b (cid:1) + 4 a − t + 4 a −
23 + O (cid:0) t − (cid:1) = 0 . Now, in order to satisfy the above expressions up to the order O (cid:0) t − (cid:1) , we impose the conditions a (6 − b ) + 36 a − , a −
23 + I a = 0 , that has four arbitrary constants a , b , a , I , from which we cannot eliminate simultaneously a and b , or I and a . This leads to a two-parameter family of solutions. We solve for a and I to keep the parameters b and a . Finally, ˙ φ ( t ) = (cid:0) − a (cid:1) λt + λ ( − a + 27 b − (1 − a ) ln t )9 t + O (cid:0) t − (cid:1) , (B16a) H ( t ) = 13 t + 27 b − ln t t + O (cid:0) t − (cid:1) , (B16b) φ ( t ) = ln (cid:16) V − a ) λ (2 αλ +1) − (cid:17) + 2 ln tλ + O (cid:0) t − (cid:1) , (B16c) a ( t ) = a √ t −
118 ( a (27 b − ln t − t − + O (cid:16) t − (cid:17) . (B16d) • We can equate the coefficient of t b − to zero, that is, to solve (4 a − b ) = 0 for b . The next termon the expansion will be (4 a − b ) t b − . Setting (4 a − b ) = 0, we obtain the constraints I = 0 and b = a , b = a , a (cid:0) a λ (cid:0) α (cid:0) λ − (cid:1) − λ (cid:1) + 1 (cid:1) = 0. • Finally, we have the fine tuned solution a = and b = , that leads to the constraints 2 αλ − αλ − λ = 0and I a + 4 a − b = 0.This point has not been obtained previously in [16] and in [17], since in these works the authors used H -normalization, which fails obviously when H = 0.For examining the stability of the fixed point P , we introduce the new variable ˜ z = z + 1. Then, the evolutionequations (16) can be written as x ′ = x (cid:18) − αλ (cid:19) + r λy + O (2) , y ′ = − y + O (2) , ˜ z ′ = − z + O (2) . (B17)For examining the stability of the the fixed point P , we introduce the new variable ˜ z = z −
1. Then, the evolutionequations (16) can be written as x ′ = − x (cid:18) − αλ (cid:19) + r λy + O (2) , y ′ = 3 y + O (2) , ˜ z ′ = 3˜ z + O (2) . (B18)From this linearized equations we extract the eigenvalues for P and P .3 Appendix C: Center Manifold for P The point P has the opposite stability behavior.For investigating the stability of the center manifold of P we introduce the new variables w = y, u = z + 1 , v = x − √ λ z. Hence, we have the evolution equations w ′ = 12 uw + λ (cid:18) √ vw − w α (cid:19) + λ w α + O (3) , (C1a) u ′ + 2 √ λuv + O (3) , (C1b) v ′ = − v + 9 uv + λ w ( u + w )4 √ α − q uw α + λ (cid:0) √ α v − αvw + √ w (cid:1) α + − λ w (cid:0) √ w − αv (cid:1) α + λ w √ α + O (3) . (C1c)The center manifold of P is then given locally by { ( w, u, v ) : u = h ( w ) , v = h ( w ) , h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 , | w | < δ } , where δ is smallenough, and the functions h , h satisfy the equations h ′ ( w ) − wh ( w ) − √ λwh ( w ) − λ (cid:0) λ − (cid:1) w α ! + h ( w ) (cid:16) √ λh ( w ) − (cid:17) + 15 h ( w ) = 0 , (C2a) h ′ ( w ) − wh ( w ) − √ λwh ( w ) − λ (cid:0) λ − (cid:1) w α ! + h ( w ) h ( w ) + (cid:0) λ − (cid:1) w √ α ! + h ( w ) λ (cid:0) λ − (cid:1) w α − ! + 3 r λh ( w ) + λw (cid:0) αλ + λ − λ + 12 (cid:1) √ α = 0 . (C2b)We have to solve this equations using Taylor series up to third order since the equations (C1) were truncatedat third order. Assuming h ( w ) = aw + O ( w ) and h ( w ) = bw + O ( w ) (we start at second order in w sincethe conditions h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 has to be satisfied), substituting back into (C2),and equating the coefficients with the same powers of w we find a = 0 , b = λ ( αλ + λ − λ +12 ) √ α . The evolution onthe center manifold is governed by w ′ = λ ( λ − ) w α + λ w ( αλ + λ − λ +12 ) α + O (cid:0) w (cid:1) , a “gradient” equation withpotential W ( w ) = − λ w ( αλ + λ − λ +12 ) α − λ ( λ − ) w α . Since the first non-zero derivative at w = 0 is the third one,then the point w = 0 is an inflection point of W ( u ). Hence, for λ ( λ − ) α >
0, the solutions with w (0) > ∞ in finite time) whereas the solutions with w (0) < λ ( λ − ) α <
0, the solutions with w (0) > w (0) < ∞ in finite time). Such an equilibrium with one-sided stability is sometimessaid to be semi-stable. However, since we have to restrict our attention to the region w ≥ y ≥ P is a saddle for λ ( λ − ) α >
0, and stable for λ ( λ − ) α <
0. But α is always non-negative, thus, P is stable for λ < −√ < λ < √
6, and saddle for −√ < λ < λ > √ Appendix D: Center Manifold for P We discuss in more detail the stability of the fixed point P . The point P has the opposite stability behavior.For investigating the stability of the center manifold of P we introduce the new variables w = z, u = y − z + 1 , v = x. Hence, we have the evolution equations w ′ = − uw + O (3) , (D3a) u ′ = − (cid:0) u + uw + u − v − w (cid:1) + O (3) , (D3b) v ′ = − v (cid:16) u + 2 (cid:16) √ αv + w + 1 (cid:17)(cid:17) + O (3) . (D3c)4The center manifold of P is then given locally by { ( w, u, v ) : u = h ( w ) , v = h ( w ) , h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 , | w | < δ } , where δ is smallenough, and the functions h , h satisfy the equations3 wh ( w ) h ′ ( w ) − h ( w ) − w + 1) h ( w ) + 3 h ( w ) + 3 w = 0 , (D4a)3 wh ( w ) h ′ ( w ) − h ( w ) h ( w ) − √ αh ( w ) − w + 1) h ( w ) = 0 . (D4b)We have to solve this equations using Taylor series up to third order since the equations (D3) were truncated atthird order. Assuming h ( w ) = aw + O ( w ) and h ( w ) = bw + O ( w ) (we start at second order in w since theconditions h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 has to be satisfied), substituting back into (D3), andequating the coefficients with the same powers of w we find a = 1 , b = 0. The evolution on the center manifold isgoverned by w ′ = − w . The equilibrium w = 0 is then asymptotically stable. Furthermore, perturbations from theequilibrium grow or decay algebraically in time, not exponentially as in the usual linear stability analysis. This is thesame result obtained previously in [16] and in [17], for the de-Sitter (constant potential) solution. Appendix E: Center manifold for the line of fixed points P ( z c ) The line of fixed points P ( z c ) : (cid:0) , z c , z c (cid:1) exists for λ = 0. The eigenvalues are (0 , − z c , − z c ). Thus, itis nonhyperbolic. For investigating the stability of the center manifold of P ( z c ) we introduce the new variables w = − ( y +1) z c + y +2 zz c − z c z c , u = ( z c − ) ( y + z c ( z c − z ))2 z c , v = x. Hence, we have the evolution equations w ′ = 3 u z c + 3 uw (cid:0) z c − (cid:1) z c z c − O (3) , (E5a) u ′ = − u (cid:0) z c + 4 z c − (cid:1) z c −
1) + u ( − w − z c ) + 32 v (cid:0) z c − (cid:1) + 32 w (cid:0) z c − (cid:1) + O (3) , (E5b) v ′ = uv (cid:0) − z c (cid:1) z c − − √ αv + v ( − w − z c ) + O (3) . (E5c)The center manifold of P is then given locally by { ( w, u, v ) : u = h ( w ) , v = h ( w ) , h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 , | w | < δ } , where δ is smallenough, and the functions h , h satisfy the equations − z c h ( w ) − w (cid:0) z c − (cid:1) z c h ( w ) z c − ! h ′ ( w ) − (cid:0) z c + 4 z c − (cid:1) h ( w ) z c − − h ( w )( w + z c )+ 32 (cid:0) z c − (cid:1) h ( w ) + 32 w (cid:0) z c − (cid:1) = 0 , (E6a) − z c h ( w ) − w (cid:0) z c − (cid:1) z c h ( w ) z c − ! h ′ ( w ) + (cid:0) − z c (cid:1) h ( w ) h ( w ) z c − − √ αh ( w ) − h ( w )( w + z c ) = 0 . (E6b)We have to solve this equations using Taylor series up to third order since the equations (E5) were truncated atthird order. Assuming h ( w ) = aw + O ( w ) and h ( w ) = bw + O ( w ) (we start at second order in w since theconditions h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 has to be satisfied), substituting back into (E5), andequating the coefficients with the same powers of w we find a = z c − z c , b = 0. The evolution on the center manifoldis governed by w ′ = − w z c (cid:0) − z c (cid:1) w . The equilibrium w = 0 is then asymptotically stable for 0 < z c ≤ − ≤ z c <
0. As before, perturbations from the equilibrium grow or decay algebraicallyin time.
Appendix F: Center manifold for the lines of fixed points P , ( z c ) The lines of fixed points P , ( z c ) : (cid:0) βz c , √ αβz c , z c (cid:1) exists for λ = 0 , β = β , , where β , = √ (cid:0) √ α ∓ √ α − (cid:1) , respectively. The eigenvalues are (0 , − z c , − z c ). Thus, they are nonhyperbolic. The twocases can be treated as one case under the choice α = β +1 √ β (this is equivalent to setting β = √ (cid:0) √ α ∓ √ α − (cid:1) ).5For investigating the stability of the center manifold of P , ( z c ) we introduce the new variables w = y − yz c β z c +2 z c − z c (cid:0) − zz c + z c + 1 (cid:1) , u = ( z c − )( y − ( β +1 ) z c (2 z − z c ) ) β +1) z c , v = x + β y ( zc − ) β + z c ( − zz c + z c − ) ! z c . Hence,we have the evolution equations w ′ = 3 u z c + 3 uw (cid:0) z c − (cid:1) z c z c − O (3) , (F7a) u ′ = 3 u (cid:0) β + z c + (cid:0) − β (cid:1) z c − (cid:1) β −
1) ( z c −
1) + u − βv (cid:0) z c − (cid:1) β − − w − z c ! + 3 v (cid:0) z c − (cid:1) β − w (cid:0) z c − (cid:1) + O (3) , (F7b) v ′ = 3 βu (cid:0) z c + 1 (cid:1) (cid:0) z c + 2 (cid:0) β − (cid:1) z c + 1 (cid:1) β −
1) ( z c − + 3 uv (cid:0) z c (cid:0) − β (cid:0) z c + 4 (cid:1)(cid:1) + 1 (cid:1) ( β −
1) ( z c −
1) + v (cid:0) β (cid:0) z c + 3 (cid:1) − (cid:1) β ( β − v ( − w − z c ) + 32 βw (cid:0) z c − (cid:1) + O (3) . (F7c)The center manifold of P is then given locally by { ( w, u, v ) : u = h ( w ) , v = h ( w ) , h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 , | w | < δ } , where δ is smallenough, and the functions h , h satisfy the equations − z c h ( w ) − w (cid:0) z c − (cid:1) z c h ( w ) z c − ! h ′ ( w ) + h ( w ) − β (cid:0) z c − (cid:1) h ( w ) β − − w + z c ) ! + 3 h ( w ) (cid:0) β + z c + (cid:0) − β (cid:1) z c − (cid:1) β −
1) ( z c −
1) + 3 (cid:0) z c − (cid:1) h ( w ) β −
1) + 32 w (cid:0) z c − (cid:1) = 0 , (F8a) − z c h ( w ) − w (cid:0) z c − (cid:1) z c h ( w ) z c − ! h ′ ( w ) + 3 h ( w ) h ( w ) (cid:0) z c (cid:0) − β (cid:0) z c + 4 (cid:1)(cid:1) + 1 (cid:1) ( β −
1) ( z c − β (cid:0) z c + 1 (cid:1) h ( w ) (cid:0) z c + 2 (cid:0) β − (cid:1) z c + 1 (cid:1) β −
1) ( z c − + 3 h ( w ) (cid:0) β (cid:0) z c + 3 (cid:1) − (cid:1) β ( β − − h ( w )( w + z c )+ 32 βw (cid:0) z c − (cid:1) = 0 . (F8b)We have to solve this equations using Taylor series up to third order since the equations (F7) were truncated at thirdorder. Assuming h ( w ) = aw + O ( w ) and h ( w ) = bw + O ( w ) (we start at second order in w since the conditions h (0) = 0 , h (0) = 0 , h ′ (0) = 0 , h ′ (0) = 0 has to be satisfied), substituting back into (F7), and equating thecoefficients with the same powers of w we find a = z c − z c , b = β ( z c − z c . The evolution on the center manifold isgoverned by w ′ = − w z c (cid:0) − z c (cid:1) w . The equilibrium w = 0 is then asymptotically stable for 0 < z c ≤ − ≤ z c <
0. As before, perturbations from the equilibrium grow or decay algebraicallyin time. [1] A. Nicolis, R. Rattazzi and E. Trincherini, Phys. Rev. D , 084003 (2009)[3] G.W. Horndeski, Int. J. Ther. Phys.
363 (1974)[4] C. Deffayet, S. Deser and G. Esposito-Farese, Phys. Rev. D , 055 (2014)[6] C. Charmousis, B. Gout´eraux and E. Kiritsis, JHEP , 11 (2012)[7] C. Charmousis and M. Tsoukals, Phys. Rev. D , 104050 (2015)[8] E. Babichev, C. Charmousis, A. Leh´ebel and T. Moskalets,JCAP , 011 (2016)[9] A. Cisterna, T. Delsate and M. Rinaldi, Phys. Rev. D 92, 044050 (2015)[10] A. Cisterna, T. Delsate, L. Ducobu and M. Rinaldi, Phys. Rev. D 93, 084026 (2016)[11] A. Barreira, B. Li, W.A. Hellwing, C.M. Baugh and S. Pascoli, JCAP , 027 (2013) [12] E. Bellini and R. Jimenez, Phys. Dark Matter , 179 (2013)[13] N. Bartolo, E. Bellini, D. Bertacca and S. Matarrese, JCAP , 034 (2014)[14] E. Babichev, C. Charmousis, A. Leh´ebel and T. Moskalets, JCAP , 011 (2016)[15] S. Bhattacharya, K.F. Dialektopoulos and T.N. Tomaras, JCAP , 036 (2016)[16] G. Leon and E.N. Saridakis, JCAP
025 (2013)[17] R. De Arcia, T. Gonzalez, G. Leon, U. Nucamendi and I. Quiros., Class. Quant. Grav. , no. 12, 125036 (2016)[18] M. Andrews, K. Hinterbichler, J. Stokes and M. Trodden, Class. Quant. Grav. , 184006 (2013)[19] C. Deffayet, S. Deser and G. Esposito-Farese, Phys. Rev. D , 044037 (2017)[21] M. Shahalam, S.K.J. Pacif and R. Myrzakulov, EPJC , 410 (2016)[22] N. Chow and J. Khoury, Phys. Rev. D , 024037 (2009)[23] A. De Felice and S. Tsujikawa, Phys. Rev. D , 124029 (2011)[24] C. Deffayet and D.A. Steer, Class. Quantum Grav. , 214006 (2013)[25] M. Tegmark et al., Astrophys. J.
702 (2004)[26] M. Kowalski et al., Astrophys. J.
749 (2008)[27] E. Komatsu et al., Astrophys. J. Suppl. Ser. , 330 (2009)[28] P. A. R. Ade et al., (Planck Collaboration), A&A.
A13 (2016)[29] P.A.R. Ade, et al. (Planck 2015 Collaboration), A&A
A20 (2016)[30] C. Burrage, C. de Rham, D. Seery and A.J. Tolley, JCAP , 014 (2011)[31] J. Ohashi and S. Tsujikawa, JCAP , 035 (2012)[32] F. Arroja, N. Bartolo, E. Dimastrogiovanni and M. Fasiello, JCAP 13, 005 (2013)[33] D. Regan, G. J. Anderson, M. Hull and D. Seery, JCAP , 015 (2015)[34] A. Barreira, B. Li, C. Baugh and S. Pascoli, JCAP , 059 (2014)[35] K.S. Jumar, J.C.B. Sanchez, C. Escamilla-Rivera, J. Marto and P.V. Moniz, JCAP , 063 (2016)[36] T. Kobayashi, M. Yamaguchi and J.’i. Yokoyama, Phys. Rev. Lett. , 124054 (2010)[39] A. De Felice and S. Tsujikawa, Phys. Rev. D , 025 (2012)[40] J. Neveu, V. Ruhlmann-Kleider, P. Astier, M. Besan¸con, J. Guy, A. M¨oller and E. Babichev, DOI: 10.1051/0004-6361/201628878[41] J. Neveu, V. Ruhlmann-Kleider, A. Conley, N. Palanque-Delabrouille, P. Astier, J. Guy, E. Babichev, A&A , A54(2013)[42] R. Gannouji and M. Sami, Phys. Rev. D , 1041 (2012)[44] N. Dimakis, A. Giacomini, S. Jamal, G. Leon and A. Paliathanasis, to appear in Phys. Rev. D, arXiv:1702.01603 [gr-qc].[45] A. De Felice and S. Tsujikawa, Phys. Rev. Lett. , 111301 (2010)[46] A. De Felice and S. Tsujikawa, JCAP , 007 (2012)[47] S. Appleby and E. V. Linder, JCAP (2012) 043[48] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos , Springer, New York (2003).[49] E. J. Copeland, A. R. Liddle and D. Wands, Phys. Rev. D57