This is a preprint of: “Global dynamics of the Horava-Lifshitz cosmology with non-zero curvature and a wide range of potentials”, Fabao Gao, Jaume Llibre, Eur. Phys. J. C Part. Fields, vol. 80, 137:1–13, 2020. DOI: [10.1140/epjc/s10052-020-7714-3] GLOBAL DYNAMICS OF THE HOŘAVA-LIFSHITZ COSMOLOGY WITH NON-ZERO CURVATURE FABAO GAO1,2 , AND JAUME LLIBRE2 Abstract. The global dynamics of a cosmological model based on Hořava-Lifshitz gravity in the presence of curvature is described by using the qualitative theory of differential equations. 1. Introduction In recent years Hořava [1] proposed a spacetime asymmetric gravi- tational theory similar to Lifshitz’s scalar field theory, also known as Hořava-Lifshitz gravity. This theory has inspired a great deal of re- search for its applications in cosmology and black hole physics (see [2]-[17] or the review articles [18], [19] and the references therein). Here we will investigate the global dynamics of the Hořava-Lifshitz scalar field cosmology under the Friedmann-Lemaı̂tre-Robertson-Walker background spacetime in the presence of curvature and no cosmological constant term. The corresponding dimensionless field equations admit the following form dx √ = x (3x2 − 2z 2 − 3) + 6s (1 − x2 + z 2 ) , dt dz (1) = z (3x2 − 2z 2 − 2) , dt ds √ = −2 6xf (s), dt Key words and phrases. Hořava-Lifshitz; Global dynamics; Cosmology; Poincaré compactification. 1 2 FABAO GAO1,2 , AND JAUME LLIBRE2 where the power law potential f (s) = −s2 /(2n) with a natural number n. See equations equations (113)-(115) of [8] or (44)-(46) of [14] for more details. The Hořava-Lifshitz cosmological model with zero curvature and without cosmological constant term was studied in [5]. Furthermore, in the references [3, 7, 8, 9, 14], the authors either studied the global dynamics of the planar case of system (1) considering only the vari- ables x and z by using the two-dimensional Poincaré compactification, or discussed only local dynamics of system (1) without investigate the dynamics close to infinity. Our study fully describes the global dy- namics of system (1) in the region G = {(x, z, s) ∈ R3 : x2 − z 2 ≤ 1} of physical interest. 2. Phase portraits on the invariant planes and surface In order to study the local phase portraits of the finite and infinite equilibrium points, and the global phase portraits of system (1) in the region G, which is the meaningful region for cosmology, see again [8] or [14]. We start discussing the phase portraits on its invariant planes and surface z = 0, s = 0, x2 − z 2 = 1. 2.1. The invariant plane z = 0. On this plane system (1) becomes dx √  = (x2 − 1) 3x − 6s , dt √ (2) ds 6 2 = xs . dt n The phase portraits of system (2) in the strip z = 0 and x2 − z 2 ≤ 1, i.e. in z = 0, −1 ≤ x ≤ 1, it has been studied in [5], and the phase 3 portraits is shown in Figure 1, where the hyperbolic equilibrium point e0 = (0, 0) is a saddle, both the semi-hyperbolic equilibrium points e1 = (1, 0) and e2 = (−1, 0) are saddle-nodes. s x Figure 1. The phase portraits of the invariant plane z = 0 in the region −1 ≤ x ≤ 1. 2.2. The invariant plane s = 0. On this plane system (1) becomes dx dz (3) = x (3x2 − 2z 2 − 3) , = z (3x2 − 2z 2 − 2) , dt dt which has three equilibrium points e0 = (0, 0), e1 = (1, 0) and e2 = (−1, 0). Here e0 is a hyperbolic stable node with eigenvalues -3 and -2, both e1 and e2 are unstable hyperbolic nodes with eigenvalues 6 and 1. 4 FABAO GAO1,2 , AND JAUME LLIBRE2 On the local chart U1 (see Chapter 5 of [20] for more details on the Poincaré compactification) system (3) becomes du dv (4) = uv 2 , = v (−3 + 2u2 + 3v 2 ) . dt dt Since this system vanishes at v = 0, all the points at infinity are equilib- rium points. Taking the transformation with respect to time dτ = vdt yields du dv (5) = uv, = −3 + 2u2 + 3v 2 . dτ dτ √ This system has two hyperbolic points at infinity, e3 = (− 6/2, 0) and √ e4 = ( 6/2, 0), both of them are unstable hyperbolic saddle points √ with eigenvalues ± 6. On the local chart U2 system (3) writes du dv (6) = −uv 2 , = v (2 − 3u2 + 2v 2 ) . dt dt Rescaling the time of system (6) by letting dτ = vdt we obtain du dv (7) = −uv, = 2 − 3u2 + 2v 2 . dτ dτ In view of (0, 0) is not an equilibrium point of system (7), we will not continue to study other infinite equilibrium points of system (7) because they have been studied in local chart U1 . Therefore the global phase portraits of system (3) can be found in Figure 2. 2.3. The invariant surface x2 − z 2 = 1. First we prove that the surface x2 − z 2 = 1 is invariant under the flow of system (1). If l = l(x, z, s) = x2 − z 2 − 1, then in order that the surface x2 − z 2 = 1 be 5 z x Figure 2. The phase portrait on the invariant plane s = 0 restricted to the region x2 − z 2 ≤ 1. invariant we must have ∂l ∂l ∂l ẋ + ż + ṡ = Kl, ∂x ∂z ∂s for some polynomial K, and this is the case with K = 2(3x2 − 2z 2 − √ 6xs). On the surface x2 − z 2 = 1 system (1) can be written as √ (8) dx ds 6 2 = x (x2 − 1) , = xs . dt dt n Then except for all the points on x = 0 which are equilibrium points, system (8) also admits two finite equilibrium points e1 = (1, 0) and e2 = (−1, 0). By using Theorem 2.19 of [20], we can find that both e1 and e2 are semi-hyperbolic saddle-nodes. 6 FABAO GAO1,2 , AND JAUME LLIBRE2 On the local chart U1 system (8) becomes √ ! du 6 dv (9) =u u + v2 − 1 , = v (v 2 − 1) . dt n dt √ It has two infinite equilibrium points e5 = (0, 0) and e6 = ( 6n/6, 0), where e5 is a hyperbolic stable node with eigenvalues −1 of multiplicity two, and e6 is a hyperbolic unstable saddle point with eigenvalues ±1. On the local chart U2 system (8) writes √ ! √ du 6 2 2 dv 6 (10) =u − u+u −v , =− uv. dt n dt n Let dτ = udt we obtain √ √ (11) du 6 2 2 dv 6 =− u+u −v , =− v. dτ n dτ n The origin e7 = (0, 0) on the local chart U2 is a hyperbolic stable node √ with eigenvalues − 6/n of multiplicity two. In short the global phase portraits of system (8) is shown in Figure 3. 2.4. The finite equilibrium points. It is noted that system (1) admits three finite equilibrium points p0 = (0, 0, 0) with eigenvalues (−3, −2, 0), p1 = (1, 0, 0) and p2 = (−1, 0, 0) with the same eigenvalues (6, 1, 0). Here p1 and p2 are located at the intersection of the invariant planes z = 0, s = 0 and the invariant surface x2 −z 2 = 1, corresponding to the equilibrium points e1 and e2 in subsections 2.1-2.3, respectively. The origin p0 of system (1) lies at the intersection of the invariant planes z = 0 and s = 0, which is the same point as the equilibrium point e0 studied in subsections 2.1 and 2.2. 7 s x Figure 3. The phase portrait of the invariant sur- face x2 − z 2 = 1. 3. Phase portraits on the Poincaré sphere at infinity In order to describe the dynamics of system (1) at infinity. We use the method of the three-dimensional Poincaré compactification (see [21] for more details) x = 1/z3 , z = z1 /z3 , s = z2 /z3 , then the analytical vector field of system (1) on the local chart U1 becomes dz1  √  = z1 z32 − 6z2 (−1 + z12 + z32 ) , dt    dz2 2 2 √ 1 2 2 (12) = z2 −3 + 2z1 + 3z3 + 6z2 + 1 − z1 − z3 , dt n dz3  √  √   = −z3 z12 −2 + 6z2 + −3 + 6z2 (−1 + z32 ) . dt In the different local charts of the Poincaré sphere, the infinity of R3 corresponds to z3 = 0, then system (12) has the equilibrium points listed in Table 1. Except that the equilibrium point u31 denotes the 8 FABAO GAO1,2 , AND JAUME LLIBRE2 origin of the local chart U3 , the rest equilibrium points lie in the local chart U1 . Moreover the straight line s = 0 of the local chart U1 is filled with the equilibrium points ua0 for all a ∈ R. Table 1. Equilibrium points on the different local charts of the Poincaré sphere at the infinity of R3 . Equilibrium points Eigenvalues u11 = (0, 0, 0) (−3, −3, 0) √ ! 6 u12 = − , 0, 0 (0, 0, 0) 2 √ ! 6 u13 = , 0, 0 (0, 0, 0) 2 √ ! 6 u14 = −1, n, 0 (−1, 1, −2n) 6 √ ! 6 u15 = 1, n, 0 (−1, 1, −2n) 6 √ !   6 n 3 3n u16 = 0, ,0 − , ,3 2 1+n 1+n 1+n  ua0 = (a, 0, 0) 0, −3 + 2a2 , −3 + 2a2 u31 = (0, 0, 0) (0, 0, 0) For the case z3 = 0 system (12) is reduced to dz1 √ = − 6z1 z2 (−1 + z12 ) , (13) dt    dz2 2 √ 1 2 = z2 −3 + 2z1 + 6z2 + 1 − z1 . dt n After changing the of time dτ = z2 dt system (13) becomes dz1 √ = − 6z1 (−1 + z12 ) , (14) dτ   dz2 2 √ 1 2 = −3 + 2z1 + 6z2 + 1 − z1 . dτ n 9 Then this system has three equilibrium points ei,1 , ei,2 and ei,3 with √ √ √ coordinates (−1, 6n/6), (1, 6n/6) and (0, 6n/(2(1 + n))), respec- tively. Here both ei,1 and ei,2 are unstable saddle points with eigenval- √ √ √ ues 6/n and −2 6. ei,3 is an unstable node with eigenvalues 6 and √ 6(1 + n)/n. The phase portrait on the Poincaré sphere at infinity on local chart U1 is shown in Figure 4. s z Figure 4. The phase portrait of system (1) at in- finity on local chart U1 . On the local chart U2 in view of Poincaré compactification x = z1 /z3 , z = 1/z3 , s = z2 /z3 , system (1) writes dz1 √ = −z1 z32 + 6z2 (1 − z12 + z32 ) , dt √ ! dz2 2 6 2 (15) = z2 2 − 3z1 + z1 z2 + 2z3 , dt n dz3 = z3 (2 − 3z12 + 2z32 ) . dt 10 FABAO GAO1,2 , AND JAUME LLIBRE2 Since we want to study the infinity we take z3 = 0, rescaling the time dτ = z2 dt system (15) is equivalent to dz1 √ = 6 (1 − z12 ) , (16) dτ √ dz2 2 6 = 2 − 3z1 + z1 z2 . dτ n For any constant b, since (b, 0) is not the equilibrium point of the system (16), we will not continue to discuss other infinite equilibrium points of the system (16). Because this has been discussed in the case U1 . On the local chart U3 the three-dimensional Poincaré compactifica- tion is x = z1 /z3 , z = z2 /z3 , s = 1/z3 , then system (1) becomes (17) " √ # dz1 6(1 + n) √ = z1 − z1 + 3z12 − 2z22 − 3z32 + 6 (z22 + z32 ) , dt n " √ # dz2 6 = z2 − z1 + 3z12 − 2 (z22 + z32 ) , dt n √ dz3 6 =− z1 z3 . dt n For the case z3 = 0 from system (17) we obtain " √ # dz1 6(1 + n) √ = z1 − z1 + 3z1 − 2z2 + 6z22 , 2 2 dt n (18) √ ! dz2 6 = z2 − z1 + 3z12 − 2z22 . 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 [20] for more details). In order to study its local phase portrait we shall use vertical blow-ups (see for instance [22]) by letting w = z2 /z1 , then 11 we obtain     dz1 2 √ 1 2 2 = z1 6 − − 1 + w + (3 − 2w ) z1 , (19) dt n dw √ = − 6z1 w (−1 + w2 ) . dt Eliminating the common factor z1 of system (19) by changing the time dτ = z1 dt it yields     dz1 √ 1 2 2 = z1 6 − − 1 + w + (3 − 2w ) z1 , (20) dτ n dw √ = − 6w (−1 + w2 ) . dτ System (20) admits three equilibrium points ei,4 = (0, −1), ei,5 = (0, 1) and ei,6 = (0, 0) on z1 = 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 Figure 5(a). Note that there is a time rescaling dτ = z1 dt between systems (19) and (20), so the direction of the trajectories in the local phase portraits of system (19) is opposite to that of Figure 5(a) when z1 < 0, see Figure 5(b) for more details. In addition, all points on the w axis, i.e. z1 = 0, are singularities of system (19). Thus the local phase portraits at (0, 0) of system (18) is shown in Figure 5(c). Then the phase portrait in the local chart U3 is shown in Figure 6. In summary joining the previous information we obtain the global phase portraits at infinity in the Poincaré sphere in Figure 7. 12 FABAO GAO1,2 , AND JAUME LLIBRE2 (a) (b) (c) Figure 5. In (a), (b) and (c) there are the local phase portraits of the equilibrium points in systems (20), (19) and (18), respectively. 4. Phase portraits inside the Poincaré ball restricted to x2 − z 2 ≤ 1 It is noted that system (1) is invariant under the symmetries (x, z, s) 7→ (−x, −z, −s) and (x, z, s) 7→ (−x, z, −s), so it is invariant under the 13 z x Figure 6. The phase portrait of system (18), i.e. the phase portrait at infinity in the local chart U3 of system (1). symmetry with respect to the origin and to the z-axis. Now we di- vide the Poincaré ball restricted to x2 − z 2 ≤ 1 into the following four regions: R1 : z ≤ 0, s ≥ 0. R2 : z ≤ 0, s ≤ 0. R3 : z ≥ 0, s ≥ 0. R4 : z ≥ 0, s ≤ 0. Then due to the symmetries we only need to study the phase portraits of system (1) in the region R1 . Combining the phase portraits of the invariant surface x2 − z 2 = 1 with the phase portraits of the planes z = 0, and s = 0, together with the phase portrait at infinity, we get the phase portraits on the boundary of the region R1 as shown in Figures 8-10. Here we explain the definition of the three-dimensional cartesian coordinate system in this paper as follows: we regard the xz-plane as the horizontal plane 14 FABAO GAO1,2 , AND JAUME LLIBRE2 s z x Figure 7. The global phase portrait at infinity in the Poincaré sphere. The North Pole of the Poincaré ball is the equilibrium point u31 . The symmetric points of u12 and u13 with respect to the center of the sphere are U12 and U13 , respectively. in R3 , where the direction of the z-axis is horizontal towards the right. If the z-axis axis rotates 90 degrees clockwise we get the x-axis. The s-axis is vertical upward, and then xzs constitutes a three-dimensional right-handed cartesian coordinate system. In order to show the phase portraits more clearly, the boundary of the region R1 is divided into six surfaces according to the orientation towards us, the back to us and the bottom area. See Figure 11 for more details. It is noted that the North Pole u31 of the Poincaré ball on the front boundary surfaces F1 and F2 is stable, and there is an elliptic sector and a stable parabolic sector segment of U3 on the boundary F3 , but it is unstable on the back boundary surfaces B1 and B2 . 15 s z x Figure 8. Phase portrait in the front boundary of the region R1 . s z x Figure 9. Phase portrait in the back boundary of the region R1 . 5. Dynamics in the interior of the region R1 Note that the original system (1) admits the three finite equilib- rium points p0 , p1 and p2 in the three-dimensional cartesian coordinate 16 FABAO GAO1,2 , AND JAUME LLIBRE2 s z x Figure 10. Phase portrait in the bottom of the re- gion R1 . u 31 u 31 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 Figure 11. The six boundary surfaces of the region R1 . system. The dynamical behavior of the system inside the region R1 de- pends on the behavior of the flow in the following surfaces and planes h(x, z, s) = 0, g(x, z) = 0, x = 0, z = 0, s = 0, 17 where √ h(x, z, s) = x (3x2 − 2z 2 − 3) + 6s (1 − x2 + z 2 ) , g(x, z) = 3x2 − 2z 2 − 2. These surfaces and planes divide the region R1 into seven different subregions R1i , i = (1, 2, . . . , 7), see Figures 12-19 for more details. It should be noted that R11 and R13 represent parallel tunnels inside R1 , respectively. It is easy to verify that h > 0 in the subregions 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 and R15 . It should be noted that the dotted and solid lines in Figures 13-19 are consistent with those in Figure 12. s z x Figure 12. There are seven subregions inside the region R1 of the Poincaré ball. As shown in the subregion R11 (see Figure 13) the front surface consists of three dashed lines and one solid line contained in the surface 18 FABAO GAO1,2 , AND JAUME LLIBRE2 s z x Figure 13. The subregion R11 . s z x Figure 14. The subregion R12 . h = 0, and the surface on the back side (opposite to the above front surface) contained in the surface g = 0. In Table 2 we describe the behavior of ẋ, ż and ṡ in the seven sub- regions R11 , . . . , R17 . From this table we obtain that the variables x and z in the subregion R11 decrease monotonically, and the variable s increases monotonically, so an orbit in the subregion R11 either crosses the boundary h = 0 and enters into the subregion R12 , or crosses the 19 s z x Figure 15. The subregion R13 . s z x Figure 16. The subregion R14 . boundary g = 0 and goes to the subregion R13 , or crosses the intersec- tion curve of the subregions R11 , R12 and R13 into the subregion R14 with very low probability, i.e. an orbit in the subregion R11 will not stay in the future in this region, but will exit through its boundaries into other subregions, and these orbits in backwards time come from the equilibrium point p1 in the subregion R11 . This process can be 20 FABAO GAO1,2 , AND JAUME LLIBRE2 s z x Figure 17. The subregion R15 . s z x Figure 18. The subregion R16 . simply summarized as follows R12 p1 R11 R13 . R14 21 s z x Figure 19. The subregion R17 . Table 2. Dynamical behavior in seven different subregions Subregions Corresponding Region Increase or decrease 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 The front boundary of R12 (see Figure 14), i.e. the surface containing the equilibrium points u12 , u14 , u31 and p1 , is contained in the invariant surface x2 − z 2 = 1 (x > 0). The surface on the back side of R12 that contains the equilibrium points p1 and u31 consists of two parts: the upper part surface containing the infinite equilibrium point u31 is the intersection of the subregions R12 and R14 , which is contained in the 22 FABAO GAO1,2 , AND JAUME LLIBRE2 surface g = 0, and the lower part surface containing the finite equi- librium point p1 , is the intersection of subregions 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 u14 at infin- ity. 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 p1 or from the subregions R11 and R14 which have common boundaries with the subregion R12 . This is due to the fact that the variable z is decreasing monotonically and the variables x and s are increasing monotonically inside the subregion R12 according to Table 2. Moreover a trajectory in the subregion R12 will tend to the equilibrium point u14 at infinity. In summary, this dynamic behavior process can be represented as p1 R11 R12 u14 . R14 In the subregion R13 (see Figure 15) the front surface composed of three dashed lines and one solid line is contained in the surface g = 0, and its opposite surface is contained in the surface h = 0. From Table 2 an orbit in the subregion R13 may only comes from the equilibrium points locate at the infinity on s = 0 in R13 , or it comes from the subregions R11 . Then the orbit passes through the intersection of subregions R13 and R14 (contained in the surface h = 0) and enters 23 the subregion R14 , or tends to the infinite equilibrium point u31 , i.e. lies in the North Pole of the Poincaré sphere. So we obtain the following dynamics s = 0 in R13 u31 R13 R11 R14 . In subregion R14 (see Figure 16) the front triangular surface is con- tained in g = 0, and its opposite triangular surface is on the back side, which is included in the plane x = 0. According to Table 2 an orbit of the subregion R14 may come from the infinite equilibrium points on x = 0, or come from the subregion R13 through the surface h = 0, or from subregion R15 traversing plane x = 0, then go through the surface g = 0 into the subregion R12 , or directly goes to the infinite equilib- rium point u31 at the North Pole of the Poincaré sphere. This can be represented as R13 R12 R15 R14 x = 0 in R14 u31 . It can be noted from Figure 17 that the surface on the left of subre- gion R15 (consisting of two dashed lines and one solid line) is contained in plane x = 0, and the opposite surface to it is included in the surface g = 0. According to Table 2 an orbit in the subregion R15 can only come from the adjacent subregion R16 , and then either enter into the 24 FABAO GAO1,2 , AND JAUME LLIBRE2 subregion R14 adjacent to it, or directly tend to the finite equilibrium point p0 located at the center of the Poincaré ball. This dynamics can be denoted as R14 R16 R15 p0 . We observe that the curved surfaces of the left and right sides of the subregion R16 (see Figure 18) are contained in the surfaces g = 0 and h = 0, respectively. It is known from Table 2 that an orbit in the subregion R16 can only come from the subregion R17 and crossing the surface g = 0, it enters into the subregion R15 . Thus we get that R17 R16 R15 . The left and right surfaces of subregion R17 (see Figure 19) are in- cluded in the surface h = 0 and in the invariant surface x2 − z 2 = 1 (x < 0), respectively. According to Table 2, it is known that all three variables in subregion R17 are decreasing monotonically, so an orbit in this subregion must come from the infinite equilibrium point u31 located at the North Pole of Poincaré sphere and eventually tends to the infinite equilibrium point U13 in the invariant plane s = 0. Then we have that u31 R17 U13 . 25 The dynamic behavior of the orbits inside the seven subregions of R1 discussed above can be summarized as p0 x = 0 in R14 R16 R15 R14 R12 u14 U13 R17 u31 R13 R11 p1 . s = 0 in R13 The above flow chart shows that the orbits of system (1) contained in the interior of the region R1 have α-limit at the finite equilibrium point p1 , and the some orbits on the boundary of the region R1 have α-limit at the x = 0 in R14 , or s = 0 in R13 . Moreover the orbits have ω-limit either at the finite equilibrium point p0 , or at the infinite equilibrium points u14 and U13 , where u14 is located on the intersection curve of the Poincaré sphere and the invariant surface x2 − z 2 = 1 (x > 0) at infinity, and U13 is located on the intersection of the Poincaré sphere and the invariant surfaces s = 0, x2 − z 2 = 1 (x < 0) at infinity (see Figure 12). Furthermore the orbits of system (1) on the boundary of R1 also have α-limit and ω-limit at the infinite equilibrium point u31 , i.e. the North Pole of the Poincaré sphere. In this way the qualitative global dynamic behavior of system (1) is described. 26 FABAO GAO1,2 , AND JAUME LLIBRE2 6. Discussion and conclusions The physical meaning of the variables x, z and s in this paper can be found in references [3, 7, 8, 9, 14]. So we have φ̇ µ V 0 (φ) (21) x= √ , z= s=− , 2 6H 4(3λ − 1)a2 H V (φ) where λ is a dimensionless constant, V (φ) is the potential, µ is a con- stant, and H = ȧ/a is the Hubble parameter, and a is the dimensionless scale factor for the expanding universe. Near infinity we note that the phase portrait described in Figure 2 has been obtained in [9], where the phase portrait of system (1) is limited to s = 0. Now we have completely described the phase portrait of system (1) with all s values in the region G as shown in Figures 8-10 and 12. By using the fact that system (1) remains invariant under two aforementioned symmetries, we provide the complete global phase portrait of system (1) in G. The phase portrait shows that the final evolution of the orbital direction of system (1) within G tends to the center of the Poincaré ball as well as the infinite equilibrium points u14 and U13 located at the waist and equator of the Poincaré sphere (except for the initial conditions s = 0), respectively. Based on Hořava- Lifshitz gravity in a Friedmann-Lemaı̂tre-Robertson-Walker space-time with non-zero curvature and without the cosmological constant term, equations (21) imply that the value of Hubble parameter H is tending to be zero in the cosmological model as time goes on. 27 Acknowledgments The first author gratefully acknowledges the support of the Na- tional Natural Science Foundation of China (NSFC) through grant Nos.11672259 and 11571301, the China Scholarship Council through grant No.201908320086, the Ministry of Land and Resources Research of China in the Public Interest through grant No.201411007. The second author gratefully acknowledges the support of the Min- isterio de Economı́a, Industria y Competitividad, Agencia Estatal de Investigación grants MTM2016-77278-P (FEDER) and MDM-2014- 0445, the Agència de Gestió d’Ajuts Universitaris i de Recerca grant 2017SGR1617, and the H2020 European Research Council grant MSCA- RISE-2017-777911. References [1] P. Hořava, Quantum gravity at a Lifshitz point, Physical Review D 79, 084008, 2009. [2] E.M.C. Abreu, A.C.R. Mendes, G. Oliveira-Neto et al., Hořava-Lifshitz cosmo- logical models with noncommutative phase space variables, General Relativity and Gravitation 51, 95, 2019. [3] S. Carloni, E. Elizalde and P.J. Silva, An analysis of the phase space of Hořava-Lifshitz cosmologies. In: S.D. Odintsov, D. Sáez-Gómez and S. Xambó- Descamps (eds) Cosmology, Quantum Vacuum and Zeta Functions, Springer Proceedings in Physics 137, 139-148, Springer-Verlag, Berlin, 2011. [4] B. Chen, On Hořava-Lifshitz cosmology, Chinese Physics C 35 (5), 429-435, 2011. [5] F.B. Gao and J. Llibre, Global dynamics of the Hořava-Lifshitz cosmological system, General Relativity and Gravitation, Accepted, 2019. [6] X. Gao,Y. Wang, R. Brandenberger and A. Riotto, Cosmological perturbations in Hořava-Lifshitz gravity, Physical Review D 81, 083508, 2010. 28 FABAO GAO1,2 , AND JAUME LLIBRE2 [7] G. Leon and C.R. Fadragas, Cosmological dynamical systems: and their ap- plications, Lambert Academic Publishing, GmbH & Co. KG, 2012. [8] G. Leon and A. Paliathanasis, Extended phase-space analysis of the Hořava- Lifshitz cosmology, The European Physical Journal C 79, 746, 2019. [9] G. Leon and E.N. Saridakis, Phase-space analysis of Hořava-Lifshitz cosmology, Journal of Cosmology and Astroparticle Physics 2009, 006, 2009. [10] S. Lepe and J. Saavedra, On Hořava-Lifshitz cosmology, Astrophysics and Space Science 350, 839-843, 2014. [11] M. Li and Y. Pang, A trouble with Hořava-Lifshitz gravity, Journal of High Energy Physics 08, 015, 2009. [12] O. Luongo, M. Muccino and H. Quevedo, Kinematic and statistical inconsis- tencies of Hořava-Lifshitz cosmology, Physics of the Dark Universe 25, 100313, 2019. [13] N.A. Nilsson and E. Czuchry, Hořava-Lifshitz cosmology in light of new data, Physics of the Dark Universe 23, 100253, 2019. [14] A. Paliathanasis and G. Leon, Cosmological solutions in Hořava-Lifshitz grav- ity, arXiv preprint arXiv: 1903.10821, 2019. [15] E.N. Saridakis, Aspects of Hořava-Lifshitz cosmology, International Journal of Modern Physics D 20 (08), 1485-1504, 2011. [16] A. Tawfik and E. Abou El Dahab, FLRW cosmology with Hořava-Lifshitz gravity: impacts of equations of state, International Journal of Theoretical Physics 56 7, 2122-2139, 2017. [17] M. Bhattacharjee, Gravitational radiation and black hole formation from grav- itational collapse in theories of gravity with broken Lorentz symmetry, Baylor University, ProQuest Dissertations Publishing 22585106, 2019. [18] S. Mukohyama, Hořava-Lifshitz cosmology: a review, Classical and Quantum Gravity 27, 223101, 2010. [19] T.P. Sotiriou, Hořava-Lifshitz gravity: a status report, Journal of Physics: Conference Series 283, 012034, 2011. [20] F. Dumortier, J. Llibre and J.C. Artés, Qualitative theory of planar differential systems, Springer-Verlag, Berlin, 2006. 29 [21] A. Cima and J. Llibre, Bounded polynomial vector fields, Transactions of the American Mathematical Society 318(2), 557-579, 1990. [22] M.J. Álvarez, A. Ferragut and X. Jarque, A survey on the blow up technique, International Journal of Bifurcation and Chaos 21(11), 3103-3118, 2011. 1 School of Mathematical Science, Yangzhou University, Yangzhou 225002, China E-mail: gaofabao@sina.com (Fabao Gao, ORCID 0000-0003-2933-1017) 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Bellaterra 08193, Barcelona, Catalonia, Spain E-mail: jllibre@mat.uab.cat (Jaume Llibre, ORCID 0000-0002-9511-5999)