Global dynamics of the Horava-Lifshitz cosmological model in a non-flat universe with non-zero cosmological constant

When the cosmological constant is non-zero the dynamics of the cosmological model based on Horava-Lifshitz gravity in a non-flat universe is characterized by using the qualitative theory of differential equations.


Conclusions 56
Appendix A 57

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 general 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 length renormalization [2], entropy argument [3], cosmology [4]- [30], dark energy [31]- [35], black holes [36]- [40], gravitational waves [41], and electromagnetics [42]- [45]. More information can also be found from the review articles [46]- [48] and the references therein.
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 between its cosmology and observations. Appignani et al. [49] showed that the huge difference between the standard predictions from quantum field theory and the observed value of Λ may have a solution in the Hořava-Lifshitz gravity framework. Akarsu et al. [50] 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 [51] 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, we have a 50:1 odds to conclude that the universe is closed if the Planck CMB data is correct [52]. Besides, Valentino et al. [53] also 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 [14]- [15], and the case of Λ = 0 with zero curvature has also been addressed in [16]. In the present paper we will discuss the global dynamics of a non-flat universe with Λ = 0.
We will provide the detailed information on obtaining cosmological equations in section 2.

The cosmological equations
In this section we first briefly recall the classic Hořava-Lifshitz gravitational theory, where the field content can be represented by the space vector N i and scalar N , and they are common 'lapse' and 'shift' variables in general relativity [1], [9], [22]. Then the full metric can be defined as where g ij (i, j = 1, 2, 3) is a spatial metric. The rescaling conversion meets the conditions t → l 3 t, x i → lx i , under which g ij and N remain unchanged, but N i is scaled to N i → l −2 N i .
Under the detailed-balance condition, the full gravitational action of Hořava-Lifshitz is written as where K ij = (ġ ij − ∇ i N j − ∇ j N i )/(2N ) denotes the extrinsic curvature, C ij = ǫ ijm ∇ k 4R j i − Rδ j i /(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 g ij . 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 and the metric N i = 0, g ij = a 2 (t)γ ij , γ ij dx i dx j = r 2 dΩ 2 2 + dr 2 /(1 − kr 2 ). 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 number one, and then we can describe the cosmological model as (4) where H =ȧ(t)/a(t) is the Hubble parameter which gives the expansion rate of the universe.
In this paper we study the global dynamics of system (7) in the physical region of interest G = (x, z, u) : 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) = x 2 − (u − z) 2 − 1 = 0 is an invariant surface because there is a polynomial The invariant surface here is essential for understanding the complex dynamic behavior of the model (7), because if a point on an orbit of system (7) is located on the invariant surface, then the whole orbit is contained in the surface. But f − (x, z, u) = x 2 − (u + z) 2 − 1 = 0 is not an invariant surface for the case k = −1. In addition it is also noted that system (7) is invariant under the three symmetries (x, z, u) → (x, −z, −u) and (x, z, u) → (−x, −z, −u), (x, z, u) → (−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 (7) into four cases taking into account the existence or not of the symmetric plane and of the invariant surface.
Case I: s = 0 and k = 1, so system (7) is symmetric with respect to the x-axis and it has the invariant surface f + (x, z, u) = 0.
Case II: s = 0 and k = −1, so system (7) is symmetric with respect to the x-axis and it has not the invariant surface.
Case III: s = 0 and k = 1, so system (7) 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 (7) 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 system (7) 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 (7) inside the Poincaré ball restricted to the region G. Based on these sections, considering the symmetry of system (7), we will study the global dynamics of system (7) 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 (7) 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 = 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 (7) in the region G (refer to [10,11,9] or [12] again), we first study its phase portraits on the invariant planes z = 0 and u = 0, as well as on the invariant surface x 2 − (u − z) 2 = 1, respectively.
Therefore the equilibrium point e 1 is a hyperbolic unstable node when s <  [54] for more details). Similarly for s = √ 6/2, e 1 = e 3 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. to draw the vector field of system (8) in the local charts U 1 and U 2 , and then we can determine how the orbits come from or go to infinity. On the local chart U 1 we let x = 1/V and u = U/V , then system (8) becomes Then all the points of system (9) at infinity V = 0 are equilibrium points, we rescale the time dτ 1 = V dt, so this system writes However this system has no equilibrium points at infinity V = 0.
Similarly on the local chart U 2 we have x = U/V and u = 1/V , then system (8) reduces to (11) On the local chart U 2 we only need to study the origin e o = (0, 0) of system (11). Obviously e o is an equilibrium point. Since its linear part is identically zero, we cannot use the usual eigenvalue method to determine the type of the equilibrium point and its local phase portrait, but we note that the straight line V = 0 of system (11) is full of equilibrium points, and if the common factor V in system (11) is eliminated, there are no other equilibrium points in V = 0. When s > 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 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 e 4 of system (11) is shown in Figure 1(a) when s > 0. Similarly the local phase portrait of e 4 is illustrated in Figure 1(b) when s < 0.
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 U 1 , then system (12) becomes Hence all the infinite points of system (13) are filled with equilibrium points at V = 0, doing the time rescaling However this system does not admit any equilibrium point on V = 0. Similarly on the local chart U 2 system (12) reduces to (15) The origin e 4 = (0, 0) of system (15) is a equilibrium point with eigenvalues 2 and 0, but it is not semihyperbolic because it is not isolated in the set of all equilibrium points. It is noted that the axis V = 0 is full of equilibrium points of system (15). For the positive semi-axis of V near e 4 , dV /dt > 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 decreases monotonically when sV < 0. Therefore the local phase portrait of the semihyperbolic equilibrium point e 4 in system (15) is illustrated in Figure 3(a) when s > 0. Similarly the local phase portrait of e 4 is shown in Figure 3(b) when s < 0.
Hence the corresponding global phase portrait of system (12) when u ≥ z, and for the case u ≤ z system (7) reads Note that |x| ≥ 1, then for the case u ≥ z, system (16) admits two equilibrium points e 7 = (1, 0) and e 8 = (−1, 0). In addition we also note that system (16) is symmetric with respect to u-axis under the symmetry (x, u) → (−x, u), so we only need to discuss the equilibrium point e 7 . However the functions on the right side of system (16) have no derivative at e 7 , which means that the local dynamics near e 7 cannot be studied by the methods in the previous sections 3.1.1 and 3.1.2. 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 invariant straight line x = 1 the solution is u(t) = ce 3t (c is an arbitrary constant), which indicates that u(t) tends to infinity in forward time and leads to e 7 in backward time.
We know the dynamics on x = 1 near e 7 , but we want to know the dynamics in a neighborhood of e 7 .
To this end we set √ x 2 − 1 = y > 0 (|x| > 1) i.e. x = ± 1 + y 2 (y > 0), considering the aforementioned symmetry, we only discuss the case x = 1 + y 2 (y > 0) here, then system (9) can be rewritten as follows (18) dy dt = (y + 2u) 1 + y 2 , 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. e 9 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 e 7 and system (18) near (0, 0) have the same local phase portrait. Then considering the symmetry (y, u) → (−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.
Note that system (16) can be transformed into a polynomial differential system by letting y = √ x 2 − 1 when |x| = 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 U 1 system (18) writes  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 (20) dU On the local chart U 2 system (18) reads Using time rescaling dτ 4 = V dt we obtain The origin e 10 = (0, 0) is an equilibrium point, which is a hyperbolic stable center with eigenvalue ±3i (i is the imaginary unit). Hence e 10 is either a weak focus or a center, but since is a first integral of system (22) defined at (0, 0), e 10 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 (  Since different values of s determine the types of these five equilibrium points, we list the relevant results more clearly in Table 3. 3.1.5. Phase portrait on the Poincaré sphere at infinity. According to the three-dimensional Poincaré compactification (see [55] for more details), we set x = 1/z 3 , z = z 1 /z 3 , u = z 2 /z 3 , then on the local chart U 1 Table 3. Finite equilibrium points for the different values of s, where p 1 = (1, 0, 0), p 1 and p 5 are non-hyperbolic equilibrium points, p 2 is an unstable node √ 6 2 < s < +∞ p 1 is a saddle, p 2 and p 5 are unstable nodes system (7) is rewritten as

Values of s
Since the infinity at the different local charts of Poincaré sphere corresponds to z 3 = 0, then for all z 1 , z 2 ∈ R system (23) has the equilibrium point (z 1 , z 2 , 0) with eigenvalues {0, 0, 2z 1 (z 1 − z 2 ) − 3}, which means that the local chart U 1 are full of equilibrium points at infinity. Note that the corresponding eigenvectors are By using normally hyperbolic submanifold theorem (see Appendix A for details), the equilibrium point (z 1 , z 2 , 0) has a one-dimensional stable manifold when 2z 1 (z 1 − z 2 ) < 3 and unstable when 2z More details are shown in region I as well as in regions II and III of Figure 11, respectively. However for the 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 U 1 . equilibrium points on the hyperbola 2z 1 (z 1 −z 2 ) = 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 equilibrium points when they tend to equilibrium points on the hyperbola 2z 1 (z 1 − z 2 ) = 3 coming from both sides of this hyperbola, there is a hyperbolic sector on the equilibrium points of the branches L 1 and L 2 of the hyperbola, with the exception of two points p 6 and p 7 , see the first two pictures in Figure 12 for details.
On the local chart U 2 , we let x = z 1 /z 3 , z = 1/z 3 , u = z 2 /z 3 , then system (7) writes By eliminating the common factor z 3 in system (27) through time scale transformation dτ 6 = z 3 dt we get Note that system (28) with z 1 = z 3 = 0 has the unique infinite equilibrium point p 8 = (0, 1, 0) with eigenvalues [56] 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 U 1 .
Similarly on the local chart U 3 , we let x = z 1 /z 3 , z = z 2 /z 3 , u = 1/z 3 , then system (7) becomes In this local chart U 3 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 U 1 and U 2 . After changing the time scale 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 U 3 .
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 x 2 − (u − z) 2 ≤ 1. Note that system (7) is invariant with respect to the x-axis due to the symmetry (x, z, u) → (x, −z, −u). Now we divide the Poincaré ball restricted to the region Due to the symmetry with respect to the x-axis, we only need to discuss the phase portrait of system (7) in the regions R 1 and R 2 .
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 R 1 and R 2 as shown in Moreover the equilibrium point p 5 is unstable on the back boundary planes B 12 and B 22 , and it is stable Figure 13. The sphere S 2 (the infinity of R 3 ) is filled up with equilibrium points. The stable and unstable manifolds of these equilibrium points are different in the regions I, II, III, lines L 1 , L 2 , and points p 6 , p 7 , p 8 , P 6 , P 7 , P 8 on the sphere defined by the closure of U 1 and its symmetric closure of V 1 with respect to the origin of R 3 . Thus P 6 , P 7 and P 8 are the symmetric points with respect to the origin of p 6 , p 7 and p 8 , respectively. The points N and S denote the north pole and south pole of the Poincaré ball, respectively.
on the bottom boundary planes B 13 and B 23 and on their intersection. In addition the properties of the remaining 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.
3.3. Dynamics in the interior of the regions R 1 and R 2 . Without loss of generality and considering the physical region of interest, we take s = √ 6/4, and the dynamics of system (7) can be studied in the same way when we take other values of s. Then the five finite equilibrium points of system (7) The surfaces above and below the long dashed line represent regions I and III respectively.
The surfaces above and below the long dashed line represent regions I and III respectively. dynamical behavior of system (7) inside the region R 1 is determined by the behavior of the flow in the following planes and surfaces where     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 lines and is on the intersection of the surface h 3 and the Poincaré sphere. Similarly in the xzu coordinate system we have , since the infinite equilibrium points P 6 and p 6 are symmetric about the origin, then we obtain P 6 = (−1/z 3 , 3/(2z 3 ), 1/(2z 3 )) (z 3 → 0), which is symmetric with p 7 with respect to the plane x = 0. For the infinite equilibrium point e 9 = (−1/2, 0) of equations (20) when 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 e 9 in the three-dimensional coordinate system xzu can be denoted as p 9 = ( , obviously the infinite equilibrium point p 9 is not in the region R 1 . But when u ≤ z we have p 9 = ( , and it is exactly located on the surface h 1 , on the invariant surface f p = 0 and in their intersection with the Poincaré sphere (see Figure 21).
As it is shown in the subregion R 11 (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 h 1 = 0, and the right-front segment surface is contained in h 2 = 0. From Table 6 we find that the orbits of system (7) increase monotonically along the positive directions of the three coordinate axes, which means that the orbits in R 11 come from the finite equilibrium points p 1 , p 4 , or from the subregion R 12 , and then go to the boundary of Poincaré sphere restricted to this subregion. In the subregion R 12 , its bottom plane is contained in the invariant plane u = 0, the left side surface is contained in the surface h 2 = 0, and the right side surface is contained in the surface h 1 = 0. Then we can see that the orbits of system (7) The Therefore the dynamic process of the orbits inside the eleven subregions of R 1 discussed above can be summarized as

RPS of R 17
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 (7) contained in the region R 1 have an ω-limit in the subregions R 11 , R 19 , R 110 and R 111 when these four subregions are restricted to the Poincaré sphere at infinity. Furthermore the orbits have an α-limit at the finite equilibrium points p 1 , p 2 and p 4 . In addition the orbits also have an ω-limit in the subregions R 12 , R 14 , R 15 and R 16 when they are confined to the Poincaré sphere at infinity. Therefore the dynamic behavior of the orbits inside the nine subregions of R 2 discussed above can be represented as

PS in R 29
The flow chart in the region R 2 discussed above shows that the orbits of system (7) contained in this region have an α-limit at the finite equilibrium points p 1 , p 2 and p 4 , as well as on the Poincaré sphere restricted to the subregions R 25 , R 26 and R 27 . Moreover the orbits have an ω-limit on the Poincaré sphere restricted to the subregions R 21 and R 29 . Therefore we fully describe all qualitative global dynamic behaviors of system Taking the same approach as in section 3.2, we also divide the Poincaré ball restricted to the region x 2 − (u + z) 2 ≤ 1 into four regions as below and we shall only study the phase portrait of system (7) in the regions S 1 and S 2 due to the symmetry as mentioned before. The phase portrait of system (7) on the boundaries of these two regions are the same as the     Table 9, we find that all the orbits in this subregion increase monotonically, which means that the orbits in this subregion start at the finite equilibrium points p 1 , p 2 and p 4 , or come from the subregion S 12 , eventually approaching the infinite equilibrium points on the Poincaré sphere.  Table 9 shows that the orbits in this subregion 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 infinite equilibrium points on the Poincaré sphere, and then cross the left surface h 2 = 0 into the subregion S 11 .
In the subregion S 13 the left and right surfaces are contained in the surfaces h 1 = 0 and h 2 = 0 respectively, the bottom plane is contained in the invariant plane u = 0, the back plane is contained in the invariant plane  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  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 S 16 start from the subregion S 15 or from the infinite equilibrium points on Poincaré sphere, then cross the surface h 3 = 0 into the subregion S 18 . The orbits in the subregion S 17 come from the subregion S 18 and finally return to this subregion, that is the orbits in the subregion S 18 can cross the subregion S 17 from the right to left. Table 9. Dynamical behavior in the nineteen subregions.

Subregions Corresponding Region Increase or decrease
The top surface of the subregion S 18 is contained in the surface h 2 = 0, the bottom surface is contained in the surface h 3 = 0, the bottom plane is contained in the invariant plane u = 0, the left side surface is contained in the surface h 1 = 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 h 3 = 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 p 2 , and then cross the surface h 2 = 0 into the subregion S 18 and eventually into the subregion S 11 .
The front surface of the subregion S 19 is contained in the surface h 1 = 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 S 19 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 p 2 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 (7) in the region S 1 can be concluded into the following flow chart

PS in S 11
Note: OS represents the outer space.
The above flow chart in the region S 1 indicates that the orbits of system (7) 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 points p 1 and p 4 , or from the neighboring subregion S 22 , they eventually tend to the infinite equilibrium points of the Poincaré sphere that is restricted to the subregion S 21 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 S 22 are contained in the invariant plane u = 0 and the surface h 1 = 0 respectively, the back plane is contained in the invariant plane z = 0, and the left side surface is contained in the surfaces h 2 = 0 and f − (x, z, u) = 0, the right side surface is contained in the surface h 3 = 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 S 24 and S 26 , and then enter in the subregion S 21 or tend to the infinite equilibrium points on the Poincaré sphere in the subregion S 22 , or cross the surface f − (x, z, u) = 0 into the outer space.
The structure of subregions S 23 and S 24 and the dynamic behavior of the orbits in them are the same as those of subregions R 23 and R 24 (see Figures 26 and 23), then the orbits in the subregion S 23 start at the finite equilibrium point p 1 and cross the surface h 2 = 0, and then enter into the subregion S 24 , and the orbits in the subregion S 24 start at the finite equilibrium point p 4 , the connecting subregion S 23 , as well as the infinite equilibrium points on the Poincaré sphere, and then cross the surfaces h 1 = 0 and h 3 = 0, and eventually tend to the subregions S 22 and S 25 respectively.
For the subregion S 25 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 R 25 . According to the Tables 6 and 9 Table 9 shows that the orbits decrease monotonically along the positive direction of the z-axis in the subregion S 26 , 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 S 27 and the outer space passing through the right side surface f − (x, z, u) = 0, and finally enter the subregion S 22 through the left side surface h 3 = 0.
The front and back surfaces of the subregion S 27 are contained in the surfaces h 3 = 0 and h 2 = 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 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 S 27 start from the adjacent subregion S 28 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 h 3 = 0 and enter the subregion S 26 .
The structure of the subregion S 28 is similar to the subregion S 27 discussed above, except that the back surface of the subregion S 27 happens to be the front surface of the subregion S 28 , and the back surface of the subregion S 28 is contained in the surface h 1 = 0. Note that Table 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  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 p 2 , and then tend to the Poincaré sphere or go through the surface f − (x, z, u) = 0 into the outer space.
Accordingly the dynamic behavior of system (7) in the region S 2 can be shown by the following flow chart This flow chart in the region S 2 implies that the orbits of system (7) have an α-limit at the finite equilibrium points p 1 , p 2 and p 4 , and at the infinite equilibrium points on the Poincaré sphere restricted to the subregion S 24 . Besides the orbits have an ω-limit in the subregions S 21 , S 22 , S 27 , S 29 and S 210 , which are restricted to Poincaré sphere.
When the field potential V (φ) takes the form of constant, i.e. s = 0. Then system (7) is reduced to 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).
By using the Poincaré compactification z = 1/V and u = U/V , we obtain that system (32) on the local chart U 1 has the form Rescaling the time via dτ 8 = 2V dt the previous system becomes It admits an infinite equilibrium point e 11 = (1, 0) with eigenvalues ±i, which indicates that e 11 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 e 11 is a center.
On the local chart U 2 we have z = U/V and u = 1/V , then system (32) can be rewritten as On the local chart U 2 we only need to examine the origin (0, 0) of system (35). Obviously e 12 is an equilibrium point with eigenvalues zero, the conventional eigenvalue method cannot be used to determine the type of e 12 and its local phase portrait. We apply horizontal blow-up by introducing the transformation W = U/V (see [58] for more details), and then we obtain Doing the time transformation dτ 9 = 2V W dt we eliminate the common factor, and we have Then the equilibrium point e 12 = (0, 0) of system (37) with eigenvalues ±1 is a hyperbolic saddle.
Eliminating the common factor 2U V of system (35) This system is exactly the same as system (9) in [15], so the global phase portrait of system (38) is shown in Figure 28. The finite equilibrium points e 13 = (1, 0) and e 14 = (−1, 0) are hyperbolic unstable nodes. Besides the line x = 0 and the infinity of the local chart U 1 are filled with equilibrium points (see [15] for more details).
Note that system (39) is the same as system (9) in [15], so the global phase portrait of system (39) is illustrated in Figure 29. The finite equilibrium point e 16 = (0, 0) is a hyperbolic stable node, e 17 = (1, 0) and e 18 = (−1, 0) are hyperbolic unstable nodes. In addition the infinity of system (39) is full of the equilibrium points (see [15] 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) On the local chart U 2 On the local chart U 3 It is noted that equations (40), (41), (42) are the same as systems (23), (27) and (29) 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 Q 1 restricted to the region x 2 − (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 Q 1 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 q 1 is unstable on the invariant boundary plane z = 0 and on the invariant boundary surface f + (x, z, u) = 0.
The signs of the functions h 0 , h 2 and h 3 in the subregion of Q 1 are given in Table 10.  According to Figure 32 it can be seen that the bottom plane of the subregion Q 11 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 h 1 = 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 q 1 and then finally tend to the equilibrium points at infinity of Poincaré sphere. Table 11. Dynamical behavior in the ten subregions.

Subregions Corresponding Region
Increase or decrease Similarly the bottom plane, the back-left plane, and the right-back plane of the subregion Q 12 are contained in the invariant planes u = 0, z = 0, and x = 0, respectively. The back surface is contained in the surface The front surface of the subregion Q 13 is contained in the invariant 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 invariant planes z = 0 and x = 0, respectively, and the middle surface of the back is contained in the surface h 1 = 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 q 1 , and may also come from the adjacent subregion Q 12 . They tend to the Poincaré sphere at infinity.  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 Q 12 , and then go to the subregion Q 16 covering the u-axis, because the entire u-axis is filled with the equilibrium points of system (31).  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 q 1 and then run towards the equilibrium points at infinity on the Poincaré sphere.
The subregion Q 18 has the same composition as in Q 17 except that the left and right surfaces are contained in the surfaces h 1 = 0 and h 2 = 0, respectively. We find from Table 11  The subregion Q 19 also has the same composition as Q 17 except that the left and right surfaces are contained in the surfaces h 2 = 0 and h 3 = 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 Q 18 , and then tend to the equilibrium points on the u-axis.
In the subregion Q 110 the left and right surfaces are contained in the surface h 3 = 0 and in the invariant plane x = 0, the top plane is contained 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 Q 1 studied above can be sketched in the following flow chart

O
The above flow chart indicates that the orbits of system (31) contained in the region Q 1 admit an α-limit at the finite equilibrium point q 1 . Moreover the orbits also have an α-limit at the equilibrium points in the subregions Q 14 , Q 15 and Q 19 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 equilibrium points restricted to the subregions Q 11 , Q 13 and Q 17 on the Poincaré sphere. Here we again divide the Poincaré ball into four regions restricted to the region x 2 − (u + z) 2 ≤ 1 as follows We only need to examine the phase portrait of system (31) in the region T 1 taking into account the symmetries  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 the positive directions of the three coordinate axes, indicating that the orbits start at the finite equilibrium point   In the subregion T 16 the front surface is contained in the surface f − (x, z, u) = 0, the right surface is contained in the surface h 0 = 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 q 1 and finally run to the equilibrium points on the Poincaré shpere on the sphere at infinity, or enter in the outer space through the surface f − (x, z, u) = 0.
For the subregion T 17 the left side surface is contained in the surface f − (x, z, u) = 0, the right side surface is contained in the surface h 0 = 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 T 16 .
In the subregion T 18 the left side surface and the right side surface are contained in the surfaces h 0 = 0 and h 2 = 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 q 1 and then enter the subregion T 19 .
The left side surface and the right side surface in the subregion T 19 are contained in the surfaces h 2 = 0 and h 3 = 0, respectively, and the composition structure of the remaining part is the same as the corresponding part of the subregion T 18 . The orbits in this subregion monotonically decrease along the positive direction of the three coordinate axes, implying that the orbits start at the infinity equilibrium points on the Poincaré sphere or come from the subregion T 18 , and finally tend to the equilibrium points on the u-axis.
The front and back parts of the left surface of the subregion T 110 are contained in the surfaces h 3 = 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.
In short the dynamic behavior of the orbits of system (31) in the region T 1 can be summarized as follows

PS in T 11
The above flow chart in the region T 1 indicates that the orbits of system (31) have an α-limit at the finite equilibrium point q 1 , and at the infinite equilibrium points located on the Poincaré sphere restricted to the subregions T 13 , T 14 , T 19 and T 110 . 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 T 11 , T 12 , T 16 and T 17 . Table 13. Dynamical behavior in the ten subregions.

Subregions Corresponding Region
Increase or decrease For the studied non-flat universe due to some simple calculation combined with the analysis of the phase portrait of system (7) in the previous sections, we conclude from the perspective of cosmology that unstable or non-hyperbolic finite equilibrium points p 1 , p 2 , p 3 , p 4 and p 5 correspond to the universe dominated by dark matter. We note that there are two finite equilibrium points P 17 and P 18 , that were considered as non-physical points in [9]- [10], because the value of their 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  (5), (6) 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. [57]. For normally hyperbolic submanifolds, one usually have smooth stable 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.