This is a preprint of: “Global dynamics of the Horava-Lifshitz cosmological model in a non-flat universe with non-zero cosmological constant”, Fabao Gao, Jaume Llibre, Universe (Basel), vol. 7(445), 1–50, 2021. DOI: [10.3390/universe7110445] GLOBAL DYNAMICS OF THE HOŘAVA-LIFSHITZ COSMOLOGICAL MODEL IN A NON-FLAT UNIVERSE WITH NON-ZERO COSMOLOGICAL CONSTANT FABAO GAO1,2 AND JAUME LLIBRE2 Abstract. In the case of non-zero cosmological constant the global dynamics of the cosmological model based on Hořava-Lifshitz grav- ity in a non-flat universe is characterized by using the qualitative theory of differential equations. Contents 1. Introduction 3 2. The cosmological equations 5 3. Case I: s 6= 0, k = 1 9 3.1. Phase portraits on the invariant planes and surface 9 3.1.1. The invariant plane z = 0 9 3.1.2. The invariant plane u = 0 12 3.1.3. The invariant surface f+ (x, z, u) = 0 17 3.1.4. The finite equilibrium points 23 3.1.5. Phase portrait on the Poincaré sphere at infinity 25 Key words and phrases. Global dynamics; Hořava-Lifshitz; non-flat universe; cos- mology constant. 1 2 F. GAO AND J. LLIBRE 3.2. Phase portrait inside the Poincaré ball restricted to the physical region of interest x2 − (u − z)2 ≤ 1 30 3.3. Dynamics in the interior of the regions R1 and R2 37 4. Case II: s 6= 0, k = −1 53 4.1. Dynamics in the interior of the regions S1 and S2 53 5. Case III: s = 0, k = 1 66 5.1. Phase portraits on the invariant planes and surface 66 5.1.1. The invariant plane x = 0 66 5.1.2. The invariant plane z = 0 69 5.1.3. The invariant plane u = 0 69 5.1.4. The invariant surface f+ (x, z, u) = 0 70 5.1.5. The finite equilibrium points 70 5.1.6. Phase portrait on the Poincaré sphere at infinity 71 5.2. Phase portrait inside the Poincaré ball restricted to the physical region of interest x2 − (u − z)2 ≤ 1 71 5.3. Dynamics in the interior of the region Q1 72 6. Case IV: s = 0, k = −1 80 6.1. Dynamics in the interior of the region T1 80 7. Conclusions 86 Appendix A 88 Acknowledgments 89 3 References 90 1. Introduction After the Newtonian era, Einstein put forward the general theory of relativity, and this made our understanding of gravity once again a huge leap forward. In 2009 Hořava proposed a non-relativistic theory of renormalizable gravity [1], which can be reduced to Einstein’s gen- eral theory of relativity on a large scale. It is named Hořava-Lifshitz gravity together with the scalar field theory of Lifshitz. This theory has inspired many studies and applications in cosmology [2]-[25], dark en- ergy [26]-[30], black hole [31]-[36], and electromagnetics [35]-[38]. More information can also be found from the review articles [39], [40] and the references therein. In the past ten years Leon et al. [7]-[10] have conducted several excellent studies on the Hořava-Lifshitz cosmological model whether the curvature k of the universe is zero and whether the cosmologi- cal constant Λ is considered. They divided the cosmological model into four types in the Friedmann-Lemaı̂tre-Robertson-Walker (FLRW) background spacetime: (1) Λ = 0, k = 0; (2) Λ 6= 0, k = 0; (3) Λ = 0, k 6= 0; (4) Λ 6= 0, k 6= 0. By using the phase-space analysis, they either discussed the two-dimensional dynamics of the Hořava-Lifshitz cosmological model under the usual exponential potential, and partially studied its three-dimensional dynamics. 4 F. GAO AND J. LLIBRE For the important cosmological constant Λ that many researchers have been paying attention to, this constant put the Hořava-Lifshitz gravitational theory with detailed balance, which led to a conflict be- tween its cosmology and observations. Appignani et al. [41] showed that the huge difference between the standard predictions from quan- tum field theory and the observed value of Λ may have a solution in the Hořava-Lifshitz gravity framework. Akarsu et al. [42] investigated the Λ in the standard cold-dark matter model by introducing graduated dark energy. Their results provided a high probability that the sign of Λ could be spontaneously converted, and inferred that the universe had transformed from anti-de Sitter vacua to de Sitter vacua and triggered late acceleration. Carlip [43] proposed that the vacuum fluctuations produce a huge Λ and produce a high curvature k on the Planck scale under the standard effective field theory. Although the debate about the shape of the universe has not yet reached an agreement or the boundaries of the universe are blurred, Valentino et al. [44] believed that the universe can be a closed three-dimensional sphere compared with the prediction in the standard Λ cold dark matter model. The curvature can be positive according to the enhanced lensing amplitude in the cosmic microwave background power spectrum confirmed by the Planck Legacy 2018. For the case Λ = 0 the global dynamics of the Hořava-Lifshitz scalar field cosmological model under the background of FLRW was described in [12]-[13], and the case of Λ 6= 0 with zero curvature has also been addressed in [14]. In the present paper we will discuss the global dy- namics of a non-flat universe with Λ 6= 0. Its field equations admit the 5 following dimensionless form dx √ = 6s [−x2 + (u − z)2 + 1] + x [3x2 + 2(u − z)z − 3] , dt dz (1) = z [3x2 + 2(u − z)z − 2] , dt du = u [3x2 + 2(u − z)z] , dt 1 dV (φ) where s = − , κ is a constant, V (φ) is the usual scalar κV (φ) dφ field potential, which admits various mathematical representations (see [9, 11, 16, 17]), it can even be presented in a constant form in scalar cosmology ([18, 19]), and s is assumed to be a constant under the usual exponential potentials. We will provide the detailed information on obtaining cosmological equations (1) in section 2. 2. The cosmological equations In this section we first briefly recall the classic Hořava-Lifshitz gravi- tational theory, where the field content can be represented by the space vector Ni and scalar N , and they are common ‘lapse’ and ‘shift’ vari- ables in general relativity [1], [7], [20]. Then the full metric can be defined as (2) ds2 = −N 2 dt2 + gij (dxi + N i dt)(dxj + N j dt), Ni = gij N j , where gij (i, j = 1, 2, 3) is a spatial metric. The rescaling conversion meets the conditions t → l3 t, xi → lxi , under which gij and N remain unchanged, but N i is scaled to N i → l−2 Ni . 6 F. GAO AND J. LLIBRE Under the detailed-balance condition, the full gravitational action of Hořava-Lifshitz is written as Z  3 √ 2 ij 2  κ2 Sg = dtd x gN K ij K − λK − Cij C ij κ2 2w4 µκ2 ijm l µ2 κ2 (3) + R ∇ √ il j k R − Rij Rij 2w2 g 8   µ2 κ2 1 − 4λ 2 2 − R + ΛR − 3Λ , 8(3λ − 1) 4 where Kij = (ġij −∇i Nj −∇j Ni )/(2N ) denotes the extrinsic curvature,  √ C ij = ijm ∇k 4Rij − Rδij /(4 g) represents the Cotton tensor, and √ ijm / g is the standard general covariant antisymmetric tensor, the indices are to raise and lower with the metric gij . The parameters λ, µ and w are constants (see [1] for more details). For the potential V (φ) we take into account the gravitational action term as follows Z ! 3 √ 3λ − 1 φ̇2 (4) S= dtd x gN − V (φ) , 4 N2 and the metric N i = 0, gij = a2 (t)γij , γij dxi dxj = r2 dΩ22 + dr2 /(1 − kr2 ). Here the function a(t) is the dimensionless rescaling factor of the expanding universe, and γij is a constant curvature metric of maximally symmetric. Without loss of generality we normalize κ2 and N to the 7 number one, and then we can describe the cosmological model as (5)  2 2  2 φ̇2 V (φ) 1 µk 2 2 2µ2 Λk H = + − +µ Λ − , 24 6(3λ − 1) 16(3λ − 1)2 a4 a2  2 2  3 2 φ̇2 V (φ) 1 µk 2 2 2µ2 Λk Ḣ + H = − + + − 3µ Λ − , 2 16 4(3λ − 1) 32(3λ − 1)2 a4 a2 2V 0 (φ) φ̈ + 3H φ̇ + = 0, 3λ − 1 where H = ȧ(t)/a(t) is the Hubble parameter which gives the expansion rate of the universe. According to [7]-[10] we perform the dimensionless transformation (6) p φ̇ V (φ) µ Λµ x= √ , y= √ √ , z= 2 , u= . 2 6H 6H 3λ − 1 4(3λ − 1)a H 4(3λ − 1)H Then it can be followed from equations (5) and (6) that H0 (7) x2 + y 2 − (u − kz)2 = 1, = 2z(z − u) − 3x2 . H Therefore the field equations become the form of the autonomous dy- namical system (1). For more details on system (1) see the equation (5.72) of [8], or equations (287)-(290) of [9], or equations (59)-(61) of [10]. In this paper we study the global dynamics of system (1) in the physical region of interest  G = (x, z, u) : x2 − (u − kz)2 ≤ 1, u, z ∈ R , where k = 1, 0, −1 corresponding to closed, flat, and open universe, respectively. For the case k = 1 we note that f+ (x, z, u) = x2 − 8 F. GAO AND J. LLIBRE (u − z)2 − 1 = 0 is an invariant surface because there is a polynomial √ P = −2 6sx + 6x2 + 4(u − z)z such that ∂f+ dx ∂f+ dz ∂f+ du · + · + · = Pf+ . ∂x dt ∂z dt ∂u dt But f− (x, z, u) = x2 − (u + z)2 − 1 = 0 is not an invariant surface for the case k = −1. In addition it is also noted that system (1) is invariant under the three symmetries (x, z, u) 7→ (x, −z, −u) and (x, z, u) 7→ (−x, −z, −u), (x, z, u) 7→ (−x, z, u) if s = 0, i.e. it is symmetric with respect to the x-axis and additionally with respect to the origin and the plane x = 0 when s = 0. Therefore we divide the study of system (1) into four cases taking into account the existence or not of the symmetric plane and of the invariant surface. Case I: s 6= 0 and k = 1, so system (1) is symmetric with respect to the x-axis and it has the invariant surface f+ (x, z, u) = 0. Case II: s 6= 0 and k = −1, so system (1) is symmetric with respect to the x-axis and it has not the invariant surface. Case III: s = 0 and k = 1, so system (1) is symmetric with respect to the origin and with respect to the x-axis, and it has the invariant surface f+ (x, z, u) = 0. Case IV: s = 0 and k = −1, so system (1) is symmetric with respect to the origin and with respect to the x-axis, and it has not the invariant surface. In section 3.1 we will investigate the phase portraits of case I of sys- tem (1) on the invariant planes and surface, as well as the local phase portraits at the finite and infinite equilibrium points. In section 3.2 we will discuss the phase portraits of case I of system (1) inside the 9 Poincaré ball restricted to the region G. Based on these sections, con- sidering the symmetry of system (1), we will study the global dynamics of system (1) adding its behavior at infinity in section 3.3. In section 4 we will study the case II in the same way as in case I. In sections 5 and 6 we will describe the global dynamics of system (1) in the closed and open universe respectively when the field potential V (φ) of the system takes the form of constant, i.e. s = 0 in cases III and IV. Moreover we will give the final discussion and summary in the last section 7. 3. Case I: s 6= 0, k = 1 3.1. Phase portraits on the invariant planes and surface. In order to analyze in detail the local phase portraits at the finite and infinite equilibrium points and the global phase portrait of system (1) in the region G (refer to [8, 9, 7] or [10] again), we first study its phase portraits on the invariant planes z = 0 and u = 0, as well as on the invariant surface x2 − (u − z)2 = 1, respectively. 3.1.1. The invariant plane z = 0. On this plane system (1) becomes dx √ = 6s (−x2 + u2 + 1) + 3x (x2 − 1) , (8) dt du = 3ux2 . dt There are three equilibrium points e1 = (1, 0), e2 = (−1, 0) and e3 = √ √ (2s/ 6, 0) of system (8), where e1 has eigenvalues 3 and 6 − 2 6s, √ e2 has eigenvalues 3 and 6 + 2 6s, and e3 has eigenvalues 2s2 and 2s2 − 3. Therefore the equilibrium point e1 is a hyperbolic unstable √ node when s < 6/2, and it is a hyperbolic unstable saddle when 10 F. GAO AND J. LLIBRE √ s > 6/2. The equilibrium point e2 is a hyperbolic unstable node √ when s > − 6/2, and it is a hyperbolic unstable saddle when s < √ − 6/2. The equilibrium point e3 is a hyperbolic unstable node when √ √ |s| > 6/2, and it is a hyperbolic unstable saddle when |s| < 6/2. √ Moreover for s = − 6/2, e2 = e3 is a semi-hyperbolic saddle-node, by using the semi-hyperbolic singular point theorem (see Theorem 2.19 √ in [45] for more details). Similarly for s = 6/2, e1 = e3 is also a semi-hyperbolic saddle-node. Since the types and stability of these three finite equilibrium points vary with the different values of s, we summarize them in Table 1. Table 1. Equilibrium points for the different √ values of s, where e1 = (1, 0), e2 = (−1, 0) and e3 = (2s/ 6, 0). Values of s Equilibrium points √ 6 e1 and e3 are unstable nodes, e2 is an unstable saddle −∞ < s < − 2 √ 6 e1 is an unstable node, e2 = e3 is an unstable saddle-node s=− 2 √ √ 6 6 e1 and e2 are unstable nodes, e3 is an unstable saddle − 0, on the positive semi-axis of V , dV /dt = −3U 2 V < 0 indicates that V decreases monotonically, and on the negative semi-axis of V , dV /dt = −3U 2 V > 0 means that V increases monotonically. Near the 12 F. GAO AND J. LLIBRE V V U U e4 e4 (a) (b) Figure 1. The local phase portrait at the equilibrium point eo = (0, 0) of system (11) in (a) when s > 0, and in (b) when s < 0. √ straight line U = 0, dU/dt = 6sV (1 + V 2 ), U increases monotonically when V > 0, and U decreases monotonically when V < 0. Therefore the local phase portrait of the equilibrium point e4 of system (11) is shown in Figure 1(a) when s > 0. Similarly the local phase portrait of e4 is illustrated in Figure 1(b) when s < 0. Therefore the corresponding global phase portrait of system (8) re- stricted to the region x2 − u2 ≤ 1 can be summarized in Figure 2. 3.1.2. The invariant plane u = 0. On this plane system (1) writes dx √ = 6s (−x2 + z 2 + 1) + x (3x2 − 2z 2 − 3) , dt (12) dz = z (3x2 − 2z 2 − 2) . dt There are five equilibrium points e1 = (1, 0), e2 = (−1, 0), e3 = √ √ √ √ √ (2s/ 6, 0), e5 = (2/( 6s), − 1 − s2 /s) and e6 = (2/( 6s), 1 − s2 /s) 13 (a) (b) (c) (d) (e) (f) u x Figure 2. The phase portrait on the invariant plane z = 0 of system (8) restricted to the region x2 − u2 ≤ 1 √ inside the Poincaré √ disc for different √ values of s: (a) s < −√6/2, − 6/2, (c) −√ 6/2 < s < 0, (d) 0 < s < 6/2, (b) s = √ (e) s = 6/2, (f) s > 6/2. Here Eo is the diametrically opposite equilibrium point of eo at infinity. 14 F. GAO AND J. LLIBRE in system (12) when s2 6= 1, and the latter three equilibrium points e3 , e5 and e6 coincide when s2 = 1. The equilibrium point e1 has √ eigenvalues 1 and 6 − 2 6s, it is a hyperbolic unstable node when √ √ s < 6/2, and it is a hyperbolic unstable saddle when s > 6/2. √ The equilibrium point e2 has eigenvalues 1 and 6 + 2 6s, it is a hy- √ perbolic unstable node when s > − 6/2, and it is a hyperbolic un- √ stable saddle when s < − 6/2. The equilibrium point e3 has eigen- values 2s2 − 2 and 2s2 − 3, it is a hyperbolic unstable node when √ |s| < 1 or |s| > 6/2, and it is a hyperbolic unstable saddle when √ 1 < |s| < 6/2. The equilibrium points e5 and e6 have the same eigen- √ √ values −1/2 − 16 − 15s2 /(2|s|) and −1/2 + 16 − 15s2 /(2|s|), they are hyperbolic unstable saddles when |s| < 1, and they are not real singularities when |s| > 1x. √ In addition, for the value s = − 6/2, the equilibrium point e2 = e3 is a semi-hyperbolic saddle-node by using the semi-hyperbolic singular √ point theorem. Similarly for the value s = 6/2, the equilibrium point e1 = e3 is also a semi-hyperbolic saddle-node. For the value s = ±1, the equilibrium point e3 = e5 = e6 is a semi-hyperbolic unstable saddle. Since the types and stability of these five finite equilibrium points change with different values of s, we summarize them in Table 2. By using the Poincaré compactification again on the local chart U1 , then system (12) becomes dU  √  = −U V −V + 6s(U 2 + V 2 − 1) , (13) dt dV  √ √  = −V U 2 ( 6sV − 2) + ( 6sV − 3)(V 2 − 1) . dt 15 Table 2. Equilibrium points for the different √ values of s, where √ e = 1 √ (1, 0), e 2 = (−1, 0), e = 3 √ (2s/√ 6, 0), e5 = (2/( 6s), − 1 − s /s) and e6 = (2/( 6s), 1 − s2 /s). 2 Values of s Equilibrium points √ 6 e1 and e3 are unstable nodes, e2 is an unstable saddle −∞ < s < − 2 √ 6 e1 is an unstable node, e2 = e3 is an unstable saddle-node s=− 2 √ 6 e1 and e2 are unstable nodes, e3 is an unstable saddle − < s < −1 2 e1 and e2 are unstable nodes, s = −1 e3 = e5 = e6 is a semi-hyperbolic unstable saddle e1 and e2 are unstable nodes, e3 is stable node, −1 < s < 1 e5 and e6 are unstable saddles e1 and e2 are unstable nodes, s=1 e3 = e5 = e6 is a semi-hyperbolic unstable saddle √ 6 e1 and e2 are unstable nodes, e3 is an unstable saddle 1 0 means that V increases monotonically, and on the negative semi-axis of V , dV /dt < 0 indicates that V decreases √ monotonically. Moreover dU/dt = 6sV (V 2 + 1) around the straight line U = 0, thus U increases monotonically when sV > 0, and U de- creases monotonically when sV < 0. Therefore the local phase portrait of the semi-hyperbolic equilibrium point e4 in system (15) is illustrated in Figure 3(a) when s > 0. Similarly the local phase portrait of e4 is shown in Figure 3(b) when s < 0. Hence the corresponding global phase portrait of system (12) re- stricted to the region x2 − z 2 ≤ 1 is illustrated in Figures 4 and 5. However we need to pay attention that the global phase portrait of √ system (12) is similar when − 6/2 < s < −1 and s = −1 (See Figures 4(c) and 4(d)). The main difference is that the equilibrium point e3 is √ a hyperbolic unstable saddle when − 6/2 < s < −1, while s = −1 it is a semi-hyperbolic unstable saddle. The same situation occurs when √ s = 1 and 1 < s < 6/2 (See Figures 5(b) and 5(c)). 17 V V U U e4 e4 (a) (b) Figure 3. The local phase portrait at the equilibrium point e4 = (0, 0) of system (15) in (a) when s > 0, and in (b) when s < 0. 3.1.3. The invariant surface f+ (x, z, u) = 0. On this surface system (1) becomes dx √  = x x2 − 1 + 2u x2 − 1 , (16) dt du √  = u x2 + 2 + 2u x2 − 1 , dt when u ≥ z, and for the case u ≤ z system (1) reads dx √  = x x2 − 1 − 2u x2 − 1 , (17) dt du √  = u x2 + 2 − 2u x2 − 1 . dt Note that |x| ≥ 1, then for the case u ≥ z, system (16) admits two equilibrium points e7 = (1, 0) and e8 = (−1, 0). In addition we also note that system (16) is symmetric with respect to u-axis under the symmetry (x, u) 7→ (−x, u), so we only need to discuss the equilibrium point e7 . However the functions on the right side of system (16) have 18 F. GAO AND J. LLIBRE e4 e4 e4 e2 e1 e2 e1 e2 e3 e1 E4 E4 E4 (a) (b) (c) e4 e4 e5 e2 e3 e1 e2 e3 e1 e6 z x E4 E4 (d) (e) Figure 4. The phase portrait on the invariant plane u = 0 restricted to the region x2 − z 2 ≤ 1 inside √ the Poincaré disc √ for different √ values of s: (a) s < − 6/2, (b) s = − 6/2, (c) − 6/2 < s < −1, (d) s = −1, (e) −1 < s < 0. Here E4 is the diametrically opposite equilibrium point of e4 at infinity. 19 e4 e4 e4 e5 e2 e3 e1 e2 e3 e1 e2 e3 e1 e6 E4 E4 E4 (a) (b) (c) e4 e4 e2 e1 e2 e1 z x E4 E4 (d) (e) Figure 5. The phase portrait on the invariant plane u = 0 restricted to the region x2 − z 2 ≤ 1 inside the Poincaré disc for different √ values√of s: (a) 0 < s√< 1, (b) s = 1, (c) 1 < s < 6/2, (d) s = 6/2, (e) s > 6/2. no derivative at e7 , which means that the local dynamics near e7 can- not be studied by the methods in the previous sections 3.1.1 and 3.1.2. 20 F. GAO AND J. LLIBRE When x = 1, the first equation in system (16) is always equal to zero, and the second equation is simplified to du/dt = 3u, i.e. on the in- variant straight line x = 1 the solution is u(t) = ce3t (c is an arbitrary constant), which indicates that u(t) tends to infinity in forward time and leads to e7 in backward time. We know the dynamics on x = 1 near e7 , but we want to know the √ dynamics in a neighborhood of e7 . To this end we set x2 − 1 = y > 0 p (|x| > 1) i.e. x = ± 1 + y 2 (y > 0), considering the aforementioned p symmetry, we only discuss the case x = 1 + y 2 (y > 0) here, then system (9) can be rewritten as follows dy = (y + 2u) (1 + y 2 ) , (18) dt du = u (3 + 2yu + y 2 ) . dt It is obvious that system (18) has a fictitious equilibrium point (0, 0) because y > 0, and its eigenvalues are 1 and 3 respectively, i.e. e9 is a fictitious hyperbolic unstable node. The first equation of system (18) shows that when y + 2u > 0, so y increases monotonically. In contrast if y + 2u < 0, then y decreases monotonically. Note that x and y have the same monotonicity when y > 0, so the local phase portrait of system (16) near e7 and system (18) near (0, 0) have the same local phase portrait. Then considering the symmetry (y, u) 7→ (−y, −u) of system (18), we can find that the local phase portrait of system (16) is shown in Figure 6. Similarly the local phase portrait of system (17) is shown in Figure 7. 21 u -1 0 1 e8 e7 x Figure 6. The phase portrait on the invariant surface f+ (x, z, u) = 0 restricted to the region u ≥ z. u -1 0 1 e8 e7 x Figure 7. The phase portrait on the invariant surface f+ (x, z, u) = 0 restricted to the region u ≤ z. 22 F. GAO AND J. LLIBRE Note that system (16) can be transformed into a polynomial differen- √ tial system by letting y = x2 − 1 when |x| =6 1. In order to study the dynamic behavior of the equilibrium points of system (16) at infinity, we first study the infinite equilibrium points of system (18). On the local chart U1 system (18) writes dU = 2U (1 − U )V 2 , (19) dt dV = −(1 + 2U )V (1 + V 2 ) . dt It is easy to find that all the infinite points of system (19) at V = 0 are equilibrium points. Doing time scale transformation dτ3 = V dt to eliminate the common factor V in system (19), then we have dU = 2U (1 − U )V, dτ3 (20) dV = −(1 + 2U ) (1 + V 2 ) . dτ3 This system has a unique equilibrium point in V = 0, e9 = (−1/2, 0) √ with eigenvalues ± 3, it is a hyperbolic saddle. On the local chart U2 system (18) reads dU = 2(1 − U )V 2 , (21) dt dV = −U (2 + U )V − 3V 3 . dt Using time rescaling dτ4 = V dt we obtain dU = 2(1 − U )V, dτ4 (22) dV = −U (2 + U ) − 3V 2 . dτ4 23 e 10 e9 u y E 10 Figure 8. The global phase portrait of system (18) with y > 0. Here E10 is the diametrically opposite equilibrium point of e10 at infinity. The origin e10 = (0, 0) is an equilibrium point, which is a hyperbolic stable center with eigenvalue ±3i (i is the imaginary unit). Hence e10 is either a weak focus or a center, but since Hf = (2U − U 2 + V 2 − 2/3)/(U − 1)3 is a first integral of system (22) defined at (0, 0), e10 is a center. Then the global phase portrait of system (18) with y > 0 is shown in Figure 8. Therefore, combining Figure 8 with Figures 6 and 7, we obtain the global phase portrait of systems (16) and (17) as shown in Figures 9 and 10, respectively. 3.1.4. The finite equilibrium points. It is easy to find that there are five finite equilibrium points of system (1) because s 6= 0 in case I. The equilibrium point p1 = (1, 0, 0) with eigenvalues 3, 1 and 6 − √ 2 6s, the equilibrium point p2 = (−1, 0, 0) with eigenvalues 3, 1 and √ √ √ 6 + 2 6s, the equilibrium points p3 = ( 6/(3s), − 1 − s2 /s, 0), p4 = 24 F. GAO AND J. LLIBRE e 10 e8 e7 E9 e9 u x E 10 Figure 9. The global phase portrait of system (16) on the invariant surface f+ (x, z, u) = 0 restricted to the region u ≥ z. Here E9 is the symmetric point of e9 with respect to the axis x = 0, and E10 is the diametrically opposite equilibrium point of e10 at infinity. u x Figure 10. The global phase portrait of system (17) on the invariant surface f+ (x, z, u) = 0 restricted to the region u ≤ z. √ √ ( 6/(3s), 1 − s2 /s, 0), p3 and p4 have the same eigenvalues 2, −1/2− √ √ 16 − 15s2 )/(2|s|) and −1/2+ 16 − 15s2 )/(2|s|), and the equilibrium 25 √ point p5 = ( 6s/3, 0, 0) with eigenvalues 2s2 , 2s2 −2, and 2s2 −3. Since different values of s determine the types of these five equilibrium points, we list the relevant results more clearly in Table 3. Table 3. Finite equilibrium points for the different val- ues √ of s, where √ p1 = (1, 0, 0), p√2 = (−1,√0, 0), p3 = ( 6/(3s), − 2 2 √ 1 − s /s, 0), p4 = ( 6/(3s), 1 − s /s, 0) and p5 = ( 6s/3, 0, 0). Values of s Equilibrium points √ 6 p1 and p5 are unstable nodes, −∞ < s < − p2 is a saddle 2 √ 6 p1 is an unstable node, s=− p2 and p5 are non-hyperbolic equilibrium points 2 √ 6 p1 and p2 are unstable nodes, − < s < −1 p5 is a saddle 2 p1 and p2 are unstable nodes, s = −1 p3 , p4 and p5 are non-hyperbolic equilibrium points p1 and p2 are unstable nodes, −1 < s < 1 p3 , p4 and p5 are saddles p1 and p2 are unstable nodes, s=1 p3 , p4 and p5 are non-hyperbolic equilibrium points √ 6 p1 and p2 are unstable nodes, 1 3. More details are shown in region I as well as in regions II and III of Figure 11, respectively. However for the equilibrium points on the hy- perbola 2z1 (z1 −z2 ) = 3 there are six local phase portraits in Figure 12. Note that there is one-dimensional stable manifold in region I and one- dimensional unstable manifold in regions II and III filled with infinite equilibrium points. Since the orbits arriving or ending at equilibrium points at infinity in the different regions cannot collide into finite equi- librium points when they tend to equilibrium points on the hyperbola 27 (a) (b) Figure 11. There is an one-dimensional stable manifold in the region I and one-dimensional unstable manifold in regions II and III on the local chart U1 . Figure 12. There are six kinds of equilibrium points on the hyperbola 2z1 (z1 − z2 ) = 3. 2z1 (z1 − z2 ) = 3 coming from both sides of this hyperbola, there is a hyperbolic sector on the equilibrium points of the branches L1 and L2 of the hyperbola, with the exception of two points p6 and p7 , see the first two pictures in Figure 12 for details. 28 F. GAO AND J. LLIBRE Doing time rescaling dτ5 = z3 dt system (23) becomes dz1 √ = − 6sz1 [(z1 − z2 )2 + z32 − 1] + z1 z3 , dτ5 dz2 √ = − 6sz2 [(z1 − z2 )2 + z32 − 1] + 3z2 z3 , (25) dτ5 dz3 √ = −3 + 2z12 − 2z1 z2 − 6sz3 [(z1 − z2 )2 − 1] dτ5 √ +3z32 − 6sz33 . At z3 = 0 system (25) has two lines |z1 − z2 | = 1 filled with infi- nite equilibrium points and these two lines intersect with hyperbola 2z1 (z1 − z2 ) = 3 at p6 = (−3/2, −1/2, 0) and p7 = (3/2, 1/2, 0), re- spectively. Both of them are hyperbolic with the same eigenvalues  √ √ √ − 3, 3, −2 6s and the corresponding eigenvectors have the fol- lowing form (26) ! ! √ ! √ √ √ √ √ 3 3 3 3 6(8s2 + 1) 6(8s2 + 5) , ,1 , − ,− ,1 , , , 1 . 2 2 2 2 16s 48s √ √ ! √ √ ! √ √ ! 3 3 3 3 6(8s2 + 1) 6(8s2 + 5) − ,− ,1 , , ,1 , − ,− , 1 . 2 2 2 2 16s 48s Then p6 and p7 have an unstable manifold of dimension one (respec- tively two) and a stable manifold of dimension two (respectively one) if s > 0 (respectively s < 0), see the last four pictures in Figure 12 for details. 29 On the local chart U2 , we let x = z1 /z3 , z = 1/z3 , u = z2 /z3 , then system (1) writes dz1  √  = z3 −z1 z3 + 6s(−z12 + (z2 − 1)2 + z32 ) , dt dz2 (27) = 2z2 z32 , dt dz3 = z3 (2 − 3z12 − 2z2 + 2z32 ) . dt By eliminating the common factor z3 in system (27) through time scale transformation dτ6 = z3 dt we get dz1 √ = −z1 z3 + 6s [−z12 + (z2 − 1)2 + z32 ] , dτ6 dz2 (28) = 2z2 z3 , dτ6 dz3 = 2 − 3z12 − 2z2 + 2z32 . dτ6 Note that system (28) with z1 = z3 = 0 has the unique infinite equilib- rium point p8 = (0, 1, 0) with eigenvalues {2i, -2i, 0} and eigenvectors {0, -i, 1}, {0, i, 1}, {1, 0, 0}, so there will be a fold-Hopf bifurcation at the infinite equilibrium point p8 , sometimes called a zero-pair bi- furcation or a Gavrilov-Guckenheimer (see Chapter 5 of [47] for more details). We will not continue to discuss other infinite equilibrium points of this system, because these are already included in the local chart U1 . 30 F. GAO AND J. LLIBRE Similarly on the local chart U3 , we let x = z1 /z3 , z = z2 /z3 , u = 1/z3 , then system (1) becomes dz1  √  = z3 −3z1 z3 + 6s(−z12 + (z2 − 1)2 + z32 ) , dt dz2 (29) = −2z2 z32 , dt dz3 = −3z12 z3 + 2(z2 − 1)z2 z3 . dt In this local chart U3 we only need to study the infinite equilibria located in its origin because all the other infinite equilibrium points have been studied in the local charts U1 and U2 . After changing the time scale dτ7 = z3 dt we obtain dz1 √ = −3z1 z3 + 6s [−z12 + (z2 − 1)2 + z32 ] , dt dz2 (30) = −2z2 z3 , dt dz3 = −3z12 + 2(z2 − 1)z2 . dt Obviously the origin (0, 0, 0) is not an equilibrium point of system (30), so we will not continue to investigate other equilibrium points at infinity in the local chart U3 . In summary the equilibrium points filling up the infinity with these different stable and unstable manifolds are summarized in Figures 11, 12 and 13. 3.2. Phase portrait inside the Poincaré ball restricted to the physical region of interest x2 − (u − z)2 ≤ 1. Note that system (1) is invariant with respect to the x-axis due to the symmetry (x, z, u) 7→ (x, −z, −u). Now we divide the Poincaré ball restricted to the region 31 Figure 13. The sphere S2 (the infinity of R3 ) is filled up with equilibrium points. The stable and unstable man- ifolds of these equilibrium points are different in the regions I, II, III, lines L1 , L2 , and points p6 , p7 , p8 , P6 , P7 , P8 on the sphere defined by the closure of U1 and its symmetric closure of V1 with respect to the origin of R3 . Thus P6 , P7 and P8 are the symmetric points with respect to the origin of p6 , p7 and p8 , respectively. The points N and S denote the north pole and south pole of the Poincaré ball, respec- tively. 32 F. GAO AND J. LLIBRE x2 − (u − z)2 − 1 ≤ 0 into four regions as follows R1 : z ≥ 0, u ≥ 0. R2 : z ≥ 0, u ≤ 0. R3 : z ≤ 0, u ≤ 0. R4 : z ≤ 0, u ≥ 0. Due to the symmetry with respect to the x-axis, we only need to discuss the phase portrait of system (1) in the regions R1 and R2 . Combining the phase portraits on the invariant surface f+ (x, z, u) = 0, on the invariant planes z = 0 and u = 0, and at infinity, we obtain the phase portrait on the boundary of the regions R1 and R2 as shown in Figures 14-17. Now as shown in Figures 18 and 19, we divided the boundary of the regions R1 and R2 into sub-surfaces B11 , B12 , · · · , B15 and B21 , B22 , · · · , B25 respectively, then the phase portrait on the boundary of R1 will be displayed more clearly. Therefore we can find from Figures 13, 15 and 17 that there is a hyperbolic sector at the north pole N on spherical boundary B11 of the Poincaré ball, and N is stable on the back boundary plane B12 . The equilibrium points p1 and p2 are un- stable on the boundary surfaces Bl1 l2 (l1 = 1, 2, l2 = 2, 3, 4). Moreover the equilibrium point p5 is unstable on the back boundary planes B12 and B22 , and it is stable on the bottom boundary planes B13 and B23 and on their intersection. In addition the properties of the remain- ing equilibrium points that are not located at the intersection of these boundary surfaces and planes have been studied in the previous section and will not be repeated here. 33 (a) (b) (c) (d) (e) Figure 14. Phase portrait √ in the boundary √ of the region √ R1 for differ- ent values of s: (a) s < − 6/2, (b) s = − 6/2, (c) − 6/2 < s < −1, (d) s = −1, (e) −1 < s < 0. The surfaces above and below the long dashed line represent regions I and III respectively. 34 F. GAO AND J. LLIBRE (a) (b) (c) (d) (e) Figure 15. Phase portrait in the boundary of the region R√1 for dif- ferent√values of s: (a) √ 0 < s < 1, (b) s = 1, (c) 1 < s < 6/2, (d) s = 6/2, (e) s > 6/2. The surfaces above and below the long dashed line represent regions I and III respectively. 35 (a) (b) (c) (d) (e) Figure 16. Phase portrait √ in the boundary √ of the region √ R2 for differ- ent values of s: (a) s < − 6/2, (b) s = − 6/2, (c) − 6/2 < s < −1, (d) s = −1, (e) −1 < s < 0. The surfaces above and below the long dashed line represent regions I and III respectively. 36 F. GAO AND J. LLIBRE (a) (b) (c) (d) (e) Figure 17. Phase portrait in the boundary of the region R√2 for dif- ferent√values of s: (a) √ 0 < s < 1, (b) s = 1, (c) 1 < s < 6/2, (d) s = 6/2, (e) s > 6/2. The surfaces above and below the long dashed line represent regions I and III respectively. 37 Figure 18. The three boundary surfaces and two bound- ary planes of the region R1 . Figure 19. The three boundary surfaces and two bound- ary planes of the region R2 . 3.3. Dynamics in the interior of the regions R1 and R2 . Without loss of generality and considering the physical region of interest, we take √ s = 6/4, and the dynamics of system (1) can be studied in the same way when we take other values of s. Then the five finite equilibrium points of system (1) have the form p1 = {1, 0, 0}, p2 = {−1, 0, 0}, 38 F. GAO AND J. LLIBRE √ √ p3 = {4/3, − 15/3, 0}, p4 = {4/3, 15/3, 0} and p5 = {1/2, 0, 0}. The dynamical behavior of system (1) inside the region R1 is determined by the behavior of the flow in the following planes and surfaces z = 0, u = 0, f+ (x, z, u) = 0, h1 (x, z, u) = 0, h2 (x, z, u) = 0, h3 (x, z, u) = 0, where √ h1 (x, z, u) = 6s [−x2 + (u − z)2 + 1] + x [3x2 + 2(u − z)z − 3] , h2 (x, z, u) = 3x2 + 2(u − z)z − 2, h3 (x, z, u) = 3x2 + 2(u − z)z. The above planes and surfaces divide the regions R1 and R2 into eleven subregions R1i , i = (1, . . . , 11) and nine subregions R2j , j = (1, . . . , 9) respectively, see Figures 20-23 for more details. The signs of the functions h1 , h2 and h3 in these subregions of R1 and R2 can be found in Tables 4 and 5 respectively. Table 4. Signs of functions h1 , h2 and h3 in the subre- gions of R1 . Functions Positive Negative h1 R11 , R12 , R16 , R17 , R18 , R19 R13 , R14 , R15 , R110 , R111 h2 R11 , R13 , R19 , R110 , R111 R12 , R14 , R15 , R16 , R17 , R18 h3 R11 , R12 , R13 , R14 , R18 , R19 , R110 , R111 R15 , R16 , R17 In addition it should also be noted that in order to avoid visual confusion, we change the dashed lines and solid lines in Figures 21 and 23 to the normal perspective, instead of corresponding to the dashed 39 Figure 20. The surfaces h1 , h2 and h3 restricted to the invariant surface f+ (x, z, u) = 0 and the region R1 of the Poincaré ball. lines and solid lines in Figures 20 and 22 respectively. What we need to pay attention here is that for any value of s, the finite equilibrium point p4 is located at the intersection of the surfaces h1 and h2 on the 40 F. GAO AND J. LLIBRE Figure 21. There are the eleven subregions inside the region R1 of the Poincaré ball. Table 5. Signs of functions h1 , h2 and h3 in the subre- gions of R2 . Functions Positive Negative h1 R21 , R22 , R26 , R27 , R28 R23 , R24 , R25 , R29 h2 R21 , R23 , R28 , R29 R22 , R24 , R25 , R26 , R27 h3 R21 , R22 , R23 , R24 , R27 , R28 , R29 R25 , R26 41 Figure 22. The surfaces h1 , h2 and h3 restricted to the invariant surface f+ (x, z, u) = 0 and the region R2 of the Poincaré ball. invariant plane u = 0, and the finite equilibrium point p5 is located at the intersection of the surface h1 and x-axis (see Figures 21 and 23). For the infinite equilibrium point p7 = (3/2, 1/2, 0) of equations (25), according to the three-dimensional Poincaré transformation in the local 42 F. GAO AND J. LLIBRE Figure 23. There are the nine subregions inside the region R2 of the Poincaré ball. chart U1 , we know that p7 = (1/z3 , 3/(2z3 ), 1/(2z3 )) (z3 → 0) is on the intersection of the surface h3 and the Poincaré sphere. Similarly in the xzu coordinate system we have p6 = (1/z3 , −3/(2z3 ), −1/(2z3 )) (z3 → 0), since the infinite equilibrium points P6 and p6 are symmetric about the origin, then we obtain P6 = (−1/z3 , 3/(2z3 ), 1/(2z3 )) (z3 → 0), which is symmetric with p7 with respect to the plane x = 0. For the infinite equilibrium point e9 = (−1/2, 0) of equations (20) when 43 p u ≥ z, we combine the relationship x = ± 1 + y 2 in equations (18) and the invariant surface f+ (x, z, u) = 0 to know that the coordinate of e9 in the three-dimensional coordinate system xzu can be denoted √ as p9 = ( 1 + V 2 /V, 1/(2V ), −1/(2V )) (V → 0), obviously the infinite equilibrium point p9 is not in the region R1 . But when u ≤ z we √ have p9 = ( 1 + V 2 /V, 3/(2V ), 1/(2V )) (V → 0), and it is exactly located on the surface h1 , on the invariant surface fp = 0 and in their intersection with the Poincaré sphere (see Figure 21). As it is shown in the subregion R11 (see Figure 21) the left side surface is contained in the invariant surface f+ (x, z, u) = 0, the bottom plane is contained in the invariant plane u = 0, the right-back segment surface is contained in h1 = 0, and the right-front segment surface is contained in h2 = 0. From Table 6 we find that the orbits of system (1) increase monotonically along the positive directions of the three coordinate axes, which means that the orbits in R11 come from the finite equilibrium points p1 , p4 , or from the subregion R12 , and then go to the boundary of Poincaré sphere restricted to this subregion. In the subregion R12 , its bottom plane is contained in the invariant plane u = 0, the left side surface is contained in the surface h2 = 0, and the right side surface is contained in the surface h1 = 0. Then we can see that the orbits of system (1) in this subregion monotonically decrease along the positive direction of the z-axis, but increase mono- tonically along the positive direction of the x-axis and u-axis, so the orbits start at the infinite equilibrium points on the Poincaré sphere in this subregion and then enter into the subregion R11 . 44 F. GAO AND J. LLIBRE Table 6. Dynamical behavior in the twenty subregions. Subregions Corresponding Region Increase or decrease R11 h1 > 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 R12 h1 > 0, h2 < 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ > 0 R13 h1 < 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż > 0, u̇ > 0 R14 h1 < 0, h2 < 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ > 0 R15 h1 < 0, h2 < 0, h3 < 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ < 0 R16 h1 > 0, h2 < 0, h3 < 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ < 0 R17 h1 > 0, h2 < 0, h3 < 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ < 0 R18 h1 > 0, h2 < 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ > 0 R19 h1 > 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 R110 h1 < 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż > 0, u̇ > 0 R111 h1 < 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż > 0, u̇ > 0 R21 h1 > 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 R22 h1 > 0, h2 < 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż < 0, u̇ < 0 R23 h1 < 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 R24 h1 < 0, h2 < 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ < 0 R25 h1 < 0, h2 < 0, h3 < 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ > 0 R26 h1 > 0, h2 < 0, h3 < 0, z > 0, u < 0 ẋ > 0, ż < 0, u̇ > 0 R27 h1 > 0, h2 < 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż < 0, u̇ < 0 R28 h1 > 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 R29 h1 < 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 The bottom plane in the subregion R13 is contained in the invariant plane u = 0, the left side surface is contained in the surface h1 = 0, the right side surface is contained in the surface h2 = 0, and the backplane is contained in the invariant plane z = 0, The upper and lower surfaces 45 of the front are composed of the intersection of the surfaces h1 = 0, h2 = 0 and the invariant surface f+ (x, z, u) = 0 on the Poincaré sphere. The orbits in this subregion increase monotonically along the positive z-axis and u-axis, but decrease monotonically along the positive x-axis, so the orbits start at the finite point p1 , and then cross the side surfaces of this subregion and eventually go to the subregion R14 respectively. The bottom plane of the subregion R14 is contained in the invariant plane u = 0, the left surface is contained in the surface composed of the surfaces h1 = 0, h2 = 0 and the invariant plane f+ (x, z, u) = 0, and the right surface is contained in the surface h3 = 0, the back surface is contained in the surface h1 = 0, the backplane is contained in the invariant plane z = 0, the front surface is contained in the surface enclosed by the intersection lines of the surfaces h1 = 0, h2 = 0, h3 = 0 and the invariant surface f+ (x, z, u) = 0 on the Poincaré sphere. The orbits in this subregion increase monotonically along the positive u-axis, and decrease monotonically along the positive two other coordinate axes, which indicates that the orbits of this subregion start from the finite point p4 or the equilibrium points on the Poincaré sphere at infinity in this subregion, or come from the subregion R13 and finally enter into subregions R15 or R18 . The bottom plane of the subregion R15 is contained in the invariant plane u = 0, the front surface is contained in the Poincaré sphere, the left surface is contained in the surface h3 = 0, and the right surface is contained in the surface h1 = 0. The dynamic behavior of the orbits in this subregion are the same as that in the subregion R14 , and they 46 F. GAO AND J. LLIBRE decrease monotonically along the positive direction of the three coordi- nate axes, which means that the orbits of start from subregion R14 or the infinite equilibrium points on the Poincaré sphere in this subregion, and then enter into subregion R16 . The bottom plane of the subregion R16 is contained in the invariant plane u = 0, the left and right surfaces are contained in the surface h3 = 0, the front left surface is contained in the surface h1 = 0 and the front right surface is contained in Poincaré sphere. The orbits in this subregion increase monotonically along the positive x-axis, and decrease monotonically along the positive two other coordinate axes, thereby the orbits originate from the infinite equilibrium points on the Poincaré sphere in this subregion, and then enter into the subregion R18 . The front surface of the subregion R17 is contained in the surface h3 = 0, and the back surface is contained in the invariant plane z = 0. The dynamic behavior of the orbits in this subregion is the same as in the subregion R16 , they come from the subregion R18 , after crossing the right part boundary surface of the subregion R17 and then from its left boundary surface back to the subregion R18 . The bottom plane of the subregion R18 is contained in the invariant plane u = 0, the left surface is contained in the surface h1 = 0, and the upper surface and the right surface are contained in the surface h2 = 0 and the invariant surface f+ = 0, the front surface is contained in the Poincaré sphere, the surface below the front surface is contained in the surface h3 = 0, and the back surface is composed of the invariant plane 47 z = 0 and the surface h3 = 0. The orbits in this subregion decrease monotonically in the positive direction along the z-axis, and decrease monotonically in the positive direction along the other two coordinate axes, for this reason the orbits start at the equilibrium points on the Poincaré sphere at infinity or come from the subregions R14 , R16 , the left side part surface of R17 , and then tend to the subregion R19 . The backplane of the subregion R19 is divided into the left and right parts by the surface h2 = 0, but both of them are contained in the invariant plane z = 0. The left and right side surfaces of R19 are con- tained in the surface composed of the surface h1 = 0 and the invariant surface f+ (x, z, u) = 0, the top surface is contained in the Poincaré sphere, the lower-right surface is contained in the surface h2 = 0, and the lower-left surface is contained in the surface h1 = 0. The orbits in this subregion are monotonically increasing along the positive direction of the three coordinate axes, so the orbits come from the finite equilib- rium points p1 and p2 , or from the subregion R18 , and finally approach the equilibrium points on the Poincaré sphere at infinity. The subregions R110 and R111 are connected by the finite equilibrium point p2 , their front surfaces are contained in the surface h1 = 0, the right surface of R110 and the back surface of R111 are contained in the invariant surface f+ (x, z, u) = 0, the back surface of R110 is contained in the invariant plane z = 0, the top surface of R110 and the surface on the right side of R111 are contained in the Poincaré sphere. The orbits in these two subregions decrease monotonically along the positive direction of the x-axis, and monotonically increase along the positive direction along the other two coordinate axes. Thus the orbits in these 48 F. GAO AND J. LLIBRE two subregions start at the equilibrium point p2 , and then eventually go to the corresponding infinite equilibrium point on the Poincaré sphere. Therefore the dynamic process of the orbits inside the eleven subre- gions of R1 discussed above can be summarized as PS in R12 R12 R11 PS in R11 PS in R110 p1 p4 PS in R14 PS in R111 p2 R13 R14 R15 PS in R15 PS in R19 R19 R18 R16 PS in R16 RPS of R17 LPS of R17 . Note: PS represents Poincaré sphere. RPS denotes the right part of surface. LPS stands for the left part of surface. The flow chart discussed above shows that the orbits of system (1) contained in the region R1 have an ω-limit in the subregions R11 , R19 , R110 and R111 when these four subregions are restricted to the Poincaré sphere at infinity. Furthermore the orbits have an α-limit at the finite equilibrium points p1 , p2 and p4 . In addition the orbits also have an ω-limit in the subregions R12 , R14 , R15 and R16 when they are confined to the Poincaré sphere at infinity. For the subregions in R2 (see Figure 23), the front and bottom sur- faces of the subregion R21 are contained in the Poincaré sphere, the back surface is contained in the invariant plane z = 0, the top sur- face is contained in the invariant plane u = 0, and the left surface is 49 contained in the invariant plane f+ (x, z, u) = 0, the lower right sur- face is contained in the surface h2 = 0, and the upper right surface is contained in the surface h1 = 0. The orbits in this subregion increase monotonically along the positive direction of the x-axis and z-axis, and the monotonically decrease along the positive direction of the u-axis, so the orbits in this subregion start from the finite equilibrium points p1 , p4 and the right surface h2 = 0, and finally approaching the point of infinity at Poincaré sphere. The front and bottom surfaces of the subregion R22 are contained in the Poincaré sphere, the top plane is divided into two parts by the surface h1 = 0, but these two parts are contained in the same invariant plane u = 0, the left side surface is contained in the surface h2 = 0, the right side surface is contained in the surface h3 = 0, and the back surface is contained in the invariant plane z = 0. The orbits in this subregion increase monotonically along the positive direction of the x-axis, and decrease along the positive direction of the other two axes. Therefore the orbits in this subregion start from the connected subregions R24 and R26 , they enter into this subregion through the surfaces h1 = 0 and h3 = 0, and then go through the surface h2 = 0 into the subregion R21 . The top and back planes of the subregion R23 are contained in the invariant plane u = 0 and the invariant plane z = 0 respectively, the left surface is contained in the surface h1 = 0, and the right surface is contained in the surface h2 = 0. The orbits in this subregion increase monotonically along the positive direction of the z-axis and decrease monotonically along the positive direction of the other two axes, which 50 F. GAO AND J. LLIBRE indicates that the orbits in this subregion start at the finite equilibrium point p1 and cross the surface h2 = 0, and then enter into the subregion R24 . The top and back planes of the subregion R24 are also contained in the invariant plane u = 0 and the invariant plane z = 0 respectively, the bottom surface is contained in the surface h1 = 0, and the left sur- face is contained in the surface h2 = 0, the right surface is contained in the surface h3 = 0, and the front surface is contained in the Poincaré sphere. Then the orbits in this subregion are monotonically decreas- ing along the positive direction of the three axes, which means that the orbits in this subregion start at the finite equilibrium point p4 , or come from the subregion R23 and the infinite equilibrium points on the Poincaré sphere, and then directly enter into the subregions R22 and R25 through the surfaces h1 = 0 and h3 = 0 respectively. The top plane of the subregion R25 is contained in the invariant plane u = 0, the bottom surface is contained in the surface h1 = 0, the left surface is contained in the surface h3 = 0, and the front surface is contained in the Poincaré sphere. The orbits in this subregion are monotonically decreasing along the positive direction of the x-axis and z-axis, but monotonically increasing along the positive direction of the u-axis, which indicates that the orbits in this subregion start from the infinity equilibrium point on the Poincaré sphere, then cross the surface h1 = 0 into the subregion R26 . The top plane of the subregion R26 is also contained in the invari- ant plane u = 0, the right and bottom surfaces are contained in the 51 Poincaré sphere, the front surface is contained in the surface h1 = 0, and the back-left and back-right surfaces are contained in the surface h3 = 0. The orbits in this subregion are monotonically decreasing along the positive direction of the z-axis, but are increasing monotonically along the positive directions of the remaining two axes, then the orbits in this subregion start from the subregion R27 and cross the back-right surface h3 = 0 or from the infinite equilibrium points on the Poincaré sphere, and then cross the back-left surface h3 = 0 into the subregion R22 . The subregion R27 consists of the top plane contained in the invariant plane u = 0, the back plane contained in the invariant plane z = 0, the back surface contained in the surface h2 = 0, the front surface contained in the surface h3 = 0, and the bottom surface contained in the Poincaré sphere. The orbits in this subregion are monotonically increasing along the positive x-axis, but are monotonically decreasing along the positive directions of the remaining two axes, so they start from the subregion R28 or the equilibrium points on the Poincaré sphere at infinity, then go through the left surface h3 = 0 to the subregion R26 . The subregion R28 consists of the top plane contained in the invariant plane u = 0, the left plane contained in the invariant plane z = 0, the front surface contained in the surface h2 = 0, the back surface contained in h1 = 0, and the right surface contained in the Poincaré sphere. All the orbits in this subregion are monotonically decreasing along the positive direction of the u-axis, and monotonically increasing along the positive direction of the other two axes, so they start from 52 F. GAO AND J. LLIBRE the finite equilibrium point p2 and then enter into the subregion R27 through the front surface h2 = 0. The top plane of the subregion R29 is contained in the invariant plane u = 0, the left plane is contained in the invariant plane z = 0, the right surface is contained in the Poincaré sphere, and the front surface is contained in the surface h1 = 0, the back surface is contained in the constant surface f+ (x, z, u) = 0. The orbits in this subregion are monotonically increasing along the positive direction of the z-axis, and monotonically decreasing along the positive direction of the other two axes, which means that the orbits start at the finite equilibrium point p2 and finally tend to the equilibrium points on Poincaré sphere at infinity. Therefore the dynamic behavior of the orbits inside the nine subre- gions of R2 discussed above can be represented as R24 R23 PS in R25 R25 p4 PS in R27 R27 R26 R22 R21 p1 . p2 R28 PS in R26 PS in R21 R29 PS in R29 The flow chart in the region R2 discussed above shows that the orbits of system (1) contained in this region have an α-limit at the finite equilibrium points p1 , p2 and p4 , as well as on the Poincaré sphere restricted to the subregions R25 , R26 and R27 . Moreover the orbits 53 have an ω-limit on the Poincaré sphere restricted to the subregions R21 and R29 . Therefore we fully describe all qualitative global dynamic behaviors of system (1). 4. Case II: s 6= 0, k = −1 For an open universe the phase portraits of system (1) on the in- variant planes z = 0 and u = 0 are the same as those in sections 3.1.1 and 3.1.2. Moreover the local phase portraits of the finite and infinite equilibrium points on the Poincaré sphere are consistent with sections 3.1.4 and 3.1.5 respectively. However the physical region of interest f− (x, z, u) = 0 is no more an invariant surface. Taking the same approach as in section 3.2, we also divide the Poincaré ball restricted to the region x2 − (u + z)2 ≤ 1 into four regions as below S1 : z ≥ 0, u ≥ 0. S2 : z ≥ 0, u ≤ 0. S3 : z ≤ 0, u ≤ 0. S4 : z ≤ 0, u ≥ 0. and we shall only study the phase portrait of system (1) in the regions S1 and S2 due to the symmetry as mentioned before. The phase portrait of system (1) on the boundaries of these two regions are the same as the closed universe except for the boundary surfaces B14 , B15 , B24 , B25 contained on the surface f− (x, z, u) = 0 (see Figures 14-19). 4.1. Dynamics in the interior of the regions S1 and S2 . For √ convenience we continue the discussion of the case s = 6/4. The invariant planes z = 0, u = 0, and the surfaces h1 , h2 , h3 as well as f− (x, z, u) = 0 divide the regions S1 and S2 into nine subregions 54 F. GAO AND J. LLIBRE S11 , S12 , · · · , S19 and ten subregions S21 , S22 , · · · , S210 respectively, see Figures 24-26 for more details. The signs of the functions h1 , h2 and h3 in these subregions of S1 and S2 can be found in Tables 7 and 8 respectively. Table 7. Signs of functions h1 , h2 and h3 in the subre- gions of S1 . Functions Positive Negative h1 S11 , S12 , S16 , S17 , S18 S13 , S14 , S15 , S19 h2 S11 , S13 , S19 S12 , S14 , S15 , S16 , S17 , S18 h3 S11 , S12 , S13 , S14 , S18 , S19 S15 , S16 , S17 Table 8. Signs of functions h1 , h2 and h3 in the subre- gions of S2 . Functions Positive Negative h1 S21 , S22 , S26 , S27 , S28 S23 , S24 , S25 , S29 , S210 h2 S21 , S23 , S28 , S29 , S210 S22 , S24 , S25 , S26 , S27 h3 S21 , S22 , S23 , S24 , S27 , S28 , S29 , S210 S25 , S26 As shown in Figure 25 the top surface of the subregion S11 is con- tained in the Poincaré sphere, the bottom surface is contained in the surfaces h1 = 0 and h2 = 0, the bottom plane is contained in the in- variant plane u = 0, and the left and right sides are contained in the surfaces f− = 0 and h1 = 0, the back plane is contained in the invari- ant plane z = 0. Noting the dynamic behavior of the orbits in Table 9, we find that all the orbits in this subregion increase monotonically, 55 (a) (b) Figure 24. The surfaces h1 , h2 and h3 restricted to the surface f− (x, z, u) = 0 and the regions S1 and S2 of the Poincaré ball, respectively. which means that the orbits in this subregion start at the finite equilib- rium points p1 , p2 and p4 , or come from the subregion S12 , eventually approaching the infinite equilibrium points on the Poincaré sphere. 56 F. GAO AND J. LLIBRE Figure 25. The nine subregions inside the region S1 of the Poincaré ball. For the subregion S12 the left and right surfaces are contained in the surfaces h2 = 0 and h1 = 0 respectively, the bottom plane is con- tained in the invariant plane u = 0, and the front plane is contained in the Poincaré sphere. Table 9 shows that the orbits in this subre- gion increase monotonically along the positive direction of the x-axis and u-axis, and decrease monotonically along the positive direction of the z-axis. Therefore the orbits in this subregion start from the infi- nite equilibrium points on the Poincaré sphere, and then cross the left surface h2 = 0 into the subregion S11 . In the subregion S13 the left and right surfaces are contained in the surfaces h1 = 0 and h2 = 0 respectively, the bottom plane is contained 57 Figure 26. The ten subregions inside the region S2 of the Poincaré ball. in the invariant plane u = 0, the back plane is contained in the invariant plane z = 0, and the front plane is contained in the Poincaré sphere. Table 9 means that the orbits in this subregion decrease monotonically along the positive direction of the x-axis, while the orbits increase monotonically along the other two axes. It means that the orbits in this subregion start at the finite equilibrium point p1 and cross the surface h2 = 0 into the neighboring subregion S14 . 58 F. GAO AND J. LLIBRE Table 9. Dynamical behavior in the nineteen subregions. Subregions Corresponding Region Increase or decrease S11 h1 > 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 S12 h1 > 0, h2 < 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ > 0 S13 h1 < 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż > 0, u̇ > 0 S14 h1 < 0, h2 < 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ > 0 S15 h1 < 0, h2 < 0, h3 < 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ < 0 S16 h1 > 0, h2 < 0, h3 < 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ < 0 S17 h1 > 0, h2 < 0, h3 < 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ < 0 S18 h1 > 0, h2 < 0, h3 > 0, z > 0, u > 0 ẋ > 0, ż < 0, u̇ > 0 S19 h1 < 0, h2 > 0, h3 > 0, z > 0, u > 0 ẋ < 0, ż > 0, u̇ > 0 S21 h1 > 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 S22 h1 > 0, h2 < 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż < 0, u̇ < 0 S23 h1 < 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 S24 h1 < 0, h2 < 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ < 0 S25 h1 < 0, h2 < 0, h3 < 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ > 0 S26 h1 > 0, h2 < 0, h3 < 0, z > 0, u < 0 ẋ > 0, ż < 0, u̇ > 0 S27 h1 > 0, h2 < 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż < 0, u̇ < 0 S28 h1 > 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 S29 h1 < 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 S210 h1 < 0, h2 > 0, h3 > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 In the subregion S14 the left side surface is contained in the surfaces h1 = 0 and h2 = 0, the right side surface is contained in the surface h3 = 0, the back plane is contained in the invariant plane z = 0, and the back surface is contained in the surface h1 = 0, the bottom plane is contained in the invariant plane u = 0, and the front surface 59 is contained in the Poincaré sphere. From Table 9 we know that the orbits in this subregion are monotonically decreasing along the positive direction of the x-axis and z-axis, while the orbits are monotonically increasing along the u-axis in the positive direction, then the orbits start at the finite equilibrium point p4 , or the equilibrium points on the Poincaré sphere at infinity in this subregion, and then cross the surfaces h1 = 0 and h3 = 0, and finally enter the adjacent subregions S15 and S18 respectively. The subregion S15 is composed of the left side surface contained in the surface h3 = 0, the right side surface contained in the surface h1 = 0, the bottom surface contained in the invariant plane u = 0, and the front surface contained in the Poincaré sphere. From Table 9 we find that the orbits are monotonically decreasing along the three directions in the positive direction in this subregion, so the orbits originate from the subregion S14 or from the infinite equilibrium points on the Poincaré sphere in the subregion S15 , and then cross the surface h1 = 0 into the connecting subregion S16 . The top surface of the subregion S16 is contained in the surface h3 = 0, the front left surface is contained in the surface h1 = 0, the front right surface is contained in the Poincaré sphere, and the bottom plane is contained in the invariant plane u = 0. The subregion S17 is composed of the front surface contained in h3 = 0 and the back plane contained in the invariant plane z = 0. Table 9 indicates that the orbits have the same dynamic behavior in these two subregions, they increase monotonically along the positive direction of the x-axis and decrease monotonically along the other two axes, so the orbits in the subregion 60 F. GAO AND J. LLIBRE S16 start from the subregion S15 or from the infinite equilibrium points on Poincaré sphere, then cross the surface h3 = 0 into the subregion S18 . The orbits in the subregion S17 come from the subregion S18 and finally return to this subregion, that is the orbits in the subregion S18 can cross the subregion S17 from the right to left. The top surface of the subregion S18 is contained in the surface h2 = 0, the bottom surface is contained in the surface h3 = 0, the bottom plane is contained in the invariant plane u = 0, the left side surface is contained in the surface h1 = 0, and the front surface is contained in the Poincaré sphere, the back surface is contained in the invariant plane z = 0 and the surface h3 = 0. Table 9 implies that the orbits decrease monotonically along the z-axis and increase monotonically along the other two axes in this subregion, so some orbits may start from infinite equilibrium points on Poincaré sphere, some orbits start from the finite equilibrium point p2 , and then cross the surface h2 = 0 into the subregion S18 and eventually into the subregion S11 . The front surface of the subregion S19 is contained in the surface h1 = 0, the bottom plane is contained in the invariant plane u = 0, the back surface part is contained in the invariant plane z = 0 and the surface f− (x, z, u) = 0. Table 9 shows that the orbits in S19 decrease monotonically along the positive direction of the x-axis, and increase monotonically along the positive direction of the other two axes, so the orbits in this subregion start at the finite equilibrium point p2 and then tend to equilibrium points at infinity on the Poincaré sphere or go through the right part of the surface f− (x, z, u) = 0 into outer space. Therefore the dynamic behavior of the orbits of system (1) in the region 61 S1 can be concluded into the following flow chart PS in S11 OS (x < 0) PS in S12 S12 S11 p2 PS in S19 p1 p4 PS in S18 S13 S14 S18 LPS of S17 PS in S14 S15 S16 RPS of S17 PS in S15 PS in S16 . Note: OS represents the outer space. The above flow chart in the region S1 indicates that the orbits of system (1) have an α-limit at the finite equilibrium points p1 , p2 and p4 , as well as on the Poincaré sphere restricted to the subregions S12 , S14 , S15 , S16 and S18 . In addition the orbits have an ω-limit in the subregions S11 and S19 , which are restricted to Poincaré sphere. For the region S2 the left surface of the subregion S21 is contained in the surface f− (x, z, u) = 0, the right surface is contained in the surfaces h1 = 0 and h2 = 0, the top plane is contained in the invariant plane u = 0, the bottom surface and the front triangle-shaped surface are contained in the Poincaré sphere, and the back plane is contained in the invariant plane z = 0. Table 9 shows that the orbits in this subregion decrease monotonically along the positive direction of u-axis, and monotonically increase along the positive directions of the other two axes, so the orbits in this subregion start at finite equilibrium 62 F. GAO AND J. LLIBRE points p1 and p4 , or from the neighboring subregion S22 , they eventually tend to the infinite equilibrium points of the Poincaré sphere that is restricted to the subregion S21 or go through the left part of the surface f− (x, z, u) = 0 into the outer space. The top plane and surface of the subregion S22 are contained in the invariant plane u = 0 and the surface h1 = 0 respectively, the back plane is contained in the invariant plane z = 0, and the left side surface is contained in the surfaces h2 = 0 and f− (x, z, u) = 0, the right side surface is contained in the surface h3 = 0, and the front and bottom surfaces are contained in the Poincaré sphere. Table 9 implies that the orbits in this subregion increase monotonically along the positive x-axis direction and decrease monotonically along the positive z-axis and u-axis directions, which means that these orbits may come from the subregions S24 and S26 , and then enter in the subregion S21 or tend to the infinite equilibrium points on the Poincaré sphere in the subregion S22 , or cross the surface f− (x, z, u) = 0 into the outer space. The structure of subregions S23 and S24 and the dynamic behavior of the orbits in them are the same as those of subregions R23 and R24 (see Figures 26 and 23), then the orbits in the subregion S23 start at the finite equilibrium point p1 and cross the surface h2 = 0, and then enter into the subregion S24 , and the orbits in the subregion S24 start at the finite equilibrium point p4 , the connecting subregion S23 , as well as the infinite equilibrium points on the Poincaré sphere, and then cross the surfaces h1 = 0 and h3 = 0, and eventually tend to the subregions S22 and S25 respectively. 63 For the subregion S25 except that the surface on the lower left side is contained in the surface f− (x, z, u) = 0, the remaining structure is the same as that of R25 . According to the Tables 6 and 9 the monotonicity of the orbits in these two subregions are also the same, so the orbits in the subregion S25 start from the infinity equilibrium points on the Poincaré sphere, or from the outer space through the left side surface f− (x, z, u) = 0, but all eventually cross the surface h1 = 0 into the subregion S26 . The top plane of the subregion S26 is contained in the invariant plane u = 0, the left and right side surfaces are contained in the surfaces h3 = 0 and f− (x, z, u) = 0, and the front surface located at the upper left corner is contained in the surface h1 = 0. The front surface containing the infinite equilibrium point e4 is contained in the Poincaré sphere. Table 9 shows that the orbits decrease monotonically along the positive direction of the z-axis in the subregion S26 , while the orbits along the positive direction of the remaining two axes is just the opposite. Therefore the orbits in this subregion start at the infinite equilibrium points on the Poincaré sphere or from the subregion S27 and the outer space passing through the right side surface f− (x, z, u) = 0, and finally enter the subregion S22 through the left side surface h3 = 0. The front and back surfaces of the subregion S27 are contained in the surfaces h3 = 0 and h2 = 0 respectively, the back plane is contained in the invariant plane z = 0, the top plane is contained in the invariant plane u = 0, and the right side surface is contained in the Poincaré sphere and the right part of the surface f− (x, z, u) = 0. Table 9 shows that the orbits decrease monotonically along the positive direction of 64 F. GAO AND J. LLIBRE the x-axis and increase monotonically along the positive direction of the z-axis and u-axis in this subregion. In this way the orbits in S27 start from the adjacent subregion S28 or come from the outer space passing through the right side surface f− (x, z, u) = 0 and the infinite equilibrium points on the Poincaré sphere, and finally pass through the front surface h3 = 0 and enter the subregion S26 . The structure of the subregion S28 is similar to the subregion S27 discussed above, except that the back surface of the subregion S27 hap- pens to be the front surface of the subregion S28 , and the back surface of the subregion S28 is contained in the surface h1 = 0. Note that Ta- ble 9 implies that the orbits monotonically decrease along the positive direction of the u-axis in this subregion, and increase monotonically along the positive direction of the other two axes. This means that the orbits start at the finite equilibrium point p2 , and then cross the surface h2 = 0 into the subregion S27 . The subregions S29 and S210 are connected by the finite equilibrium point p2 . Their front surfaces are contained in the surface h1 = 0, the bottom surface of S29 and the right surface of S210 are contained in the Poincaré sphere, the right side surface of S29 and the back surface of S210 are contained in the right part of the surface f− (x, z, u) = 0, and the back surface plane of S29 is contained in the invariant plane z = 0. From Table 9 we find that the dynamic behavior of the orbits in these two subregions is the same, both decrease monotonically along the positive direction of the x-axis and u-axis, and increase monotonically along the positive direction of the z-axis. Thus the orbits in these two subregions start at the finite equilibrium point p2 , and then tend to 65 the Poincaré sphere or go through the surface f− (x, z, u) = 0 into the outer space. Accordingly the dynamic behavior of system (1) in the region S2 can be shown by the following flow chart PS in S21 S21 p4 p1 OS (x > 0) S23 PS in S22 S22 S24 PS in S24 S26 S25 PS in S27 S27 OS (x < 0) S28 p2 S29 PS in S29 S210 PS in S210 . This flow chart in the region S2 implies that the orbits of system (1) have an α-limit at the finite equilibrium points p1 , p2 and p4 , and at the infinite equilibrium points on the Poincaré sphere restricted to the subregion S24 . Besides the orbits have an ω-limit in the subregions S21 , S22 , S27 , S29 and S210 , which are restricted to Poincaré sphere. 66 F. GAO AND J. LLIBRE When the field potential V (φ) takes the form of constant, i.e. s = 0. Then system (1) is reduced to dx = x [3x2 + 2(u − z)z − 3] , dt dz (31) = z [3x2 + 2(u − z)z − 2] , dt du = u [3x2 + 2(u − z)z] . dt Now we continue to describe the global dynamics of system (31) in the closed and open universe in sections 5 and 6 respectively. 5. Case III: s = 0, k = 1 5.1. Phase portraits on the invariant planes and surface. In this section we will investigate the local and global phase portraits of the finite and infinite equilibrium points of system (31). The phase portraits on the invariant planes x = 0, z = 0 and u = 0, and on the invariant surface f+ (x, z, u) = 0 are studied in what follows. 5.1.1. The invariant plane x = 0. On this plane system (31) writes dz du (32) = 2z(uz − z 2 − 1), = 2uz(u − z). dt dt It is easy to check that this system has no other finite equilibrium points except the straight line z = 0 which is full of the equilibrium points of system (32). 67 By using the Poincaré compactification z = 1/V and u = U/V , we obtain that system (32) on the local chart U1 has the form dU dV (33) = 2U V 2 , = 2V (1 − U + V 2 ). dt dt Rescaling the time via dτ8 = 2V dt the previous system becomes dU dV (34) = U V, = 1 − U + V 2. dτ8 dτ8 It admits an infinite equilibrium point e11 = (1, 0) with eigenvalues ±i, which indicates that e11 may be either a center or a weak focus of this system. Note that (1 − 2U + V 2 )/U 2 = C is a first integral of system (34), then e11 is a center. On the local chart U2 we have z = U/V and u = 1/V , then system (32) can be rewritten as dU dV (35) = −2U V 2 , = 2(U − 1)U V. dt dt On the local chart U2 we only need to examine the origin (0, 0) of system (35). Obviously e12 is an equilibrium point with eigenvalues zero, the conventional eigenvalue method cannot be used to determine the type of e12 and its local phase portrait. We apply horizontal blow-up by introducing the transformation W = U/V (see [49] for more details), and then we obtain dV dW (36) = 2V 2 W (V W − 1), = −2V W (V − W + V W 2 ). dt dt 68 F. GAO AND J. LLIBRE (a) (b) Figure 27. (a) The local phase portrait at the infinite equilibrium point e12 . (b) The global phase portrait of system (32). E11 and E12 are the diametrally opposite equilibrium points of e11 and e12 in the Poincaré disc, re- spectively. Doing the time transformation dτ9 = 2V W dt we eliminate the common factor, and we have dV dW (37) = V (V W − 1), = −V + W − V W 2 . dτ9 dτ9 Then the equilibrium point e12 = (0, 0) of system (37) with eigenvalues ±1 is a hyperbolic saddle. Eliminating the common factor 2U V of system (35) by taking dτ10 = 2U V dt yields dU/dτ10 = −V and dV /dτ10 = U − 1, which implies that the orbits of the local phase portrait of the infinite equilibrium point e12 decreases monotonically along the V -axis, and increases monotonically along the U -axis when V < 0, or decreases monotonically along the U -axis when V > 0. Therefore the local phase portrait of e12 and the global phase portrait of system (32) are shown in Figure 27. 69 u x Figure 28. The global phase portrait of system (38) on the invariant plane z = 0 restricted to the region x2 − u2 ≤ 1. E15 is the diametrically opposite point of e15 on Poincaré disc. 5.1.2. The invariant plane z = 0. On this plane system (31) becomes dx du (38) = 3x(x2 − 1), = 3ux2 . dt dt This system is exactly the same as system (9) in [13], so the global phase portrait of system (38) is shown in Figure 28. The finite equilibrium points e13 = (1, 0) and e14 = (−1, 0) are hyperbolic unstable nodes. Besides the line x = 0 and the infinity of the local chart U1 are filled with equilibrium points (see [13] for more details). 5.1.3. The invariant plane u = 0. On this plane system (31) is reduced to dx dz (39) = x(3x2 − 2z 2 − 3), = z(3x2 − 2z 2 − 2). dt dt 70 F. GAO AND J. LLIBRE z x Figure 29. The global phase portrait of system (39) on the invariant plane u = 0 restricted to the region x2 − z 2 ≤ 1. E19 and E20 are the diametrically opposite points of e19 and e20 on Poincaré disc respectively. Note that system (39) is the same as system (9) in [13], so the global phase portrait of system (39) is illustrated in Figure 29. The finite equilibrium point e16 = (0, 0) is a hyperbolic stable node, e17 = (1, 0) and e18 = (−1, 0) are hyperbolic unstable nodes. In addition the in- finity of system (39) is full of the equilibrium points (see [13] for more details). 5.1.4. The invariant surface f+ (x, z, u) = 0. On this surface system (31) is rewritten in the same form as systems (16) and (17) in section 3.1.3, so the global phase portrait on the invariant surface f+ (x, z, u) = 0 can be found in Figures 9 and 10 respectively. 5.1.5. The finite equilibrium points. It is easy to find that there are two finite equilibrium points q1 = (1, 0, 0) and q2 = (−1, 0, 0) of system (31), q1 and q2 with the same eigenvalues 6, 3, and 1 are hyperbolic 71 unstable nodes. Besides u-axis is full of equilibrium points of system (31). 5.1.6. Phase portrait on the Poincaré sphere at infinity. According to the three-dimensional Poincaré compactification, system (31) can be rewritten as follows. On the local chart U1 dz1 dz2 dz3 (40) = z1 z32 , = 3z2 z32 , = z3 (2z12 − 2z1 z2 + 3z32 − 3). dt dt dt On the local chart U2 dz1 dz2 dz3 (41) = −z1 z32 , = 2z2 z32 , = z3 (2 − 3z12 − 2z2 + 2z32 ) . dt dt dt On the local chart U3 dz1 dz2 dz3 (42) = −3z1 z32 , = −2z2 z32 , = −3z12 z3 + 2(z2 − 1)z2 z3 . dt dt dt It is noted that equations (40), (41), (42) are the same as systems (23), (27) and (29) when s = 0 in section 3.1.5 respectively, so the Poincaré sphere is filled with equilibrium points of system (31) at infinity. 5.2. Phase portrait inside the Poincaré ball restricted to the physical region of interest x2 − (u − z)2 ≤ 1. As mentioned in sec- tion 2 system (31) is invariant under the three symmetries (x, z, u) 7→ (−x, −z, −u), (x, z, u) 7→ (x, −z, −u) and (x, z, u) 7→ (−x, z, u). Here we divide the Poincaré ball into four regions as follows Q1 : x ≥ 0, z ≥ 0. Q2 : x ≤ 0, z ≥ 0. Q3 : x ≤ 0, z ≤ 0. Q4 : x ≥ 0, z ≤ 0. 72 F. GAO AND J. LLIBRE Figure 30. Phase portrait on the boundary of the region Q1 , and O denotes the center of the Poincaré ball. Due to the above symmetries with respect to the origin, the x-axis, and the invariant plane x = 0, we only need to discuss the phase portrait of system (31) in the region Q1 restricted to the region x2 − (u − z)2 ≤ 1. Joining the phase portraits on the invariant planes x = 0, z = 0 and u = 0, as well as on the invariant surface f+ (x, z, u) = 0, and on the Poincaré sphere at infinity, the phase portrait on the boundary of the region Q1 is displayed in Figure 30. It is noted that all the equilibrium points on the u-axis are stable along the two intersecting boundary planes x = 0 and z = 0, and the finite equilibrium point q1 is unstable on the invariant boundary plane z = 0 and on the invariant boundary surface f+ (x, z, u) = 0. 5.3. Dynamics in the interior of the region Q1 . The dynamics of system (31) inside the region Q1 is governed by the behavior of the 73 orbits in the planes and surfaces as follows x = 0, z = 0, u = 0, f+ (x, z, u) = 0, h0 (x, z, u) = 0, h2 (x, z, u) = 0, h3 (x, z, u) = 0, where h0 (x, z, u) = h2 (x, z, u)−1 = h3 (x, z, u)−3 = 3x2 +2(u−z)z −3. Then the interior of region Q1 is separated into ten subregions Q1n , n = (1, . . . , 10) by the above planes and surfaces, from these there are six subregions above the invariant plane u = 0 and four subregions below. See Figures 31 and 32 for more details. The signs of the functions h0 , h2 and h3 in the subregion of Q1 are given in Table 10. Table 10. Signs of functions h0 , h2 and h3 in the sub- region of Q1 . Functions Positive Negative h0 Q11 , Q13 , Q17 Q12 , Q14 , Q15 , Q16 , Q18 , Q19 , Q110 h2 Q11 , Q12 , Q13 , Q17 , Q18 Q14 , Q15 , Q16 , Q19 , Q110 h3 Q11 , Q12 , Q13 , Q14 , Q17 , Q18 , Q19 Q15 , Q16 , Q110 According to Figure 32 it can be seen that the bottom plane of the subregion Q11 is contained in the invariant plane u = 0, the left surface is contained in the invariant surface f+ (x, z, u) = 0, the right surface is contained in the surface h1 = 0, and the front surface is contained in the Poincaré sphere of the subregion. Table 11 shows that the orbits of system (31) increase monotonically along the positive directions of the three coordinate axes in this subregion, so the orbits in this subregion start at the finite equilibrium point q1 and then finally tend to the equilibrium points at infinity of Poincaré sphere. 74 F. GAO AND J. LLIBRE Figure 31. The surfaces h0 , h2 and h3 restricted to the surface f+ (x, z, u) = 0 and the region Q1 of the Poincaré ball. Similarly the bottom plane, the back-left plane, and the right-back plane of the subregion Q12 are contained in the invariant planes u = 0, z = 0, and x = 0, respectively. The back surface is contained in the sur- face h2 = 0, the top-back surface and the front quadrilateral surface are 75 Figure 32. The ten subregions inside the region Q1 of the Poincaré ball. contained in the Poincaré sphere, the upper and lower surfaces on the left side are contained in the surface h1 = 0 and in the invariant plane f+ = 0, respectively. Then the orbits monotonically decrease along the positive direction of the x-axis in this subregion, and monotonically increase along the positive direction of the other two coordinate axes, that is, the orbits start at the finite equilibrium point q1 and then go to the subregions Q13 and Q14 . 76 F. GAO AND J. LLIBRE Table 11. Dynamical behavior in the ten subregions. Subregions Corresponding Region Increase or decrease Q11 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 Q12 h0 < 0, h2 > 0, h3 > 0, x > 0, z > 0, u > 0 ẋ < 0, ż > 0, u̇ > 0 Q13 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 Q14 h0 < 0, h2 < 0, h3 > 0, x > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ > 0 Q15 h0 < 0, h2 < 0, h3 < 0, x > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ < 0 Q16 h0 < 0, h2 < 0, h3 < 0, x > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ < 0 Q17 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 Q18 h0 < 0, h2 > 0, h3 > 0, x > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 Q19 h0 < 0, h2 < 0, h3 > 0, x > 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ < 0 Q110 h0 < 0, h2 < 0, h3 < 0, x > 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ > 0 The front surface of the subregion Q13 is contained in the invari- ant plane f+ (x, z, u) = 0, the top surface is contained in the Poincaré sphere, the left and right planes of the back are contained in the invari- ant planes z = 0 and x = 0, respectively, and the middle surface of the back is contained in the surface h1 = 0. Table 11 shows that the orbits in this subregion increase monotonically along the positive direction of the three coordinate axes, this means that the orbits in this subregion start from the finite equilibrium point q1 , and may also come from the adjacent subregion Q12 . They tend to the Poincaré sphere at infinity. The left and right planes of the back of the subregion Q14 are con- tained in the invariant planes z = 0 and x = 0 respectively, the back surface is contained in the surface h3 = 0, the top and front surfaces are contained in the Poincaré sphere, and the left surface is contained in 77 the surface h2 = 0 and in the invariant surface f+ (x, z, u) = 0, the right side surface is contained in the surface h3 = 0, and the bottom plane is contained in the invariant plane u = 0. Note Table 11 states that the orbits in this subregion decrease monotonically along the positive direction of the x-axis and z-axis, and increase monotonically along the positive direction of the u-axis, so the orbits in this subregion start from the infinite equilibrium points in the Poincaré sphere or in the subregion Q12 , and then go to the subregion Q16 covering the u-axis, because the entire u-axis is filled with the equilibrium points of system (31). The left side surface of the subregion Q15 is contained in the surface h3 = 0, the back plane and the bottom plane are contained in the invariant planes x = 0 and u = 0, respectively, and the front surface is contained in the Poincaré sphere. The front surface of the subregion Q16 is contained in the surface h3 = 0, and the left and right planes of the back are contained in the invariant planes z = 0 and x = 0, respectively, and the section line of the invariant planes z = 0 and x = 0 in this subregion is the u-axis. The subregions Q15 and Q16 are connected together through the origin of system (31). According to Table 11 the orbits in the two subregions decrease monotonically along the positive directions of the three coordinate axes, indicating that the orbits in the subregion Q15 start in the Poincaré sphere at infinity and after enter in the subregion Q14 or tend to the origin O. The orbits in the subregion Q16 come from the subregion Q14 , and finally go to the finite and infinite equilibrium points of the u-axis. 78 F. GAO AND J. LLIBRE The right and left surfaces of the subregion Q17 are contained in the surface h1 = 0 and in the invariant surface f+ (x, z, u) = 0, the top plane is contained in the invariant plane u = 0, the front and bottom surfaces are contained in the Poincaré sphere, and the back plane is contained in the invariant plane z = 0. Table 11 shows that the orbits in this subregion increase monotonically along the positive direction of the x-axis and z-axis, and decrease monotonically along the positive direction of the u-axis, which indicate that the orbits actually originate from the finite equilibrium point q1 and then run towards the equilibrium points at infinity on the Poincaré sphere. The subregion Q18 has the same composition as in Q17 except that the left and right surfaces are contained in the surfaces h1 = 0 and h2 = 0, respectively. We find from Table 11 that the orbits monoton- ically increase along the positive direction of the z-axis and decrease monotonically along the positive directions of the other two coordi- nate axes. The orbits start at the finite equilibrium point q1 in this subregion and finally enter in the subregion Q19 . The subregion Q19 also has the same composition as Q17 except that the left and right surfaces are contained in the surfaces h2 = 0 and h3 = 0, respectively. Table 11 shows that the orbits monotonically decrease along the positive directions of the three coordinate axes. The orbits begin at the infinite equilibrium points on the Poincaré sphere or come from the subregion Q18 , and then tend to the equilibrium points on the u-axis. 79 In the subregion Q110 the left and right surfaces are contained in the surface h3 = 0 and in the invariant plane x = 0, the top plane is con- tained in the invariant plane u = 0, and the front and bottom surfaces are contained in the Poincaré sphere. It is noted from Table 11 that the orbits in this subregion decrease monotonically along the positive direction of the x-axis and z-axis, and increase monotonically along the positive direction of the u-axis, indicating that the orbits actually start at the equilibrium points on the Poincaré sphere at infinity, and eventually tend to the equilibrium points on the u-axis. In summary the dynamics of the orbits in the ten subregions inside the region Q1 studied above can be sketched in the following flow chart PS in Q13 Q13 PS in Q14 O PS in Q11 Q11 q1 Q12 Q14 Q15 PS in Q15 . PS in Q17 Q17 Q18 Q19 u-axis Q16 PS in Q19 The above flow chart indicates that the orbits of system (31) con- tained in the region Q1 admit an α-limit at the finite equilibrium point q1 . Moreover the orbits also have an α-limit at the equilibrium points in the subregions Q14 , Q15 and Q19 when they are restricted to the Poincaré sphere at infinity. In addition the orbits have an ω-limit at the equilibrium points located on the u-axis and at the infinite equi- librium points restricted to the subregions Q11 , Q13 and Q17 on the Poincaré sphere. 80 F. GAO AND J. LLIBRE 6. Case IV: s = 0, k = −1 In this case system (31) has the same phase portraits on the three invariant planes x = 0, z = 0 and u = 0 as in sections 5.1.1-5.1.3. In addition the local phase portraits of the finite and infinite equilibrium points on the Poincaré sphere are consistent in sections 5.1.4 and 5.1.5 respectively. However the physical region of interest f− (x, z, u) = 0 is no more an invariant surface. Here we again divide the Poincaré ball into four regions restricted to the region x2 − (u + z)2 ≤ 1 as follows T1 : x ≥ 0, z ≥ 0. T2 : x ≤ 0, z ≥ 0. T3 : x ≤ 0, z ≤ 0. T4 : x ≥ 0, z ≤ 0. We only need to examine the phase portrait of system (31) in the region T1 taking into account the symmetries (x, z, u) 7→ (−x, z, u), (x, z, u) 7→ (x, −z, −u) and (x, z, u) 7→ (−x, −z, −u). Then system (31) has the same phase portrait on the boundaries of this region as in section 5.2 except on the non-invariant boundary surface f− (x, z, u) = 0. 6.1. Dynamics in the interior of the region T1 . In order to exam- ine the orbital behavior inside the region T1 , we note that the invariant planes x = 0, z = 0, u = 0, and the surfaces h0 , h2 , h3 as well as f− (x, z, u) = 0 divide the region into ten subregions T11 , T12 , · · · , T110 , see Figures 33-34 for more details. The signs of the functions h0 , h2 and h3 in these subregion of T1 are shown in Table 12. 81 Figure 33. The surfaces h0 , h2 and h3 restricted to the surface f− (x, z, u) = 0 and the region T1 of the Poincaré ball. Table 12. Signs of functions h0 , h2 and h3 in the sub- region of T1 . Functions Positive Negative h0 T11 , T12 , T16 , T17 T13 , T14 , T15 , T18 , T19 , T110 h2 T11 , T12 , T16 , T17 , T18 T13 , T14 , T15 , T19 , T110 h3 T11 , T12 , T13 , T16 , T17 , T18 , T19 T14 , T15 , T110 The left and right surfaces of the subregion T11 in Figure 34 are contained in the surfaces f− (x, z, u) = 0 and h0 = 0, the left and right planes of the back are contained in the invariant planes z = 0 and x = 0, respectively, and the bottom plane is contained in the invariant plane u = 0, the front surface is contained in the Poincaré sphere. Table 13 shows that the orbits in this subregion increase monotonically along 82 F. GAO AND J. LLIBRE Figure 34. The ten subregions inside the region T1 of the Poincaré ball. the positive directions of the three coordinate axes, indicating that the orbits start at the finite equilibrium point q1 or from the adjacent subregion T12 , and then go to the infinite equilibrium points on the Poincaré sphere. The left and right surfaces of the subregion T12 are contained in the surfaces h0 = 0 and h2 = 0, respectively, and the left back plane, right 83 back plane, and bottom plane are contained in the invariant plane z = 0, x = 0, and u = 0 respectively, the top surface and the front surface are contained in the Poincaré sphere. The orbits increase monotonically along the positive direction of the three coordinate axes, indicating that the orbits start at the finite equilibrium point q1 and finally enter the subregion T11 and tend to the infinity equilibrium points on the Poincaré sphere. The composition structure of the subregion T13 is the same as Q14 in Figure 32 except that the left surface in this subregion is completely contained in the surface h2 = 0. In addition the subregions T14 and T15 have the same structure as the subregions Q15 and Q16 in Figure 32, respectively. Moreover the dynamic behavior of the orbits in the subregions T13 , T14 , and T15 is the same as that in the subregions Q14 , Q15 , and Q16 , respectively. That is, the orbits in the subregions T13 and T14 originate from their respective infinite equilibrium points on the Poincaré sphere, and then the orbits in T13 enter the subregion T15 and eventually run to the equilibrium points located on the u-axis, and the orbits in T14 go to the origin O or enter the subregion T13 . In the subregion T16 the front surface is contained in the surface f− (x, z, u) = 0, the right surface is contained in the surface h0 = 0, the back surface is contained in the invariant plane z = 0, and the bottom surface is contained in the Poincaré sphere. Table 13 implies that the orbits increase monotonically along the positive direction of the x-axis and z-axis, and decrease monotonically along the positive direction of the u-axis, so that the orbits in this subregion start from the finite equilibrium point q1 and finally run to the equilibrium points on the 84 F. GAO AND J. LLIBRE Poincaré shpere on the sphere at infinity, or enter in the outer space through the surface f− (x, z, u) = 0. For the subregion T17 the left side surface is contained in the surface f− (x, z, u) = 0, the right side surface is contained in the surface h0 = 0, the top plane is contained in the invariant plane u = 0, and the front surface is contained in the Poincaré sphere. However it can be followed from Table 13 that the dynamic behavior of the orbits in this subregion is consistent with that in the subregion T16 . In the subregion T18 the left side surface and the right side surface are contained in the surfaces h0 = 0 and h2 = 0, respectively, the top plane is contained in the invariant plane u = 0, and the surfaces on the front part and the bottom part are contained in the Poincaré sphere, the middle part of the front surface is contained in the surface f− (x, z, u) = 0. In this subregion the orbits monotonically decrease along the positive direction of the x-axis and z-axis, and increase monotonically along the positive direction of the third axis, which means that the orbits start from the finite equilibrium point q1 and then enter the subregion T19 . The left side surface and the right side surface in the subregion T19 are contained in the surfaces h2 = 0 and h3 = 0, respectively, and the composition structure of the remaining part is the same as the corresponding part of the subregion T18 . The orbits in this subregion monotonically decrease along the positive direction of the three coor- dinate axes, implying that the orbits start at the infinity equilibrium points on the Poincaré sphere or come from the subregion T18 , and finally tend to the equilibrium points on the u-axis. 85 The front and back parts of the left surface of the subregion T110 are contained in the surfaces h3 = 0 and f− (x, z, u) = 0, respectively, the right and top planes are contained in the invariant planes x = 0 and u = 0, respectively, and the front and bottom surfaces are contained in the Poincaré sphere. Note that Table 13 states that the orbits in this subregion decrease monotonically along the positive direction of the x- axis and z-axis, and increase monotonically along the positive direction of the other coordinate axis. Then the orbits in this subregion start at the infinite equilibrium points on the Poincaré sphere and eventually tend to the finite and infinite equilibrium points on the u-axis. Table 13. Dynamical behavior in the ten subregions. Subregions Corresponding Region Increase or decrease T11 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 T12 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u > 0 ẋ > 0, ż > 0, u̇ > 0 T13 h0 < 0, h2 < 0, h3 > 0, x > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ > 0 T14 h0 < 0, h2 < 0, h3 < 0, x > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ < 0 T15 h0 < 0, h2 < 0, h3 < 0, x > 0, z > 0, u > 0 ẋ < 0, ż < 0, u̇ < 0 T16 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 T17 h0 > 0, h2 > 0, h3 > 0, x > 0, z > 0, u < 0 ẋ > 0, ż > 0, u̇ < 0 T18 h0 < 0, h2 > 0, h3 > 0, x > 0, z > 0, u < 0 ẋ < 0, ż > 0, u̇ < 0 T19 h0 < 0, h2 < 0, h3 > 0, x > 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ < 0 T110 h0 < 0, h2 < 0, h3 < 0, x > 0, z > 0, u < 0 ẋ < 0, ż < 0, u̇ > 0 86 F. GAO AND J. LLIBRE In short the dynamic behavior of the orbits of system (31) in the region T1 can be summarized as follows PS in T12 T12 PS in T19 T110 PS in T110 PS in T14 PS in T11 T11 q1 T18 T19 u-axis T15 T13 T14 PS in T16 T16 OS T17 PS in T17 PS in T13 O. The above flow chart in the region T1 indicates that the orbits of system (31) have an α-limit at the finite equilibrium point q1 , and at the infinite equilibrium points located on the Poincaré sphere restricted to the subregions T13 , T14 , T19 and T110 . Furthermore the orbits have an ω- limit at the equilibrium points on the u-axis, as well as the equilibrium points at infinity of the Poincaré sphere in the subregions T11 , T12 , T16 and T17 . 7. Conclusions By taking advantage of the fact that the cosmological equations (1) remain unchanged under the symmetry (x, z, u) 7→ (x, −z, −u)when s 6= 0, and it remains unchanged under the additional two symmetries (x, z, u) 7→ (−x, −z, −u), (x, z, u) 7→ (−x, z, u) when s = 0. Then for a wide range of s in the present paper we completely describe the global phase portrait of Hořava-Lifshitz cosmology in the non-flat universe in the case of non-zero cosmological constant, all of these are located in the physical region of interest G with the invariant boundary surface f+ (x, z, u) = 0 and non-invariant boundary surface f− (x, z, u) = 0, respectively. 87 For the case s = 0 the phase portrait of system (31) restricted to the region G shows that the orbits ultimately move to the finite and infinite equilibrium points located on the u-axis or tend to the infinite equilibrium points on the Poincaré sphere in the subregions Q11 , Q13 and Q17 when the universe is closed. Moreover the orbits of system (31) finally go to the equilibrium points that lie on the u-axis or shift to the infinite equilibrium points restricted to the subregions T11 , T12 , T16 and T17 on the Poincaré sphere when the universe is open. Furthermore we apply the aforementioned symmetries of system (31) to perform simple calculations and find that the unstable equilibrium points q1 and q2 correspond to the universe ruled by dark matter. For the case s 6= 0 in addition to the initial conditions on the invari- ant planes z = 0 and u = 0, the phase portrait shows that the orbital evolution of system (1) in region G eventually tends to the equilibrium points at infinity, which are restricted to the subregions R11 , R19 , R110 and R111 , R21 and R29 on the Poincaré sphere when the universe is closed. For an open universe the orbits of system (1) will go to the infinite equilibrium points restricted to the subregions S11 , S19 , S21 , S22 , S27 , S29 and S210 on the Poincaré sphere. For the studied non-flat universe due to some simple calculation com- bined with the analysis of the phase portrait of system (1) in the pre- vious sections, we conclude from the perspective of cosmology that unstable or non-hyperbolic finite equilibrium points p1 , p2 , p3 , p4 and p5 correspond to the universe dominated by dark matter. We note that there are two finite equilibrium points P17 and P18 , that were considered as non-physical points in [7]-[8], because the value of their 88 F. GAO AND J. LLIBRE corresponding dark matter equation-of-state parameter ωM is 2, but we found that this value is 1/3 and with the notation of our paper corresponds to the equilibrium points p3 and p4 , and of course 1/3 sat- isfies the physical range (0, 1). Furthermore the previous flow charts show that the orbits of the cosmological model in the region G tend to the equilibrium points at infinity, which can be the late-time state of the universe. Moreover the finite equilibrium point p5 can also be the late-time state of the universe when the initial conditions are on the invariant boundary plane u = 0. Besides the orbits will spend a finite lapse of time near the finite equilibrium point p4 on the invariant plane u = 0 before reaching the late-time state p5 or the infinite equilibrium points on the Poincaré sphere. Based on the Hořava-Lifshitz gravity in the non-flat FLRW spacetime with Λ 6= 0, equations (6) shows that the Hubble parameter H tends to 0 in the forward time in this cosmological model. For the late-time state p5 of the universe on the invariant plane u = 0, the equations (6), (7) and H = ȧ(t)/a(t) indicate that the Hubble parameter H is an exponential function, and that the scale factor a(t) of the expanding universe is a double exponential function with respect to the time t, which expands much more quickly than an usual exponential function. Appendix A We recall the results on normally hyperbolic submanifold according to the monograph of Hirsch et al. [48]. Definition 1. It is assumed that φt is a smooth flow on a manifold M and C is a submanifold of M , where C is completely composed of 89 the equilibrium points of the flow. If the tangent bundle of M on C is divided into three invariant subbundles T C, E s , E u under φt satisfying the conditions (A1) dφt contracts E s exponentially, (A2) dφt expands E u exponentially, (A3) T C = tangent bundle of C, then C is called normally hyperbolic submanifold. For normally hyperbolic submanifolds, one usually have smooth sta- ble and unstable manifolds as well as the permanence of these invariant manifolds under small perturbations. To be more precise, we present the following theorem. Theorem 1. If C is a normally hyperbolic submanifold consisting of equilibrium points for a smooth flow φt , then there are smooth stable and unstable manifolds, which tangent along C to E s ⊕T C and E u ⊕T C respectively. In addition, the submanifold C as well as the stable and unstable manifolds are persist under small perturbations of the flow. Acknowledgments The first author gratefully acknowledges the support of the Na- tional Natural Science Foundation of China (NSFC) through grant No. 11672259, and the China Scholarship Council (CSC) through grant No. 201908320086. The second author gratefully acknowledges the support of the Minis- terio de Economı́a, Industria y Competitividad, Agencia Estatal de In- vestigación grants MTM2016-77278-P (FEDER), the Agència de Gestió 90 F. GAO AND J. LLIBRE 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] R. Brandenberger, Matter bounce in Hořava-Lifshitz cosmology, Physical Re- view D 80, 043516, 2009. [4] T. Kobayashi, Y. Urakawa and M. Yamaguchi, Large scale evolution of the curvature perturbation in Hořava-Lifshitz cosmology, Journal of Cosmology and Astroparticle Physics 2009 (11), 015, 2009. [5] S. Carloni, E. Elizalde and P.J. Silva, An analysis of the phase space of Hořava- Lifshitz cosmologies, Classical and Quantum Gravity 27 (4), 045004, 2010. [6] B. Chen, On Hořava-Lifshitz cosmology, Chinese Physics C 35 (5), 429-435, 2011. [7] G. Leon and E.N. Saridakis, Phase-space analysis of Hořava-Lifshitz cosmology, Journal of Cosmology and Astroparticle Physics 2009 (11), 006, 2009. [8] G. Leon and C.R. Fadragas, Cosmological dynamical systems: and their appli- cations, Lambert Academic Publishing, GmbH & Co. KG, Saarbrücken, 2012. [9] G. Leon and A. Paliathanasis, Extended phase-space analysis of the Hořava- Lifshitz cosmology, The European Physical Journal C 79, 746, 2019. [10] A. Paliathanasis and G. Leon, Cosmological solutions in Hořava-Lifshitz scalar field theory, Zeitschrift für Naturforschung A, doi: 10.1515/zna-2020-0003, 2020. [11] C.R. Fadragas, G. Leon and E.N. Saridakis, Dynamical analysis of anisotropic scalar-field cosmologies for a wide range of potentials, Classical and Quantum Gravity 31, 075018, 2014. 91 [12] F.B. Gao and J. Llibre, Global dynamics of the Hořava-Lifshitz cosmological system, General Relativity and Gravitation 51, 152, 2019. [13] F.B. Gao and J. Llibre, Global dynamics of Hořava-Lifshitz cosmology with non-zero curvature and a wide range of potentials, The European Physical Journal C 80, 137, 2020. [14] F.B. Gao and J. Llibre, Global dynamics of the Hořava-Lifshitz cosmology in the presence of non-zero cosmological constant and space curvature zero, preprint, 2019. [15] X. Gao and Y. Wang, R. Brandenberger and A. Riotto, Cosmological pertur- bations in Hořava-Lifshitz gravity, Physical Review D 81, 083508, 2010. [16] D. Escobar, C.R. Fadragas, G. Leon and Y. Leyva, Asymptotic behavior of a scalar field with an arbitrary potential trapped on a Randall-Sundrum’s braneworld: the effect of a negative dark radiation term on a Bianchi I brane, Astrophysics and Space Science 349, 575-602, 2014. [17] A. Alho, J. Hell and C. Uggla, Global dynamics and asymptotics for monomial scalar field potentials and perfect fluids, Classical and Quantum Gravity 32, 145005, 2015. [18] H.C. Kim, Inflation as an attractor in scalar cosmology, Modern Physics Letters A 28 (20), 1350089, 2013. [19] S. Chervon, I. Fomin, V. Yurov and A. Yurov, Scalar field cosmology, World Scientific Publishing Co. Pte. Ltd., 2019. [20] E. Kiritsis and G. Kofinas, Hořava-Lifshitz cosmology, Nuclear Physics B 821, 467-480, 2009. [21] S. Lepe and J. Saavedra, On Hořava-Lifshitz cosmology, Astrophysics and Space Science 350, 839-843, 2014. [22] R. Cordero, H. Garcı́a-Compeán and F.J. Turrubiates, A phase space descrip- tion of the FLRW quantum cosmology in Hořava-Lifshitz type gravity, General Relativity and Gravitation 51, 138, 2019. [23] N.A. Nilsson and E. Czuchry, Hořava-Lifshitz cosmology in light of new data, Physics of the Dark Universe 23, 100253, 2019. 92 F. GAO AND J. LLIBRE [24] E.N. Saridakis, Aspects of Hořava-Lifshitz cosmology, International Journal of Modern Physics D 20 (08), 1485-1504, 2011. [25] 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. [26] E.N. Saridakis, Hořava-Lifshitz dark energy, The European Physical Journal C 67, 229-235, 2010. [27] M. Jamil and E.N. Saridakis, New agegraphic dark energy in Hořava-Lifshitz cosmology, Journal of Cosmology and Astroparticle Physics 2010 (07), 028, 2010. [28] M.R. Setare and M. Jamil, Holographic dark energy with varying gravitational constant in Hořava-Lifshitz cosmology, Journal of Cosmology and Astroparticle Physics 2010 (02), 010, 2010. [29] M.R. Setare, D. Momeni and S.K. Moayedi, Interacting dark energy in Hořava- Lifshitz cosmology, Astrophysics and Space Science 338 (2), 405-410, 2011. [30] A. Sheykhi, S. Ghaffari and H. Moradpour, Ghost dark energy in the de- formed Hořava-Lifshitz cosmology, International Journal of Modern Physics D 28 (06), 1950080, 2019. [31] R.G. Cai, L.M. Cao and N. Ohta, Topological black holes in Hořava-Lifshitz gravity, Physical Review D 80, 024003, 2009. [32] 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. [33] B. Pourhassan, PV criticality of the second order quantum corrected Hořava- Lifshitz black hole, The European Physical Journal C 79, 740, 2019. [34] Z. Nekouee, J. Sadeghi and M. Shokri, The Lagrangian of charged test particle in Hořava-Lifshitz black hole and deformed phase space, International Journal of Modern Physics A 35, 2040056, 2020. [35] K. Lin, N. Yang and J. Li, Electromagnetic quasinormal modes in Hořava- Lifshitz gravity, International Journal of Theoretical Physics 50, 48-55, 2011. 93 [36] A. Borzou, K. Lin and A. Wang, Static electromagnetic fields and charged black holes in general covariant theory of Hořava-Lifshitz gravity, Journal of Cosmology and Astroparticle Physics 2012 (02), 025, 2012. [37] J. Bellorı́n, A. Restuccia and F. Tello-Orti, Anisotropic coupling of gravity and electromagnetism in Hořava-Lifshitz theory, Physical Review D 98, 104018, 2009. [38] A. Restucciaa and F. Tello-Ortizb, Pure electromagnetic-gravitational inter- action in Hořava-Lifshitz theory at the kinetic conformal point, The European Physical Journal C 80, 86, 2020. [39] S. Mukohyama, Hořava-Lifshitz cosmology: a review, Classical and Quantum Gravity 27, 223101, 2010. [40] T.P. Sotiriou, Hořava-Lifshitz gravity: a status report, Journal of Physics: Conference Series 283, 012034, 2011. [41] C. Appignani, R. Casadio and S. Shankaranarayanan, The cosmological con- stant and Hořava-Lifshitz gravity, Journal of Cosmology and Astroparticle Physics 2010 (04), 006, 2010. [42] Ö. Akarsu, J.D. Barrow, L.A. Escamilla and J.A. Vazquez, Graduated dark energy: Observational hints of a spontaneous sign switch in the cosmological constant, Physical Review D 101, 063528, 2020. [43] S. Carlip, Hiding the cosmological constant, Physical Review Letters 123, 131302, 2019. [44] E. Di Valentino, A. Melchiorri and J. Silk, Planck evidence for a closed universe and a possible crisis for cosmology, Nature Astronomy 4, 196-203, 2020. [45] F. Dumortier, J. Llibre and J.C. Artés, Qualitative theory of planar differential systems, Springer-Verlag, Berlin, 2006. [46] A. Cima and J. Llibre, Bounded polynomial vector fields, Transactions of the American Mathematical Society 318(2), 557-579, 1990. [47] Y.A. Kuznetsov, Elements of applied bifurcation theory, Springer-Verlag New York, 2004. 94 F. GAO AND J. LLIBRE [48] M.W. Hirsch, C.C. Pugh and M. Shub, Invariant manifolds, Springer-Verlag Berlin Heidelberg, 1977. [49] 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)