This is a preprint of: “Global dynamics of the Horava-Lifshitz cosmological system”, Fabao Gao, Jaume Llibre, General Relativity and Gravitation, vol. 51, 152:1–15, 2019. DOI: [10.1007/s10714-019-2635-x] GLOBAL DYNAMICS OF THE HOŘAVA-LIFSHITZ COSMOLOGICAL SYSTEM FABAO GAO1,2 , AND JAUME LLIBRE2 Abstract. Using the qualitative theory, the global dynamics of the cosmological model based on Hořava-Lifshitz gravity is de- scribed. 1. Introduction Ten years ago Hořava [1] put forward a theory of spacetime asym- metric gravity, which is similar to Lifshitz’s scalar field theory. If the spatial dimension in Lifshitz’s scalar field theory has a weight of one, then the time dimension has a weight of three. Therefore this theory is also known as Hořava-Lifshitz gravity. It ignited a great deal of re- search on the possible application of this theory in cosmology and black hole physics (see [2]-[5] or the review articles [6], [7] and the references therein). With or without detailed-balance conditions, Leon and Saridakis [8] carried out a detailed phase space analysis of Hořava-Lifshitz cosmol- ogy, and found that the universe governed by Hořava gravity had late- time solutions compatible with observations. They also presented sev- eral results on the stability of de Sitter solutions in Hořava-Lifshitz cosmology by using the central manifold theory [9]. Furthermore Sari- dakis [10] reviewed some general aspects of Hořava-Lifshitz cosmology and extracted cosmological equations from its basic version. He proved Key words and phrases. Hořava-Lifshitz, Global dynamics, cosmology. 1 2 FABAO GAO1,2 , AND JAUME LLIBRE2 that bouncing and cycling can occur naturally by phase space analy- sis. However, he also showed that Hořava-Lifshitz gravity is affected by instability in its basic version with detailed perturbation analysis. For the stability of gravitational scalar mode Chen [11] briefly summa- rized the Hořava-Lifshitz theory of gravity and its modifications and its implications in cosmology. In recent years Lepe and Saavedra [12] discussed some aspects of Hořava-Lifshitz cosmology with emphasis on some cosmological solu- tions that exist in general relativity (Friedmann cosmology), especially in the flat case that dust-driven evolution is the same in both cos- mological theories. For Hořava-Lifshitz theory of gravitation Abreu et al. [13] explored a non-commutative version of the Friedmann- Robertson-Walker cosmological model, in which material content is de- scribed by ideal fluids, and the constant curvature of the spatial sections can be positive, negative or zero. Under the background spacetime of Friedmann-Lemaı̂tre-Robertson-Walker, Paliathanasis and Leon [14] divided the integrability of Hořava-Lifshitz scalar field into four cases according to the existence of cosmological constant term and the dis- appearance of space curvature. They followed the singularity analysis method to determine integrability. It is believed that their work will be very helpful to the integrability of the gravitational field equation in cosmology. In addition to the cosmological solutions in [14], two versions of Hořava-Lifshitz gravity were discussed for finding and an- alyzing the plane symmetric, static (non-static) solutions in Hořava- Lifshitz gravity [15]. They also studied the plane symmetry of the new modified version of Hořava gravity. 3 Here we shall describe the global dynamics of the Hořava-Lifshitz cosmological model in a Friedmann-Lemaı̂tre-Robertson-Walker space- time with zero curvature and without the cosmological constant term. This dynamics is provided by the gravitational field equations in di- mensionless variables given by dy1 √  = (y12 − 1) 3y1 − 6y3 , dt dy2 (1) = (3y12 − 2) y2 , dt dy3 √ = −2 6y1 f (y3 ), dt with the power law potential f (y3 ) = −y32 /(2n), where n is a natural number. See equations (17)-(19) of [14] for more details. The description of dynamics of system (1) is given in section 4. 2. Phase portraits on the invariant planes In order to study the phase portraits of system (1), we start studying the phase portraits of its invariant planes y1 = ±1, y2 = 0, y3 = 0. After we will study the local phase portraits of the finite and infinite equilibrium points, and finally the global phase portraits in the region −1 ≤ y1 ≤ 1, which is the interest region for cosmology. 2.1. The invariant plane y1 = 1. On this plane system (1) becomes √ (2) dy 2 dy 3 6 2 = y2 , = y . dt dt n 3 4 FABAO GAO1,2 , AND JAUME LLIBRE2 The unique finite equilibrium point is qy1 ,1 = (0, 0). It is a semi- hyperbolic equilibrium point and using Theorem 2.19 of [16] qy1 ,1 is a saddle-node. Based on the Poincaré transformation y2 = 1/v, y3 = u/v, on the local chart U1 (see for more details on the Poincaré compactification Chapter 5 of [16]) equations (2) become √ (3) 6 2 u̇ = u − uv, v̇ = −v. n This system has the infinite semi-hyperbolic equilibrium point py1 ,1 = (0, 0). Applying to it Theorem 2.19 of [16] we obtain the saddle-node shown in Figure 1(a). Similarly on the local chart U2 system (2) is √ √ (4) 6 6 u̇ = − u + uv, v̇ = − v. n n Then the equilibrium point py1 ,2 = (0, 0) of system (4) is a hyperbolic √ stable node with eigenvalues − 6/n of multiplicity two. Finally joining the previous information on the studied three equi- librium points qy1 ,1 , py1 ,1 and py1 ,2 , together with the diametrically opposite equilibrium points Py1 ,1 and Py1 ,2 at infinity of py1 ,1 and py1 ,2 , we obtain the global phase portrait of system (2) in the Poincaré disc of the invariant plane y1 = 1 in Figure 1(b). Here we have also used that the straight lines y2 = 0 and y3 = 0 are invariant by the flow of system (2). 5 py , 1 2 Py , 1 1 py ,1 1 qy , 1 1 y3 y2 Py , 1 2 (a) (b) Figure 1. In (a) there is the local phase portraits of the saddle-node py1 ,1 = (0, 0) in the Poincaré disc. In (b) there is the phase portraits of the invariant plane y1 = 1. 2.2. The invariant plane y1 = −1. On this plane system (1) reduces to √ (5) dy2 dy3 6 2 = y2 , =− y . dt dt n 3 System (5) is similar to system (2), so the equilibrium point qy1 ,2 = (0, 0) is a saddle-node. On the local chart U1 system (5) is √ (6) 6 2 u̇ = − u − uv, v̇ = −v. n The semi-hyperbolic equilibrium point py1 ,3 = (0, 0) is a saddle-node by Theorem 2.19 of [16], and its local phase portrait in the Poincaré disc is described in Figure 2(a). On the local chart U2 system (5) becomes √ √ (7) 6 6 u̇ = u + uv, v̇ = v. n n 6 FABAO GAO1,2 , AND JAUME LLIBRE2 py , 1 4 qy , 1 2 Py , 1 3 py ,1 3 y3 y2 Py , 1 4 (a) (b) Figure 2. In (a) there is the local phase portraits of the saddle-node py1 ,3 = (0, 0) in the Poincaré disc. In (b) there is the phase portraits of the invariant plane y1 = −1. Its equilibrium point py1 ,4 = (0, 0) is a hyperbolic unstable node with √ eigenvalue 6/n of multiplicity two. In short joining the previous information as we did for the plane y1 = 1, we get the global phase portraits in the Poincaré disc of the invariant plane y1 = −1 in Figure 2(b). Furthermore by using the method of Poincaré compactification y1 = z1 /z3 , y2 = 1/z3 , y3 = z2 /z3 and the change of time dt = z33 dτ in R3 , then the analytical vector field of system (1) in the local chart U2 becomes √ √ z10 = − 6z12 z2 − z1 z32 + 6z2 z32 , √ ! 6 (8) z20 = z2 −3z12 + z1 z2 + 2z32 (−1 + z32 ), n z30 = −3z12 z3 + 2z33 , where the prime denotes derivative with respect to the time τ . Then this system has the equilibrium points (z1 , z2 , z3 ) = (0, a, 0) for all a. So the infinity of Figure 3 is filled up with equilibrium points when y1 = 0. 7 Therefore combining Figures 1(b) and 2(b) we obtain the phase por- trait of system (1) inside the Poincaré ball as shown in Figure 3. All the points on the dashed circle in the plane y1 = 0 are the equilibrium points. (1 ,0 ,0 ) (-1,0,0) y3 y1 y2 Figure 3. Phase portraits on the planes y1 = ±1 inside the Poincaré ball. 2.3. The invariant plane y2 = 0. On this plane system (1) can be rewritten as √ dy1 √  dy3 6 (9) 2 = (y1 − 1) 3y1 − 6y3 , = y1 y32 . dt dt n This system has three equilibrium points Oy2 ,0 = (0, 0), qy2 ,1 = (1, 0), qy2 ,2 = (−1, 0), which are semi-hyperbolic. In order to apply Theorem 2.19 of [16], the linear part at each equi- librium point of system (9) must be written into its real Jordan normal 8 FABAO GAO1,2 , AND JAUME LLIBRE2 form. So we do the following transformation to the linear part of Oy2 ,0 √ 6 y1 = (y − x), y3 = −x. 3 Then equations (9) become 2 2 ẋ = − x3 + x2 y, (10) n n   2 3 1 1 ẏ = 3y − x + − 1 x2 y + 4xy 2 − 2y 3 . n 2 n Let y = g(x) = 2x3 /3n + 4(n − 1)x5 /(9n2 ) + O(x5 ) be the solution of the second equation of system (10). Substituting y = g(x) into the first equation of system (10), we obtain ẋ = −2x3 /n + O(x5 ). So by Theorem 2.19 of [16], Oy2 ,0 = (0, 0) is a saddle point. Applying again Theorem 2.19 of [16] to the equilibrium points qy2 ,1 = (1, 0) and qy2 ,2 = (−1, 0), we obtain that they are saddle-nodes. On the local chart U1 system (9) becomes " √ # 6 u̇ = u 3 (v 2 − 1) + u (1 + n − nv 2 ) , n (11) √  v̇ = v 3 − 6u (v 2 − 1) . √  Its equilibrium points are py2 ,1 = (0, 0), py2 ,2 = 6n/(2(1 + n)), 0 . The equilibrium point py2 ,1 is a hyperbolic stable node with eigenvalues −3 of multiplicity two, and py2 ,2 is a hyperbolic unstable saddle with eigenvalues 3 and −3 + 3n/(1 + n). On the local chart U2 system (9) writes √ √ 2 √  6(1 + n) 2 3 6 (12) u̇ = v 6 − 3u − u + 3u , v̇ = − uv. n n 9 The py2 ,3 = (0, 0) is a linearly zero equilibrium point, i.e. its linear part is identically zero. Its topological index is zero by the Poincaré-Hopf theory (see Theorem 6.30 of [16] for more details). In order to study the local phase portraits of the equilibrium point (0, 0) of system (12) we apply the blow-up techniques (see [17] for more details). We do a vertical blow-up by introducing the transformation w = v/u, then we have " √ # 6 √  u̇ = u2 − + 6 − 3u (w2 − 1) , n (13) √  ẇ = −uw 6 − 3u (w2 − 1) . The common factor u of system (13) can be eliminated by rescaling the time udt = dτ . So we get " √ # 6 √  u0 = u − + 6 − 3u (w2 − 1) , n (14) √  w0 = −w 6 − 3u (w2 − 1) , where the prime represents the derivative with respect to the time τ . Note that system (14) has three equilibrium points py2 ,4 = (0, −1), py2 ,5 = (0, 0) and py2 ,6 = (0, 1) on u = 0. Therefore both the equilib- rium points py2 ,4 and py2 ,6 are hyperbolic stable nodes with eigenvalues √ √ − 6/n and −2 6, the equilibrium point py2 ,5 is a hyperbolic unstable √ √ saddle with eigenvalues − 6(1 + n)/n and 6. The local phase por- traits around these three equilibrium points are shown in Figure 4(a). Considering that there is a time rescaling dτ = udt between systems (13) and (14), so the local phase portraits of system (13) is shown in 10 FABAO GAO1,2 , AND JAUME LLIBRE2 Figure 4(b). Moreover, all points on the w-axis, i.e. u = 0, are singu- larities of system (13). Thus the local phase portraits of system (12) is shown in Figure 4(c), and then the local phase portrait at the origins of U2 and V2 for y2 = 0 can be found in Figure 4(d). In summary joining the previous information we obtain the global phase portraits in Figure 5 in the Poincaré disc of the plane y2 = 0 restricted to the strip −1 ≤ y1 ≤ 1. 2.4. The invariant plane y3 = 0. On this plane system (1) writes dy1 dy2 (15) = 3y1 (y12 − 1) , = y2 (3y12 − 2) . dt dt This system has three equilibrium points Oy3 ,0 = (0, 0), qy3 ,1 = (1, 0) and qy3 ,2 = (−1, 0). Then Oy3 ,0 is a hyperbolic stable node with eigen- values -3 and -2, both qy3 ,1 and qy3 ,2 are unstable hyperbolic nodes with eigenvalues 1 and 6. On the local chart U1 system (15) becomes (16) u̇ = uv 2 , v̇ = 3v (v 2 − 1) . So all the infinity of this system is filled of equilibrium points. We note that removing v of system (16) doing the change of time dτ = vdt, we get the system (17) u0 = uv, v 0 = 3 (v 2 − 1) , which has no infinite equilibrium points. On the local chart U2 system (15) writes (18) u̇ = −uv 2 , v̇ = v (−3u2 + 2v 2 ) . 11 w w (0,1) (0,1) (0,0) (0,0) u u (0,-1) (0,-1) (a) (b) py , 2 3 v (0,0) u Py , 2 3 (c) (d) Figure 4. In (a), (b) and (c) there are the local phase portraits of the equilibrium points in system (14), (13) and (12), respectively. In (d) there is the local phase portraits at the origin of U2 and V2 for y2 = 0. Again doing the change of variable dτ = vdt system (18) becomes (19) u0 = −uv, v 0 = −3u2 + 2v 2 . 12 FABAO GAO1,2 , AND JAUME LLIBRE2 -1 1 y3 y1 Figure 5. The phase portraits of the invariant plane y2 = 0 restricted to −1 ≤ y1 ≤ 1. The py3 ,1 = (0, 0) is also a linearly zero equilibrium point of system (19). Applying again the blow-up techniques of [17] to py3 ,1 , then we obtain (20) u0 = −u2 w, w0 = 3u(w2 − 1), after doing the change of variable w = v/u. Rescaling the time of system (20) we get (21) u̇ = −uw, ẇ = 3(w2 − 1). 13 Hence system (21) has two equilibrium points py3 ,2 = (0, −1) and py3 ,3 = (0, 1) on u = 0. Therefore both the equilibrium points py3 ,2 and py3 ,3 are hyperbolic unstable saddles with eigenvalues 1, −6 and −1, 6, respectively. The local phase portraits around these two equi- librium points is shown in Figure 6(a). Considering that there is time dτ = udt to rescale between systems (20) and (21), the local phase por- traits of system (20) is shown in Figure 6(b). In addition, all points on the w-axis, i.e. u = 0 are the singularities of the system (20). Therefore the local phase portraits of system (19) is shown in Figure 6(c), and then the local phase portraits at the origins of U2 and V2 of y3 = 0 can be found in Figure 6(d). In conclusion from the previous information and taking into account that the straight lines y1 = 0 and y2 = 0 are invariant under the flow of system (15), we obtain the global phase portraits restricted to the band −1 ≤ y1 ≤ 1 in Figure 7 in the Poincare disc of the plane y3 = 0. 3. Phase portraits inside the Poincaré ball restricted to −1 ≤ y1 ≤ 1 We divide the Poincaré ball restricted to −1 ≤ y1 ≤ 1 into four regions: R1 : y2 ≤ 0, y3 ≥ 0. R2 : y2 ≤ 0, y3 ≤ 0. R3 : y2 ≥ 0, y3 ≥ 0. R4 : y2 ≥ 0, y3 ≤ 0. Since system (1) is invariant under the symmetry with respect to the y2 -axis, i.e. (y1 , y2 , y3 ) = (−y1 , y2 , −y3 ), we only need to study the phase portraits in the regions R1 and R2 . 14 FABAO GAO1,2 , AND JAUME LLIBRE2 w w (0,1) (0,1) u u (0,-1) (0,-1) (a) (b) py , 3 1 v u Py , 3 1 (c) (d) Figure 6. In (a), (b) and (c) there are the local phase portraits of the equilibrium points in system (21), (20) and (19), respectively. In (d) there is the local phase portraits at the origins of U2 and V2 for y3 = 0. Putting together the phase portraits of the invariant planes y1 = ±1, y2 = 0 and y3 = 0 we obtain the phase portraits in the boundary of 15 -1 1 y2 y1 Figure 7. The phase portrait of the invariant plane y3 = 0 restricted to −1 ≤ y1 ≤ 1. the regions R1 and R2 in Figures 8 and 9, respectively. The three- dimensional cartesian coordinate system in this paper is defined as follows: we consider the y1 y2 -plane as the horizontal plane in R3 , in which the direction of the y2 -axis is horizontally to the right. If the y2 -axis is rotated 90 degrees counterclockwise we get the y1 -axis. The y3 -axis is vertically upward, then y1 y2 y3 constitutes a three dimensional left-handed cartesian coordinate system. These five planes divide the regions R1 and R2 into twelve differ- ent subregions Gi , i = (1, 2, · · · , 12) (see Figures 10 and 11 for more details). 16 FABAO GAO1,2 , AND JAUME LLIBRE2 N (1,0,0) (0,0,0) y3 y1 y2 (-1,0,0) Figure 8. Phase portrait in the boundary of the region R1 . N denotes the North Pole of the Poincaré ball. Note that the original system (1) admits three finite equilibrium points (−1, 0, 0), (0, 0, 0) and (1, 0, 0). The dynamical behavior of the system inside the regions R1 and R2 depends on the behavior of the flow in the following five planes y1 = ay3 , y1 = ±a, y1 = 0, y3 = 0, √ where a = 6/3. 4. Dynamics in the interior of the regions R1 and R2 In Table 1 we describe the behavior of y˙1 , y˙2 and y˙3 in the twelve subregions Gi for i = 1, 2, · · · , 12. From this table we get that the orbits from the subregion G1 must go to the equilibrium point N , this 17 y3 y1 y2 S Figure 9. Phase portrait in the boundary of the region R2 . S denotes the South Pole of the Poincaré ball. N G1 y 1 =a y 3 G3 G2 G5 G4 (1,0,0) a G6 (0,0,0) y3 y1 -a y2 (-1,0,0) Figure 10. The six subregions of R1 in the Poincaré ball. 18 FABAO GAO1,2 , AND JAUME LLIBRE2 a G 12 -aG9 G 11 G8 G 10 y 1= ay 3 G 7 y3 y1 y2 S Figure 11. The six subregions of R2 in the Poincaré ball. is represented as follows G1 N. In a similar way and taking into account that y2 = 0 and y3 = 0 are invariant planes, we obtain that G1 N G2 and G6 G5 G3 N, G4 G3 N G7 S G8 and G12 G11 G10 S. G9 G10 S 19 Table 1. Dynamical behavior in twelve different subregions Subregions Corresponding Region Increase or decrease √ G1 y1 < y3 , 6/3 < y1 < 1, y3 > 0 y˙1 > 0, y˙2 < 0, y˙3 > 0 √ G2 y1 > y3 , 6/3 < y1 < 1, y3 > 0 y˙1 < 0, y˙2 < 0, y˙3 > 0 √ G3 y1 < y3 , 0 < y1 < 6/3, y3 > 0 y˙1 > 0, y˙2 > 0, y˙3 > 0 √ G4 y1 > y3 , 0 < y1 < 6/3, y3 > 0 y˙1 < 0, y˙2 > 0, y˙3 > 0 √ G5 y1 < y3 , − 6/3 < y1 < 0, y3 > 0 y˙1 > 0, y˙2 > 0, y˙3 < 0 √ G6 y1 < y3 , −1 < y1 < − 6/3, y3 > 0 y˙1 > 0, y˙2 > 0, y˙3 < 0 √ G7 y1 > y3 , −1 < y1 < − 6/3, y3 < 0 y˙1 < 0, y˙2 < 0, y˙3 < 0 √ G8 y1 < y3 , −1 < y1 < − 6/3, y3 < 0 y˙1 > 0, y˙2 < 0, y˙3 < 0 √ G9 y1 < y3 , − 6/3 < y1 < 0, y3 < 0 y˙1 > 0, y˙2 > 0, y˙3 < 0 √ G10 y1 > y3 , − 6/3 < y1 < 0, y3 < 0 y˙1 < 0, y˙2 > 0, y˙3 < 0 √ G11 y1 > y3 , 0 < y1 < 6/3, y3 < 0 y˙1 < 0, y˙2 > 0, y˙3 > 0 √ G12 y1 > y3 , 6/3 < y1 < 1, y3 < 0 y˙1 < 0, y˙2 < 0, y˙3 > 0 If we define the subregions R1+ = {(y1 , y2 , y3 ) ∈ R1 : 0 < y1 < 1}, R1− = {(y1 , y2 , y3 ) ∈ R1 : −1 < y1 < 0}, R2+ = {(y1 , y2 , y3 ) ∈ R2 : 0 < y1 < 1}, R2− = {(y1 , y2 , y3 ) ∈ R2 : −1 < y1 < 0}, the orbits obtained from Table 1 says that the orbits of system (1) con- tained in the region R1+ have α-limit at the equilibrium point (1, 0, 0) and ω-limit at the infinite equilibrium point N ; the orbits of system (1) contained in the region R1− have α-limit and ω-limit at the infinite equilibrium point N ; the orbits of system (1) contained in the region R2+ have α-limit and ω-limit at the infinite equilibrium point S; finally the orbits of system (1) contained in the region R2− have α-limit at 20 FABAO GAO1,2 , AND JAUME LLIBRE2 the equilibrium point (−1, 0, 0) and ω-limit at the infinite equilibrium point S. This completes the description of all the qualitative dynamics of the system (1). 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] 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. [3] M. Li and Y. Pang, A trouble with Hořava-Lifshitz gravity, Journal of High Energy Physics 08, 015, 2009. 21 [4] 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. [5] N.A. Nilsson and E. Czuchry, Hořava-Lifshitz cosmology in light of new data, Physics of the Dark Universe 23, 100253, 2019. [6] S. Mukohyama, Hořava-Lifshitz cosmology: a review, Classical and Quantum Gravity 27, 223101, 2010. [7] T.P. Sotiriou, Hořava-Lifshitz gravity: a status report, Journal of Physics: Conference Series 283, 012034, 2011. [8] G. Leon, and E.N. Saridakis, Phase-space analysis of Hořava-Lifshitz cosmol- ogy, Journal of Cosmology and Astroparticle Physics 2009, 006, 2009. [9] G. Leon, and A. Paliathanasis, Extended phase-space analysis of the Hořava- Lifshitz cosmology, arXiv preprint arXiv: 1902.09961v2, 2019. [10] E.N. Saridakis, Aspects of Hořava-Lifshitz cosmology, International Journal of Modern Physics D 20 (08), 1485-1504, 2011. [11] B. Chen, On Hořava-Lifshitz cosmology, Chinese Physics C 35 (5), 429-435, 2011. [12] S. Lepe, J. Saavedra, On Hořava-Lifshitz cosmology, Astrophysics and Space Science 350, 839-843, 2014. [13] 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. [14] A. Paliathanasis, and G. Leon, Cosmological solutions in Hořava-Lifshitz grav- ity, arXiv preprint arXiv: 1903.10821, 2019. [15] M.R. Setare and D. Momeni, Plane symmetric solutions in Hořava-Lifshitz theory, International Journal of Modern Physics D 19(13), 2079-2094, 2010. [16] F. Dumortier, J. Llibre, and J.C. Artés, Qualitative theory of planar differential systems, Springer-Verlag, Berlin, 2006. [17] 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. 22 FABAO GAO1,2 , AND JAUME LLIBRE2 1 School of Mathematical Science, Yangzhou University, Yangzhou 225002, China E-mail: gaofabao@sina.com (Fabao Gao) 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Bellaterra 08193, Barcelona, Catalonia, Spain E-mail: jllibre@mat.uab.cat (Jaume Llibre)