Eur. Phys. J. C (2020) 80:137 https://doi.org/10.1140/epjc/s10052-020-7714-3 Regular Article - Theoretical Physics Global dynamics of Hořava–Lifshitz cosmology with non-zero curvature and a wide range of potentials Fabao Gao1,2,a , Jaume Llibre2,b 1 School of Mathematical Science, Yangzhou University, Yangzhou 225002, China 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Bellaterra, 08193 Barcelona, Catalonia, Spain Received: 6 November 2019 / Accepted: 4 February 2020 / Published online: 17 February 2020 © The Author(s) 2020 Abstract The global dynamics of a cosmological model Lifshitz cosmology with k = 0 and  = 0 was studied in based on Hořava–Lifshitz gravity in the presence of curvature [7]. is described by using the qualitative theory of differential In this paper we will investigate the global dynamics equations. of the Hořava–Lifshitz scalar field cosmology under the Friedmann–Lemaître–Robertson–Walker background space- time in the presence of curvature and no cosmological con- 1 Introduction stant term, i.e., k = 0 and  = 0. In recent years Hořava [1] proposed a spacetime asymmet- ric gravitational theory similar to Lifshitz’s scalar field the- ory, also known as Hořava–Lifshitz gravity. This theory has 2 The cosmological equations inspired a great deal of research for its applications in cos- mology and black hole physics (see [2–18] or the review To describe the cosmological equations in this section we articles [19,20] and the references therein). first give a brief review of the Hořava–Lifshitz gravitational With or without detailed-balance condition, few theory as was proposed in [1]. In one common version of researchers have published several papers on Hořava–Lifshitz this theory [5,21], its field content can be given by a spatial gravity through phase space analysis, see for instance [3– scalar N , and a spatial vector Ni . They are the ‘lapse’ and 7,17]. Among these papers there are also some of them ‘shift’ variables usually found in general relativity. Then the about the case of non-zero space curvature k. When the terms full metric can be written as related to k are dominant, Friedmann solution will be gener- ated. On the contrary, when the effective Hořava radiation ds 2 = −N 2 dt 2 + gi j (d x i + N i dt)(d x j + N j dt), dominates, the corresponding cosmology gives a solution Ni = gi j N j , (1) of “radiation-dominated like” [6]. According to whether the cosmological constant  is zero in flat (k = 0) and non-flat where gi j (i, j = 1, 2, 3) is the spatial metric. The scal- (k = 0) universe, Leon et al. [3–5,14] divided the Hořava– ing transformation of the coordinates is in the form of Lifshitz gravity into four cases: (1) k = 0,  = 0; (2) t → l 3 t, x i → lx i , under which both N and gi j are invari- k = 0,  = 0; (3) k = 0,  = 0; (4) k = 0,  = 0, ant, but N i is scaled to N i → l −2 Ni . and analyzed these phase spaces. More precisely, they either Based on the detailed-balance condition [1], the full grav- studied the global dynamics of the planar case of Hořava– itational action of Hořava–Lifshitz is represented by Lifshitz gravity with exponential potential by using the two-   2   dimensional Poincaré compactification, or discussed only √ κ2 local dynamics of this cosmological model with power- Sg = dtd 3 x g N K i j K ij − λK 2 − Ci j C i j κ2 2w 4 law potential but without investigate the dynamics close to infinity. The three-dimensional global dynamics of Hořava– μκ 2  i jm μ2 κ 2 + 2 √ Ril ∇ j Rkl − Ri j R i j 2w g 8   a e-mail: gaofabao@sina.com (corresponding author) μ2 κ 2 1 − 4λ 2 b e-mail: − R + R − 3 2 , (2) jllibre@mat.uab.cat 8(3λ − 1) 4 123 137 Page 2 of 13 Eur. Phys. J. C (2020) 80:137   √ j j with the Cotton tensor C i j =  i jm ∇k 4Ri − Rδi /(4 g) monomial potential V (φ) = (μφ)2n /2n with a positive con- and the extrinsic curvature K i j = (ġi j −∇i N j −∇ j Ni )/(2N ). stant μ and a natural number n. Thus the power-law √ poten- The standard general covariant antisymmetric tensor is rep- tial is f (s) = −s 2 /(2n), and we obtain ds/dt = 6xs 2 /n. √ Furthermore it can be followed from the first three transfor- resented by  i jm / g, and the indices are to raise and lower with the metric gi j . Besides λ is a dimensionless constant, mations in Eq. (6) and the first two equations in system (5) μ, κ and w are also constants. that x 2 + y 2 − z 2 = 1 and 3x 2 − 2z 2 = −H  /H . There- For the potential V (φ), we consider the gravitational fore the cosmological equations are reduced to the following action term as follows autonomous dynamical system    dx   √   3 √ 3λ − 1 φ̇ 2 = x 3x 2 − 2z 2 − 3 + 6s 1 − x 2 + z 2 , S = dtd x g N − V (φ) , (3) dt 4 N2   dz = z 3x 2 − 2z 2 − 2 , (7) and the metric N i = 0, gi j = a 2 (t)γi j , γi j d x i d x j = dt √ r 2 d 22 + dr 2 /(1 − kr 2 ), here a(t) is the dimensionless scale ds 6 2 factor for the expanding universe, γi j is a constant curvature = xs . dt n metric of maximally symmetric. The curvature k = 1, 0, −1 Our study will fully describes the global dynamics of sys- corresponding to the closed, flat, and open space, respec- tem (7) in the region G = (x, z, s) ∈ R3 : x 2 − z 2 ≤ 1 of tively. physical interest adding its boundary at infinity. It should be For simplicity we normalize κ 2 and N to the number one, noted from the derivation of system (7) that the space cur- then the cosmological equations can be described as vature k is implicit in the variables x and z. More details of 4V (φ) + (3λ − 1)φ̇ 2 system (7) can be found from Eqs. (113)–(115) of [4] or Eqs. H2 = 24(3λ − 1) (44)–(46) of [14]. μ  a 4 − 2μ2 ka 2 + μ2 k 2 2 2 − , 16a 4 (3λ − 1)2 3 Phase portraits on the invariant planes and surface 3 4V (φ) − (3λ − 1)φ̇ 2 Ḣ + H 2 = (4) 2 16(3λ − 1) In order to study the local phase portraits of the finite and 3μ  a − 2μ2 ka 2 + μ2 k 2 2 2 4 infinite equilibrium points, and the global phase portraits of + , 32a 4 (3λ − 1)2 system (7) in the region G, which is the meaningful region 2V  (φ) for cosmology, see again [4] or [14]. We start discussing the φ̈ + 3H φ̇ + = 0, phase portraits on its invariant planes and surface 3λ − 1 where H = ȧ/a is the Hubble parameter. z = 0, s = 0, x 2 − z 2 = 1. Note that the cosmological constant term is neglected in the present paper, i.e.,  = 0, so system (4) reduces to 3.1 The invariant plane z = 0 φ̇ 2 V (φ) μ2 k 2 H2 = + − , On this plane system (7) becomes 24 6(3λ − 1) 16a 4 (3λ − 1)2   dx √  φ̇ 2 μ2 k 2 = x 2 − 1 3x − 6s , Ḣ + = 4 , (5) dt 8 8a (3λ − 1)2 √ ds 6 2 2V  (φ) = xs . (8) φ̈ + 3H φ̇ + = 0. dt n 3λ − 1 The phase portraits of system (8) in the strip z = 0 and For simplicity the dimensionless transformation of the fol- x 2 − z 2 ≤ 1, i.e., in z = 0, −1 ≤ x ≤ 1, it has been lowing variables [4,5] is given as studied in [7], and the phase portraits is shown in Fig. 1, where √ φ̇ V (φ) μ the hyperbolic equilibrium point e0 = (0, 0) is a saddle, x= √ , y=√ √ , z= , both the semi-hyperbolic equilibrium points e1 = (1, 0) and 2 6H 6H 3λ − 1 4(3λ − 1)a 2 H e2 = (−1, 0) are saddle-nodes. V  (φ) V  (φ) V  (φ)2 s=− , f (s) ≡ − . (6) In order to describe the dark energy or dark matter in V (φ) V (φ) V (φ)2 the Hořava–Lifshitz universe more clearly, we introduce the Since the usual scalar field potential V (φ) admits multiple density parameters of dark energy and dark matter (see [3,5] representations, the potential function f (s) also has various for more details): D E = −k 2 z 2 −u 2 , M = x 2 +y 2 , as well mathematical expressions (see the papers [4,17,22,23] and as the corresponding equation-of-state parameters: ω D E = the references therein). In this paper we will investigate the (k 2 z 2 − 3u 2 )/(3k 2 z 2 + 3u 2 ), ω M = (x 2 − y 2 )/(x 2 + y 2 ). 123 Eur. Phys. J. C (2020) 80:137 Page 3 of 13 137 z x s Fig. 2 The phase portrait on the invariant plane s = 0 restricted to the region x 2 − z 2 ≤ 1 x On the local chart U1 (see Chapter 5 of [24] for more details on the Poincaré compactification) system (9) becomes Fig. 1 The phase portraits of the invariant plane z = 0 in the region du dv   −1 ≤ x ≤ 1 = uv 2 , = v −3 + 2u 2 + 3v 2 . (10) dt dt Since this system vanishes at v = 0, all the points at infin- Table 1 Equilibrium points and their characteristics on the invariant ity are equilibrium points. Taking the transformation with plane z = 0 for all s ∈ R respect to time dτ = vdt yields Equilibrium x s Stability ωM DE ωD E points du dv = uv, = −3 + 2u 2 + 3v 2 . (11) dτ dτ e0 0 0 Unstable −1 0 Arbitrary e1 1 0 Semi-hyperbolic 1 0 Arbitrary This√ system has two hyperbolic √ points at infinity, e3 = e2 −1 0 Semi-hyperbolic 1 0 Arbitrary (− 6/2, 0) and e4 = ( 6/2, 0), both of them√are unsta- ble hyperbolic saddle points with eigenvalues ± 6. On the local chart U2 system (9) writes du dv   Therefore the equilibrium points and their characteristics on = −uv 2 , = v 2 − 3u 2 + 2v 2 . (12) this invariant plane can be summarized in Table 1. dt dt Rescaling the time of system (12) by letting dτ = vdt we obtain 3.2 The invariant plane s = 0 du dv = −uv, = 2 − 3u 2 + 2v 2 . (13) On this plane system (7) becomes dτ dτ In view of (0, 0) is not an equilibrium point of system (13), dx   dz   we will not continue to study other infinite equilibrium points = x 3x 2 − 2z 2 − 3 , = z 3x 2 − 2z 2 − 2 , (9) of system (13) because they have been studied in local chart dt dt U1 . which has three equilibrium points e0 = (0, 0), e1 = (1, 0) Therefore the global phase portraits of system (9) can be and e2 = (−1, 0). Here e0 is a hyperbolic stable node with found in Fig. 2. In addition, similar to the previous section, eigenvalues −3 and −2, both e1 and e2 are unstable hyper- the equilibrium points on this invariant plane and their char- bolic nodes with eigenvalues 6 and 1. acteristics are shown in Table 2. 123 137 Page 4 of 13 Eur. Phys. J. C (2020) 80:137 Table 2 Equilibrium points and their characteristics on the invariant plane s = 0 Equilibrium x z Stability ωM DE ωD E points e0 0 0 Stable −1 0 Arbitrary e1 1 0 Unstable 1 0 Arbitrary e2 −1 0 Unstable 1 0 Arbitrary e3 ∞ ∞ Unstable – −∞ Arbitrary e4 ∞ ∞ Unstable – −∞ Arbitrary 3.3 The invariant surface x 2 − z 2 = 1 s First we prove that the surface x 2 − z 2 = 1 is invariant under x the flow of system (7). If l = l(x, z, s) = x 2 − z 2 − 1, then in order that the surface x 2 − z 2 = 1 be invariant we must Fig. 3 The phase portrait of the invariant surface x 2 − z 2 = 1 have ∂l ∂l ∂l Table 3 Equilibrium points and their characteristics on the invariant ẋ + ż + ṡ = K l, ∂x ∂z ∂s plane x 2 − z 2 = 1 for all s ∈ R for some polynomial K , and this is the case with K = Equilibrium x z Stability ωM DE ωD E √ points 2(3x 2 − 2z 2 − 6xs). On the surface x 2 − z 2 = 1 system (7) can be written as e1 1 0 Semi-hyperbolic 1 0 Arbitrary   ds √ e2 −1 0 Semi-hyperbolic 1 0 Arbitrary dx 6 2 = x x −1 , 2 = xs . (14) e5 ∞ ∞ Stable 1 −∞ Arbitrary dt dt n e6 ∞ ∞ Unstable 1 −∞ Arbitrary Then except for all the points on x = 0 which are equilibrium points, system (14) also admits two finite equilibrium points e1 = (1, 0) and e2 = (−1, 0). By using Theorem 2.19 of 3.2, the equilibrium points and their characteristics on this [24], we can find that both e1 and e2 are semi-hyperbolic invariant surface are collated in Table 3. saddle-nodes. On the local chart U1 system (14) becomes √ 3.4 The finite equilibrium points du 6 dv   =u u + v2 − 1 , = v v2 − 1 . (15) dt n dt It is noted that system (7) admits three finite equilibrium points p0 = (0, 0, 0) with eigenvalues {−3, −2, 0}, p1 = It√has two infinite equilibrium points e5 = (0, 0) and e6 = (1, 0, 0) and p2 = (−1, 0, 0) with the same eigenvalues ( 6n/6, 0), where e5 is a hyperbolic stable node with eigen- {6, 1, 0}. Here p1 and p2 are located at the intersection of values −1 of multiplicity two, and e6 is a hyperbolic unstable the invariant planes z = 0, s = 0 and the invariant surface saddle point with eigenvalues ±1. x 2 − z 2 = 1, corresponding to the equilibrium points e1 and On the local chart U2 system (14) writes e2 in Sects. 3.1–3.3, respectively. The origin p0 of system √ √ (7) lies at the intersection of the invariant planes z = 0 and du 6 dv 6 =u − u+u −v , 2 2 =− uv. (16) s = 0, which is the same point as the equilibrium point e0 dt n dt n studied in Sects. 3.1 and 3.2. Let dτ = udt we obtain √ √ du 6 2 dv 6 =− u+u −v , 2 =− v. (17) 4 Phase portraits on the Poincaré sphere at infinity dτ n dτ n The origin e7 = (0, 0) on the local √ chart U2 is a hyperbolic In order to describe the dynamics of system (7) at infinity. We stable node with eigenvalues − 6/n of multiplicity two. use the method of the three-dimensional Poincaré compactifi- In short the global phase portraits of system (14) is shown cation (see [25] for more details) x = 1/z 3 , z = z 1 /z 3 , s = in Fig. 3. Furthermore, similar to the previous Sects. 3.1 and z 2 /z 3 , then the analytical vector field of system (7) on the 123 Eur. Phys. J. C (2020) 80:137 Page 5 of 13 137 Table 4 Equilibrium points on the different local charts of the Poincaré sphere at the infinity of R3 Equilibrium points Eigenvalues u 11 = (0, 0, 0) {−3, −3, 0} √ 6 u 12 = − , 0, 0 {0, 0, 0} 2 √ 6 u 13 = , 0, 0 {0, 0, 0} 2 √ 6 u 14 = −1, n, 0 {−1, 1, −2n} 6 √ s 6 u 15 = 1, n, 0 {−1, 1, −2n} 6 √ z   6 n 3 3n u 16 = 0, ,0 − , ,3 2 1+n 1+n 1+n Fig. 4 The phase portrait of system (7) at infinity on local chart U1 u b0 = (b, 0, 0) 0, −3 + 2b2 , −3 + 2b2 u 31 = (0, 0, 0) {0, 0, 0} n)/n. The phase portrait on the Poincaré sphere at infinity on local chart U1 is shown in Fig. 4. local chart U1 becomes On the local chart U2 in view of Poincaré compactification dz 1 √   x = z 1 /z 3 , z = 1/z 3 , s = z 2 /z 3 , system (7) writes = z 1 z 32 − 6z 2 −1 + z 12 + z 32 , √   dt dz 1    = −z 1 z 32 + 6z 2 1 − z 12 + z 32 , dz 2 √ 1 dt = z 2 −3 + 2z 12 + 3z 32 + 6z 2 + 1 − z 12 − z 32 , √ dt n dz 2 6 dz 3  √   √   = z 2 2 − 3z 1 + 2 z 1 z 2 + 2z 32 , (21) = −z 3 z 12 −2+ 6z 2 + −3+ 6z 2 −1 + z 32 . dt n dt dz 3   (18) = z 3 2 − 3z 12 + 2z 32 . dt In the different local charts of the Poincaré sphere, the Since we want to study the infinity we take z 3 = 0, rescal- infinity of R3 corresponds to z 3 = 0, then system (18) has ing the time dτ = z 2 dt system (21) is equivalent to the equilibrium points listed in Table 4. Except that the equi- dz 1 √   librium point u 31 denotes the origin of the local chart U3 , the = 6 1 − z 12 , dτ rest equilibrium points lie in the local chart U1 . Moreover √ (22) the straight line s = 0 of the local chart U1 is filled with the dz 2 6 = 2 − 3z 1 +2 z1 z2 . equilibrium points u b0 for all b ∈ R. dτ n For the case z 3 = 0 system (18) is reduced to For any constant c, since (c, 0) is not the equilibrium point of dz 1 √   the system (22), we will not continue to discuss other infinite = − 6z 1 z 2 −1 + z 12 , equilibrium points of the system (22). Because this has been dt    dz 2 √ 1 discussed in the case U1 . = z 2 −3 + 2z 12 + 6z 2 + 1 − z 12 . (19) On the local chart U3 the three-dimensional Poincaré com- dt n pactification is x = z 1 /z 3 , z = z 2 /z 3 , s = 1/z 3 , then sys- After changing the of time dτ = z 2 dt system (19) becomes tem (7) becomes dz 1 √    √  = − 6z 1 −1 + z 12 , dz 1 6(1 + n) dτ   = z1 − z 1 + 3z 1 − 2z 2 − 3z 3 2 2 2 dz 2 √ 1 dt n = −3 + 2z 1 + 6z 2 + 1 − z1 . √   2 2 (20) dτ n + 6 z 22 + z 32 , Then this system has three equilibrium points  √  √ √ ei,1 , ei,2 dz 2 6   (23) and√ei,3 with coordinates (−1, 6n/6), (1, 6n/6) and = z2 − z 1 + 3z 1 − 2 z 2 + z 3 , 2 2 2 (0, 6n/(2(1 + n))), respectively. Here both dt n √ ei,1 and ei,2√are √ unstable saddle points with eigenvalues √ 6/n and√−2 6. dz 3 6 ei,3 is an unstable node with eigenvalues 6 and 6(1 + =− z1 z3. dt n 123 137 Page 6 of 13 Eur. Phys. J. C (2020) 80:137 For the case z 3 = 0 from system (23) we obtain  √  dz 1 6(1 + n) √ = z1 − z 1 + 3z 1 − 2z 2 + 6z 22 , 2 2 dt n √ (24) dz 2 6 = z2 − z 1 + 3z 12 − 2z 22 . dt n It is noted that (0, 0) is a linearly zero equilibrium point, the topological index is zero from the Poincaré–Hopf theory (see Theorem 6.30 of [24] for more details). In order to study its local phase portrait we shall use vertical blow-ups (see for (a) (b) instance [26]) by letting w = z 2 /z 1 , then we obtain       dz 1 √ 1 = z 12 6 − − 1 + w 2 + 3 − 2w 2 z 1 , dt n   (25) dw √ = − 6z 1 w −1 + w . 2 dt Eliminating the common factor z 1 of system (25) by changing the time dτ = z 1 dt it yields       dz 1 √ 1 = z1 6 − − 1 + w 2 + 3 − 2w 2 z 1 , (c) dτ n (26) dw √   Fig. 5 In a–c there are the local phase portraits of the equilibrium = − 6w −1 + w . 2 dτ points in systems (26), (25) and (24), respectively System (26) admits three equilibrium points ei,4 = (0, −1), ei,5 = (0, 1) and ei,6 = (0, 0) on z 1 = 0, where ei,4 and√ ei,5 are two√hyperbolic stable nodes with eigenvalues − 6/n and −2 6, √ ei,6 is a hyperbolic √ unstable saddle point with eigenvalues 6 and − 6(1+n)/n. The local phase portraits around ei,4 , ei,5 and ei,6 are shown in Fig. 5a. Note that there is a time rescaling dτ = z 1 dt between systems (25) and (26), so the direction of the trajectories in the local phase portraits of system (25) is opposite to that of Fig. 5a when z 1 < 0, see Fig. 5b for more details. In addition, all points on the w axis, i.e., z 1 = 0, are singularities of system (25). Thus the local phase portraits at (0, 0) of system (24) is shown in Fig. 5c. Then the phase portrait in the local chart U3 is shown z in Fig. 6. In summary joining the previous information we obtain x the global phase portraits at infinity in the Poincaré sphere in Fig. 7. Fig. 6 The phase portrait of system (24), i.e., the phase portrait at infinity in the local chart U3 of system (7) 5 Phase portraits inside the Poincaré ball restricted to x 2 − z2 ≤ 1 R3 : z ≥ 0, s ≥ 0. R4 : z ≥ 0, s ≤ 0. It is noted that system (7) is invariant under the symmetries Then due to the symmetries we only need to study the phase (x, z, s) → (−x, −z, −s) and (x, z, s) → (−x, z, −s), so portraits of system (7) in the region R1 . it is invariant under the symmetry with respect to the origin Combining the phase portraits of the invariant surface x 2 − and to the z-axis. Now we divide the Poincaré ball restricted z = 1 with the phase portraits of the planes z = 0, and 2 to x 2 − z 2 ≤ 1 into the following four regions: s = 0, together with the phase portrait at infinity, we get the phase portraits on the boundary of the region R1 as shown R1 : z ≤ 0, s ≥ 0. R2 : z ≤ 0, s ≤ 0. in Figs. 8, 9 and 10. Here we explain the definition of the 123 Eur. Phys. J. C (2020) 80:137 Page 7 of 13 137 s s z z x x Fig. 7 The global phase portrait at infinity in the Poincaré sphere. Fig. 9 Phase portrait in the back boundary of the region R1 The North Pole of the Poincaré ball is the equilibrium point u 31 . The symmetric points of u 12 and u 13 with respect to the center of the sphere are U12 and U13 , respectively s s z x z Fig. 10 Phase portrait in the bottom of the region R1 x Fig. 8 Phase portrait in the front boundary of the region R1 parabolic sector segment of U3 on the boundary F3 , but it is unstable on the back boundary surfaces B1 and B2 . three-dimensional cartesian coordinate system in this paper as follows: we regard the x z-plane as the horizontal plane in R3 , where the direction of the z-axis is horizontal towards the right. If the z-axis axis rotates 90◦ clockwise we get the x- 6 Dynamics in the interior of the region R1 axis. The s-axis is vertical upward, and then x zs constitutes a three-dimensional right-handed cartesian coordinate system. Note that the original system (7) admits the three finite equi- In order to show the phase portraits more clearly, the librium points p0 , p1 and p2 in the three-dimensional carte- boundary of the region R1 is divided into six surfaces accord- sian coordinate system. The dynamical behavior of the sys- ing to the orientation towards us, the back to us and the bottom tem inside the region R1 depends on the behavior of the flow area. See Fig. 11 for more details. It is noted that the North in the following surfaces and planes Pole u 31 of the Poincaré ball on the front boundary surfaces F1 and F2 is stable, and there is an elliptic sector and a stable h(x, z, s) = 0, g(x, z) = 0, x = 0, z = 0, s = 0, 123 137 Page 8 of 13 Eur. Phys. J. C (2020) 80:137 Fig. 11 The six boundary u 31 u 31 surfaces of the region R1 u 31 u 31 u 31 B1 B2 U 13 U 13 U 13 p2 F1 u 14 u 14 F3 p 2 p2 F2 p0 C p0 p1 p1 u 12 u 12 p1 u 12 s z s x z Fig. 13 The subregion R11 x Fig. 12 There are seven subregions inside the region R1 of the Poincaré ball where   √   h(x, z, s) = x 3x 2 − 2z 2 − 3 + 6s 1 − x 2 + z 2 , g(x, z) = 3x 2 − 2z 2 − 2. s These surfaces and planes divide the region R1 into seven different subregions R1i , i = (1, 2, . . . , 7), see Figs. 12, z 13, 14, 15, 16, 17, 18 and 19 for more details. It should be noted that R11 and R13 represent parallel tunnels inside R1 , x respectively. It is easy to verify that h > 0 in the subregions Fig. 14 The subregion R12 R12 , R14 , R15 and R16 , and h < 0 in the subregions R11 , R13 and R17 . Similarly, we can find that g > 0 in the subregions R11 , R12 , R16 and R17 , and g < 0 in the subregions R13 , R14 (opposite to the above front surface) contained in the surface and R15 . It should be noted that the dotted and solid lines in g = 0. Figs. 13, 14, 15, 16, 17, 18 and 19 are consistent with those In Table 5 we describe the behavior of ẋ, ż and ṡ in the in Fig. 12. seven subregions R11 , . . . , R17 . From this table we obtain As shown in the subregion R11 (see Fig. 13) the front that the variables x and z in the subregion R11 decrease mono- surface consists of three dashed lines and one solid line con- tonically, and the variable s increases monotonically, so an tained in the surface h = 0, and the surface on the back side orbit in the subregion R11 either crosses the boundary h = 0 123 Eur. Phys. J. C (2020) 80:137 Page 9 of 13 137 s z x s Fig. 15 The subregion R13 z x Fig. 18 The subregion R16 s z x s Fig. 16 The subregion R14 z x Fig. 19 The subregion R17 these orbits in backwards time come from the equilibrium point p1 in the subregion R11 . This process can be simply summarized as follows s R12 z x p1 R11 R13 . Fig. 17 The subregion R15 R14 and enters into the subregion R12 , or crosses the boundary g = 0 and goes to the subregion R13 , or crosses the inter- The front boundary of R12 (see Fig. 14), i.e. the surface section curve of the subregions R11 , R12 and R13 into the containing the equilibrium points u 12 , u 14 , u 31 and p1 , is subregion R14 with very low probability, i.e. an orbit in the contained in the invariant surface x 2 − z 2 = 1 (x > 0). The subregion R11 will not stay in the future in this region, but surface on the back side of R12 that contains the equilibrium will exit through its boundaries into other subregions, and points p1 and u 31 consists of two parts: the upper part surface 123 137 Page 10 of 13 Eur. Phys. J. C (2020) 80:137 Table 5 Dynamical behavior in Subregions Corresponding region Increase or decrease seven different subregions R11 h < 0, g > 0, x > 0, z < 0, s > 0 ẋ < 0, ż < 0, ṡ > 0 R12 h > 0, g > 0, x > 0, z < 0, s > 0 ẋ > 0, ż < 0, ṡ > 0 R13 h < 0, g < 0, x > 0, z < 0, s > 0 ẋ < 0, ż > 0, ṡ > 0 R14 h > 0, g < 0, x > 0, z < 0, s > 0 ẋ > 0, ż > 0, ṡ > 0 R15 h > 0, g < 0, x > 0, z < 0, s > 0 ẋ > 0, ż > 0, ṡ < 0 R16 h > 0, g > 0, x > 0, z < 0, s > 0 ẋ > 0, ż < 0, ṡ < 0 R17 h < 0, g > 0, x > 0, z < 0, s > 0 ẋ < 0, ż < 0, ṡ < 0 containing the infinite equilibrium point u 31 is the intersec- s = 0 in R13 u 31 tion of the subregions R12 and R14 , which is contained in the surface g = 0, and the lower part surface containing the finite equilibrium point p1 , is the intersection of subre- R13 gions R11 and R12 that is contained in the surface h = 0. The equilibrium points in the subregion R12 are the finite equilibrium point p1 , the infinite equilibrium points on the Poincaré sphere at s = 0 and the equilibrium point u 14 at R11 R14 . infinity. However an orbit in the subregion R12 does not start from these infinite equilibrium points at s = 0 or return to these points, but it comes from the finite equilibrium point In subregion R14 (see Fig. 16) the front triangular surface p1 or from the subregions R11 and R14 which have common is contained in g = 0, and its opposite triangular surface boundaries with the subregion R12 . This is due to the fact that is on the back side, which is included in the plane x = 0. the variable z is decreasing monotonically and the variables According to Table 5 an orbit of the subregion R14 may come x and s are increasing monotonically inside the subregion from the infinite equilibrium points on x = 0, or come from R12 according to Table 5. Moreover a trajectory in the subre- the subregion R13 through the surface h = 0, or from subre- gion R12 will tend to the equilibrium point u 14 at infinity. In gion R15 traversing plane x = 0, then go through the surface summary, this dynamic behavior process can be represented g = 0 into the subregion R12 , or directly goes to the infi- as nite equilibrium point u 31 at the North Pole of the Poincaré sphere. This can be represented as p1 R13 R12 R11 R12 u 14 . R15 R14 R14 x = 0 in R14 u 31 . In the subregion R13 (see Fig. 15) the front surface com- posed of three dashed lines and one solid line is contained in the surface g = 0, and its opposite surface is contained It can be noted from Fig. 17 that the surface on the left of in the surface h = 0. From Table 5 an orbit in the subregion subregion R15 (consisting of two dashed lines and one solid R13 may only comes from the equilibrium points locate at line) is contained in plane x = 0, and the opposite surface to the infinity on s = 0 in R13 , or it comes from the subregions it is included in the surface g = 0. According to Table 5 an R11 . Then the orbit passes through the intersection of sub- orbit in the subregion R15 can only come from the adjacent regions R13 and R14 (contained in the surface h = 0) and subregion R16 , and then either enter into the subregion R14 enters the subregion R14 , or tends to the infinite equilibrium adjacent to it, or directly tend to the finite equilibrium point point u 31 , i.e., lies in the North Pole of the Poincaré sphere. p0 located at the center of the Poincaré ball. This dynamics So we obtain the following dynamics can be denoted as 123 Eur. Phys. J. C (2020) 80:137 Page 11 of 13 137 R14 have α-limit and ω-limit at the infinite equilibrium point u 31 , i.e. the North Pole of the Poincaré sphere. R16 R15 For the sake of summary, similar to the previous Sects. 3.1– 3.3, the equilibrium points on the region R1 and their char- p0 . acteristics within this region are integrated in Table 6. In this way the qualitative global dynamic behavior of system (7) is described. We observe that the curved surfaces of the left and right sides of the subregion R16 (see Fig. 18) are contained in the surfaces g = 0 and h = 0, respectively. It is known from 7 Discussion and conclusions Table 5 that an orbit in the subregion R16 can only come from the subregion R17 and crossing the surface g = 0, it From a cosmological point of view, the unstable finite equi- enters into the subregion R15 . Thus we get that librium points p1 and p2 as shown in Figs. 8, 9 and 10 or R17 R16 R15 . 12 are consistent with the results of case 1 in the cosmolog- ical model in [3,5], corresponding to complete dark matter The left and right surfaces of subregion R17 (see Fig. 19) domination. It is known from the last flow chart in Sect. 6 are included in the surface h = 0 and in the invariant surface that the physical importance of the finite equilibrium point x 2 − z 2 = 1 (x < 0), respectively. According to Table 5, it is p0 is a late-time state of the universe if the initial condition known that all three variables in subregion R17 are decreas- is not on the invariant plane z = 0 or on the invariant sur- ing monotonically, so an orbit in this subregion must come face x 2 − z 2 = 1. We also note that√this late-time√solution from the infinite equilibrium point u 31 located at the North was already obtained in the strip − 6/2 < s < 6/2 for Pole of Poincaré sphere and eventually tends to the infinite the cosmological model with exponential potential in [3–5], equilibrium point U13 in the invariant plane s = 0. Then we but our results hold for all values of s in the region G of the have that three-dimensional space. In addition, the main differences between the Ref. [5] and the present paper are: u 31 R17 U13 . 1. The former specifically studied the case of exponential The dynamic behavior of the orbits inside the seven sub- potential, our research is based on a wide range of poten- regions of R1 discussed above can be summarized as tials, and specifically considering the case of power-law p0 x = 0 in R14 potential. 2. The former is based on the planar cosmological model, and we are considering a three-dimensional one. 3. The former did not discussed the dynamics at infinity, but R16 R15 R14 R12 u 14 we studied the infinite equilibrium points, and found that the infinite equilibrium points u 14 and U13 may also be the late-time state solutions of the universe. U13 R17 u 31 R13 R11 p1 . See the last flow chart in Sect. 6 again, the equilibrium points u 14 and U13 will be future attractors located at infinity. s = 0 in R13 When the initial conditions are inside the region R1 , they can be the late-time states of the universe. It is noted that only The above flow chart shows that the orbits of system (7) u 14 and U13 are studied in the region R1 because system (7) contained in the interior of the region R1 have α-limit at admits two aforementioned symmetries. Thus the results of the finite equilibrium point p1 , and the some orbits on the u 15 and u 13 (see Fig. 7) are identical with those of u 14 and boundary of the region R1 have α-limit at the x = 0 in R14 , U13 , respectively. It should also be noted that Fig. 5.5 in [3] or s = 0 in R13 . Moreover the orbits have ω-limit either at or Fig. 3 in [5] can be obtained in Fig. 7 in the present paper the finite equilibrium point p0 , or at the infinite equilibrium when s = 0.6. Then it is easy to find that u 14 corresponds points u 14 and U13 , where u 14 is located on the intersec- to P6 in the Fig. 5.5 in [3] for the plane case when s → 0. tion curve of the Poincaré sphere and the invariant surface However the infinite equilibrium point u 14 can be the late- x 2 − z 2 = 1 (x > 0) at infinity, and U13 is located on the time state of the universe from the above analysis, while P6 intersection of the Poincaré sphere and the invariant surfaces corresponds to the universe dominated by dark matter due to s = 0, x 2 − z 2 = 1 (x < 0) at infinity (see Fig. 12). Fur- the absence of  and the effect of k will decrease with the thermore the orbits of system (7) on the boundary of R1 also increasing of the scale factor a(t). 123 137 Page 12 of 13 Eur. Phys. J. C (2020) 80:137 Table 6 Equilibrium points and Equilibrium points x z s Stability ωM DE ωD E their characteristics in the interior of the region R1 for all p0 0 0 0 Stable −1 0 Arbitrary s∈R p1 1 0 0 Unstable 1 0 Arbitrary p2 −1 0 0 Unstable 1 0 Arbitrary u 12 ∞ ∞ 0 Unstable – −∞ Arbitrary U13 ∞ ∞ 0 Stable – −∞ Arbitrary u 14 ∞ ∞ ∞ Stable – −∞ Arbitrary u b0 ∞ 0 0 Unstable – 0 Arbitrary u 31 0 0 ∞ Unstable −1 0 Arbitrary As a conclusion, based on Hořava–Lifshitz gravity in References a Friedmann–Lemaître–Robertson–Walker space–time with non-zero curvature and without the cosmological constant 1. P. Hořava, Quantum gravity at a Lifshitz point. Phys. Rev. D 79, 084008 (2009) term, and combined with the Refs. [7,27], it is found that if 2. E.M.C. Abreu, A.C.R. Mendes, G. Oliveira-Neto et al., Hořava– the universe is non-flat (i.e., k = 0), it may also be closed, Lifshitz cosmological models with noncommutative phase space then the finite equilibrium point p0 , and the infinite equilib- variables. Gen. Relat. Gravit. 51, 95 (2019) rium points u 14 and U13 can be the late-time states of the 3. G. Leon, C.R. Fadragas, Cosmological dynamical systems: and their applications (Lambert Academic Publishing, GmbH & Co., universe. Besides Eq. (6) also imply that the value of Hub- KG, Saarbrucken, 2012) ble parameter H tends to be zero (i.e., the final state is the 4. G. Leon, A. Paliathanasis, Extended phase-space analysis of the Minkowski point) in the cosmological model as time goes Hořava–Lifshitz cosmology. Eur. Phys. J. C 79, 746 (2019) on. 5. G. Leon, E.N. Saridakis, Phase-space analysis of Hořava–Lifshitz cosmology. J. Cosmol. Astropart. Phys. 2009, 006 (2009) 6. S. Carloni, E. Elizalde, P.J. Silva, An analysis of the phase space Acknowledgements The authors thank the anonymous reviewer whose of Hořava–Lifshitz cosmologies, in Cosmology, Quantum Vacuum comments and suggestions helped improve and clarify this manuscript. and Zeta Functions, Springer Proceedings in Physics, vol. 137, ed. The first author gratefully acknowledges the support of the National by S.D. Odintsov, D. Sáez-Gómez, S. Xambó-Descamps (Springer- Natural Science Foundation of China (NSFC) through Grant Nos. Verlag, Berlin, 2011), pp. 139–148 11672259 and 11571301, the China Scholarship Council through Grant 7. F.B. Gao, J. Llibre, Global dynamics of the Hořava–Lifshitz cos- No. 201908320086, the Ministry of Land and Resources Research of mological system. Gen. Relat. Gravit. 51(11), 152 (2019) China in the Public Interest through Grant No. 201411007. 8. B. Chen, On Hořava–Lifshitz cosmology. Chin. Phys. C 35(5), The second author gratefully acknowledges the support of the Min- 429–435 (2011) isterio de Economı́a, Industria y Competitividad, Agencia Estatal de 9. X. Gao, Y. Wang, R. Brandenberger, A. Riotto, Cosmological per- Investigación grants MTM2016-77278-P (FEDER) and MDM-2014- turbations in Hořava–Lifshitz gravity. Phys. Rev. D 81, 083508 0445, the Agència de Gestió d’Ajuts Universitaris i de Recerca grant (2010) 2017SGR1617, and the H2020 European Research Council Grant 10. S. Lepe, J. Saavedra, On Hořava–Lifshitz cosmology. Astrophys. MSCA-RISE-2017-777911. Space Sci. 350, 839–843 (2014) Data Availability Statement This manuscript has no associated data 11. M. Li, Y. Pang, A trouble with Hořava–Lifshitz gravity. J. High or the data will not be deposited. [Authors’ comment: This manuscript Energy Phys. 08, 015 (2009) mainly involves qualitative analysis, not quantitative and numerical 12. O. Luongo, M. Muccino, H. Quevedo, Kinematic and statistical analysis, so there is no associated data.] inconsistencies of Hořava–Lifshitz cosmology. Phys. Dark Univ. 25, 100313 (2019) Open Access This article is licensed under a Creative Commons Attri- 13. N.A. Nilsson, E. Czuchry, Hořava–Lifshitz cosmology in light of bution 4.0 International License, which permits use, sharing, adaptation, new data. Phys. Dark Univ. 23, 100253 (2019) distribution and reproduction in any medium or format, as long as you 14. A. Paliathanasis, G. Leon, Cosmological solutions in Hořava– give appropriate credit to the original author(s) and the source, pro- ifshitz gravity, arXiv preprint arXiv:1903.10821 (2019) vide a link to the Creative Commons licence, and indicate if changes 15. E.N. Saridakis, Aspects of Hořava–Lifshitz cosmology. Int. J. Mod- were made. The images or other third party material in this article ern Phys. D 20(08), 1485–1504 (2011) are included in the article’s Creative Commons licence, unless indi- 16. A. Tawfik, E. Abou El Dahab, FLRW cosmology with Hořava– cated otherwise in a credit line to the material. If material is not Lifshitz gravity: impacts of equations of state. Int. J. Theoret. Phys. included in the article’s Creative Commons licence and your intended 56(7), 2122–2139 (2017) use is not permitted by statutory regulation or exceeds the permit- 17. C.R. Fadragas, G. Leon, E.N. Saridakis, Dynamical analysis of ted use, you will need to obtain permission directly from the copy- anisotropic scalar-field cosmologies for a wide range of potentials. right holder. To view a copy of this licence, visit http://creativecomm Classical Quant. Gravity 31, 075018 (2014) ons.org/licenses/by/4.0/. 18. M. Bhattacharjee, Gravitational radiation and black hole forma- Funded by SCOAP3 . tion from gravitational collapse in theories of gravity with broken Lorentz symmetry, Baylor University. ProQuest Dissertations Pub- lishing 22585106 (2019) 123 Eur. Phys. J. C (2020) 80:137 Page 13 of 13 137 19. S. Mukohyama, Hořava–Lifshitz cosmology: a review. Classical 24. F. Dumortier, J. Llibre, J.C. Artés, Qualitative Theory of Planar Quant. Gravity 27, 223101 (2010) Differential Systems (Springer-Verlag, Berlin, 2006) 20. T.P. Sotiriou, Hořava–Lifshitz gravity: a status report. J. Phys. Conf. 25. A. Cima, J. Llibre, Bounded polynomial vector fields. Trans. Am. Series 283, 012034 (2011) Math. Soc. 318(2), 557–579 (1990) 21. E. Kiritsis, G. Kofinas, Hořava–Lifshitz cosmology. Nuclear Phys. 26. M.J. Álvarez, A. Ferragut, X. Jarque, A survey on the blow up B 821, 467–480 (2009) technique. Int. J. Bifurc. Chaos 21(11), 3103–3118 (2011) 22. D. Escobar, C.R. Fadragas, G. Leon, Y. Leyva, Asymptotic behav- 27. E. Di Valentino, A. Melchiorri, J. Silk, Planck evidence for a closed ior of a scalar field with an arbitrary potential trapped on a Randall– universe and a possible crisis for cosmology. Nat. Astron. (2019). Sundrum’s braneworld: the effect of a negative dark radiation term https://doi.org/10.1038/s41550-019-0906-9 on a Bianchi I brane. Astrophys. Space Sci. 349, 575–602 (2014) 23. A. Alho, J. Hell, C. Uggla, Global dynamics and asymptotics for monomial scalar field potentials and perfect fluids. Classical Quant. Gravity 32, 145005 (2015) 123