This is a preprint of: “Qualitative study of a model with Rastall gravity”, Jaume Llibre, Chara Pantazi, Classical Quantum Gravity, vol. 37(24), 245010:1–25, 2020. DOI: [10.1088/1361-6382/abc188] QUALITATIVE STUDY OF A MODEL WITH RASTALL GRAVITY JAUME LLIBRE1 AND CHARA PANTAZI2 Abstract. We consider the Rastall theory for the flat Friedmann–Robertson–Walker universe filled with a perfect fluid that satisfies a linear equation of state. The corresponding dynamical system is a two dimensional system of polynomial differential equations depending on four parameters. We show that this differential system is always Darboux integrable. In order to study the global dynamics of this family of differential systems we classify all their non– topological equivalent phase portraits in the Poincaré disc and we obtain 16 different dynamical situations for our spacetime. 1. Introduction One characteristic of Eistein’s theory of gravity is the conservation of the energy momentum tensor. So the total rest energy of a system is conserved however there is no experimental evidence for this. Rastall in 1972 proposed a modified theory of gravity where the matter source is described by the energy momentum tensor as in general relativity and also by the metric of the external space. Basically, Rastall suggests that the stress energy tensor of the source of the gravitational field should not conserved and so there is a coupling in a non minimal way between matter and geometry. Rastall theory is considered as an extended theory of gravity, see for example [1] and provide answers to some questions relating to observational cosmology (see additionally [2, 3]) and quantum gravity [4, 5]. There is a kind of similarity between the particle creation process [6, 7] and Rastall theory since both of them do not respect the conservation of the energy momentum tensor. Additionally it has been extensively used in different contexts, see for instance [2, 8, 9, 10, 11, 12, 13]. There is a theory similar to Rastall theory, the so called curvature matter theory of gravity, where the matter and geometry are also coupled to each other in a non minimal way and the standard energy momentum conservation law is also violated, see [14, 15, 16]. There is a discussion about when Rastall theory of gravity coincides with the classical Ein- stein gravity, see for example the work of Visser [17] who supports that the Rastall theory is completely equivalent to the Einstein theory. There is an arguing of Darabi et al. [1] where explain that Rastall [18] did not define a new energy momentum tensor in contrast the claims of Visser [17]. Rastall only assumed that the conservation of the energy momentum tensor is not always valid in a curved spacetime and for this reason suggests a new relation between the energy momentum tensor and the geometry and explain that are coupled to each other in a non minimal way. Hansraj et al. [19] consider a variation of the Rastall parameter in the context of perfect fluids spheres. Visser [17] support that the geometrical part of the field equations are identical in both theories whereas in [19] defends that there is no evidence for this. Hansraj et al. [19] show that the equation of pressure isotropy remains the same in Rastall and in Eistein theory in the sense that the metric potentials for both theories are identical. Additionally using Tolman metrics they proved that there are several differences in the physical behaviour of the corresponding models in both theories, see [19]. In order to have a full knowledge of the dynamical system we should be able to control the orbits that came or escape to infinity. For this reason we use the Poincaré compactification. 1 Thinking about 2–dimensional systems we have that the Poincaré disc D is the closed disc centered at the origin of coordinates of radius one. The interior of D is identified with R2 , and its boundary the circle S1 with the infinity of R2 , in the plane we can go to infinity in as many directions as points has the circle S1 . It is well known that any polynomial vector field X can be extended to an analytic vector field p(X ) defined in the whole closed disc D, called the Poincaré compactification of X . For more details on this compactification see for instance chapter 5 of [20]. Whenever the dynamical system has a first integral then we can reduce the dimension of the system using the expression of that integral. In particular, for a 2–dimensional differential system, the existence of a first integral H̃ enclose a lot of information about the behaviour of the dynamical system. This is due to the fact that the level sets H̃ = c, with c a constant define the orbits of the system. Additionally, using the first integral we can distinguish when an equilibrium point is a center and a focus, otherwise this distinction can be difficult. In general for the differential systems depending on parameters it is extremely hard to decide about the existence or non–existence of a first integral. There are many techniques in the literature about the search of first integrals of a differential system, like Lax pairs [21], Lie symmetries [22], Noether symmetries [23], Darboux theory of integrability [24], Painlevé analysis [25], Differential Galois Theory [26], ... This variety of techniques illustrate the difficulty of finding first integrals. There are also works in the literature that consider generalized Rastall gravity with variable Rastall parameter, see for example [12, 27]. The aims of this paper are the following ones: • We present the equations of our model applying Rastall theory for the flat Friedmann– Robertson–Walker universe filled with a perfect fluid that satisfies a linear equation of state following the work of Babichev and collaborators [28]. Note that the correspond- ing dynamical system is a two dimensional system of polynomial differential equations depending on four parameters, see section 2. We apply affine transformations in the variables and we achieve to reduce our system into a two parametric family maintaining the polynomial structure of our system. • Since our dynamical system is polynomial we apply the Darboux theory of integrability in order to decide about the integrability of the system. Our result is given in Theorem 4 where we prove that for all the values of the parameters the dynamical system has a first integral that can be constructed easily using the invariant curves of the differential system. This is an important information because we have explicitly the expressions of the first integrals for all the values of the parameters. Therefore all the trajectories of the system are contained in the level sets of such first integrals and additionally helps us to decide about some singular points (to distinguish between center and focus). • Using the Poincaré compactification, which allows to study the solutions of the differ- ential system near the infinity, we are able to present the global dynamics of our two dimensional family when the parameters vary. Note that there is no numerical approxi- mation in the parameters and the study is purely analytic. Additionally we present the topological classification of the global phase portraits in the Poincaré disk and we find 16 dynamical different situations in the spacetime, see Theorem 5 and Figure 2. • The physical and cosmological consequences of our study are presented in Theorem 7. • We also study the existence of singular/no singular cosmologies when the Rastall pa- rameter λ is not a constant but it depends on the time. For the radiation era and for a choose of λ we find a singular universe, see section 3. For a pressureless matter we prove that Rastall parameter does not affect the evolution of the universe and we can obtain a non singular cosmology. 2 2. The dynamical system and statement of the main results We recall that Rastall theory [18, 29] suggest that the stress energy tensor of the source of the gravitational field should not conserved, namely it should be T;νµν 6= 0. Rastall proposed to take T;νµν = λR,µ , where R is the Ricci scalar and λ is the Rastall constant parameter. Then the gravitational field equations can be written as (see also [18]) Gµν + kλgµν R = kTµν , (1) The dynamical systems are used for studying the global properties of several cosmologies not only in general relativity [30] but also in other higher order theories of gravity [31, 32, 33]. Using such techniques some qualitative properties are obtained related to the singular points of the system and their stability, but also the invariant surfaces or invariant curves. The importance of such techniques is that there is no need of solving the system explicitly in order to know the qualitative behavior of the solutions, i.e. where they born or they die, if they are periodic or not, ..., see for example [34, 35]. When the system depends on parameters then a careful classification of all the features of the system is needed according to the parameters. where Gµν is the Einstein tensor and k is the gravitational constant in Rastall theory. Rastall showed that (1) yields to R(4kλ − 1) = T , where T is the trace of the energy momentum tensor. Since T is not always zero and λ is a constant, Rastall [18] indicates that the expression 4kλ − 1 cannot be zero. Moradpour et al. [36] using the concept of the Newtonian limit showed that also 6kλ − 1 6= 0. The case λ = 0 and k = 8π gives rise to the classical Einsten’s field equation of general relativity whenever T;νµν = 0, see [18]. For the Friedmann–Robertson–Walker universe filled with a perfect fluid with Tνµ = diag(ρ, −p, −p, −p), the gravitational equations (1) can be written according to Yuan and Huang [37] like κ 3(1 − 4kλ)H 2 − 6kλḢ + 3(1 − 2kλ) = kρ, a2 κ (2) 3(1 − 4kλ)H 2 + 2(1 − 3kλ)Ḣ + (1 − 6kλ) 2 = −kp, a and H = ȧ/a is the Hubble parameter. Here the dot represents time derivative. The a(t) is the scale factor (dimensionless). Here κ = −1, 0, +1 for open, flat and closed spatial sections respectively. As usual ρ is the energy density and p is the pressure. Here we consider the case of flat spatial sections and so we consider that κ = 0. Hence, from the first equation of (2) and considering that λ 6= 0 and k 6= 0 we obtain 4kλ − 1 2 1 Ḣ = − H − ρ. (3) 2kλ 6λ Now the contracted Bianchi identities Gµν;µ = 0 yields to the equation of continuity (see [38]) (3kλ − 1)ρ̇ + 3kλṗ + 3(4kλ − 1)(ρ + p)H = 0. (4) A perfect fluid with linear equation of state has the following form (see the work of Babichev and collaborators [28]) p = αρ − ρ0 , (5) where α and ρ0 are constant parameters. (Note that equation p = αρ − ρ0 = α(ρ − ρ̃0 ) with ρ0 = αρ̃0 is in agreement with reference [28].) Now equation (4) using relation (5) and for 3kλ(1 + α) − 1 6= 0 becomes 3(1 − 4kλ) ρ̇ = ((1 + α)ρ − ρ0 ) H. (6) 3kλ(1 + α) − 1 3 Note that equations (3) and (6) defines a polynomial quadratic differential system in the vari- ables H and ρ. Set H = x and ρ = y. In these new variables, the system writes 1 − 4kλ 2 1 ẋ = x − y, 2kλ 6λ (7) 3(1 − 4kλ) ẏ = ((1 + α)y − ρ0 ) x, 3kλ(1 + α) − 1 and note that is a four parametric family. Remark 1. For system (7) we assume that λk 6= 0 and 3kλ(1 + α) − 1 6= 0. Note that from a mathematical point of view we can allow the case that 1 − 4kλ = 0 and 1 − 6kλ = 0 unlike the Rastall theory [18] and see also [36]. Also note that for α = 1 we have that 3kλ(1 + α) − 1 = 6kλ − 1. The validity of the case 4kλ − 1 = 0 has been also investigated in [12]. Due to equation (5) we have that ρ + p = (1 + α)ρ − ρ0 . Null energy condition holds when ρ + p ≥ 0, and so holds when (1 + α)ρ − ρ0 ≥ 0, or equivalently (1 + α)y − ρ0 ≥ 0, and it is violated in the phantom region (ρ + p < 0). Hence, for α > −1 the energy density is increasing and for α < −1 is decreasing. Note that the condition ρ + p = 0 represents the straight line g1 = (1 + α)y − ρ0 = 0 and this line is invariant for system (7) with cofactor 3 (1 − 4 kλ) (1 + α) Kg1 = x. 3 kλ (1 + α) − 1 For more details about invariant algebraic curves, cofactors and Darboux first integrals see Appendix 3.1. System (7) for 4 kλ−1 6= 0 has the invariant conic g2 = 3 (4 kλ − 1) x2 −k (3 kλ(1 + α) − 1) y+ 3 k 2 λ ρ0 = 0 with cofactor Kg2 = ((1 − 4 kλ)/(kλ))x. System (7) has the three (whenever are defined) finite singular points: Q0 = (0, 0) and s ! k ρ0 ρ0 Q± = ± , . 3 (1 − 4 kλ) (1 + α) 1 + α The points Q± belong to the straight line g1 = 0 and also to the conic g2 = 0. The point Q0 belongs to the phantom region when ρ0 > 0. For ρ0 = 0 the equation of state (5) becomes p = αρ. In this case the points Q± collide with the origin Q0 and belong to the phantom dividing line g1 = 0. For α 6= −1 system (7) has the Darboux first integral 1−3kλ(1+α) H̃(x, y) = (3 (4 kλ − 1) x2 − k (3 kλ(1 + α) − 1) y + 3 k 2 λ ρ0 )((1 + α)y − ρ0 ) 3kλ(1+α) . For α = −1 the equation of state (5) becomes p = −ρ − ρ0 , and so p + ρ = −ρ0 . In this case system (7) has only one singular point, the origin Q0 and admits the Darboux first integral   2 2 y H̃(x, y) = (3 (4 kλ − 1) x + ky + 3 k λ ρ0 )exp − . 3ρ0 kλ Remark 2. For 4kλ − 1 = 0 system (7) is the linear system ẋ = −y/(6λ), ẏ = 0 and has the polynomial first integral H̃(x, y) = y. Moreover, y(t) = C and x(t) = −Ct/(6λ) + C1 with C, C1 real constants. After the change of time ds = (−y/(6λ))dt system becomes ẋ = 1 and ẏ = 0. The phase portrait of this linear system in the Poincaré disc is given in Figure 1. After Remark 2 we assume that 4kλ − 1 6= 0. Consider the change of coordinates 2kλ 12kλ2 x=− X, y= Y. (8) 4 kλ − 1 4 kλ − 1 4 Figure 1. The phase portrait of system (7) for the case 1 − 4kλ = 0. Note that the system is a linear differential system. The line y = 0 is filled of singular points. According to Rastall this case has no physical meaning, see Remark 1. Then system (7) is rewritten Ẋ = Y + X 2 , (4 kλ − 1) ρ0 kλ (1 + α) (9) Ẏ = − X +6 XY. 2λ (3 kλ (α + 1) − 1) 3 kλ (α + 1) − 1 Setting (4 kλ − 1) ρ0 6kλ (1 + α) a=− , and b = , (10) 2λ (3 kλ (α + 1) − 1) 3 kλ (α + 1) − 1 system (9) becomes Ẋ = Y + X 2 , Ẏ = X(a + bY ). (11) 6 0, λ 6= 0, 3 kλ (α + 1) − 1 6= 0, 4kλ − 1 6= 0, ρ0 and α are real constants. We recall that k = There is an abuse of notation relating to the parameter a. Note that the parameter a that appears in the differential system (11) is a parameter, and it is not the scale factor relating with the Hubble parameter. We rename the variables of system (11) and we can rewrite that system as ẋ = y + x2 , ẏ = x(a + by), (12) with a, b real constants. Remark 3. For 4kλ − 1 6= 0 system (12) is the same as system (7) after the scaling of coordinates (8). However is more easy to work with system (12) because it depends only on two parameters whereas system (7) depends on four parameters. Since system (12) is a polynomial differential system we can apply the Darboux theory of integrability [24] for finding a first integral of this system. This is a constructive method that can provide the explicit expression of a first integral using the invariant algebraic curves of the system, see for more details the Appendix 6.1. Our first result is about the Darboux integrability of system (12). Theorem 4. System (12) is Darboux integrable for all the values of the parameters a and b. More precisely system has the following first integrals. (a) If b 6= 0, 2 system (12) has the first integral ((2 − b)x2 + 2y + a)b H̃(x, y) = . (13) (by + a)2 (b) If b = 2 system (12) has the first integral   1 2 x2 − a H̃(x, y) = exp . (14) 2y + a 2y + a 5 (c) If b = 0 and a 6= 0 system (12) has the first integral a/2 H̃(x, y) = 2 x2 + 2 y + a exp(−y). (15) (d) If a = b = 0 system (12) has the first integral H(x, y) = y. The proof of Theorem 4 is given in section 4. Note that if a 6= 0 then the first integral of system (7) is always well defined in the neighbor- hood of the singular point Q0 . b . A12 . S=13, R=4 . A11 A1 . . . . . . . ... . . . . . S=10, R=2 A13 . A10 S=11, R=2 . A2 S=17, R=4 S=10, R=2 . . . . . . . ... . b=2 . A14 . .. . . A9 S=17, R=4 A3 . . . . . . S=13, R=4 . ... . . S=11, R=3 . B3 B2 B1 . . . . S=19, R=6 . . . b=1 S=10, R=2 . . A15 A8 . . . . .A4 . . . . . . S=13, R=4 . . . . S=17, R=4 . A16 B0 A5 . . . . . . . . . a . . S=13, R=4 S=10, R=2 . . . A6 A17 A7 . ... . . . . . . . . . . S=17, R=5 S=13, R=4 S=13, R=4 Figure 2. The phase portraits of system (12) according with the values of the parameters a and b. In each phase portrait S denotes the number of separatrices and R denotes the number of canonical regions. Our next goal is to study the dynamics of systems (7) depending on its parameters a and b, i.e. we shall describe all the phase portraits of that system in the Poincaré disc. Therefore we also study the dynamics of systems (7) in a neigborhood of the infinity. For this we consider the two cases, when 4kλ − 1 6= 0 and so we study system (12) and the case when 4kλ − 1 = 0 which is a linear system, see Remark 2. Thus our main result is the following. 6 Theorem 5. System (12) has 16 non–topological equivalent phase portraits in the Poincaré disc, see Figure 2. The proof of Theorem 5 is given in section 4. The phase portraits of system (7) are these 16 phase portraits that appear in Theorem 5 and we should add the phase portrait when 4kλ − 1 = 0 which appears in Figure 1, see also Remark 2. Note that this additional phase portrait of Figure 1 is non–topological equivalent to the ones of Figure 2. So we have the following result. Corollary 6. System (7) has 17 non–topological equivalent phase portraits in the Poincaré disc. The physical and cosmological consequences of our dynamical study of the model with Rastall gravity are summarized in the next theorem, which is proved along the section 4. Theorem 7. Let a and b the parameters defined in (10), H the Hubble constant and ρ the energy density. Then the flat de-Sitter model is (a) empty and unstable if a > 0 and (H, ρ) = (0, 0), (b) empty and stable if a < 0 and (H, ρ) =p (0, 0), (c) unstable if a > 0, b > 0 and (H, ρ) = (p a/b, −a/b), (d) stable if a > 0, b > 0 and (H, ρ) = (− a/b, p−a/b), and (e) unstable if if a < 0, b < 0 and (H, ρ) = (± a/b, −a/b). There are some interesting works in the literature relating gravitational collapsing and Rastall theory, see for exemple [39]. Remark 8. Here we consider ρ0 = 0 and from equations (2) we obtain 3(1 + α)(1 − 4kλ) 2Ḣ + bH 2 = 0, b= , (16) 1 − 3kλ(1 + α) and since H = ȧ/a the general solution of this equation is 2 a(t) = C2 (bt + 2C1 ) b with C1 , C2 are arbitrary constants. Using the initial conditions a(t0 ) = a0 and a(ts ) = 0 we have that  2 t − ts b 2 a(t) = a0 , H(t) = , t0 − ts b(t − ts ) and with [t0 , ts ] we denote the time interval of the collapse process, so t ∈ [t0 , ts ]. We need that the scale factor a(t) vanishes at a finite time, so we need to assume that b > 0. 3. The generalized Rastall theory with varying Rastall parameter In this section we investigate the case where the Rastall parameter λ is not a constant. So we choose λ as a function of time and we check if there is non-singular cosmological solution. From equations (2) (we recall that κ = 0) and using equation (5) we obtain 3(1 + α)(1 − 4kλ) 2 ρ0 2Ḣ + H = −k . (17) 1 − 3kλ(1 + α) 1 − 3kλ(1 + α) First, for simplicity, we consider the case α = 1/3. We choose the Rastall parameter Hb1 + a1 λ= , 4b1 kH with b1 , a1 constants. Then (17) becomes 2a1 Ḣ + 4a1 H 2 − b1 kρ0 H = 0, 7 and for ρ0 6= 0 has the general solution b1 k ρ0 H (t) =    . C b1 k ρ0 exp − b1 2a k p0 t 1 + 4 a 1 Let H0 be the Hubble parameter at an arbitrary fixed time t0 , namely H(t0 ) = H0 . Then we have H (t) b1 k ρ0 =   , H0 (b1 k ρ0 − 4 H0 a1 ) exp − b1 k ρ0 (t−t0 ) + 4 H0 a1 2a1 and presents a singularity at    1 H0 a 1 tc = − −b1 kρ0 t0 + 2a1 ln 4 . b1 kρ0 −b1 kρ0 + 4 H0 a1 Additionally, for a(t0 ) = a0 we have     b1 kρ0 (t−t0 ) 2 a (t) 4 H0 a1 exp 2a1 + b1 kρ 0 − 4 H0 a1 2 = . a0 b1 kρ0 So we have a big bang universe. After the big bang the universe evolves in an exponential manner. Now lets consider the particular case ρ0 = 0 (and still λ = (Hb1 + a1 )/(4b1 kH) and α = 1/3, this is the radiation era). Then (17) becomes Ḣ + 2H 2 = 0, and the solution satisfying the condition H(t0 ) = H0 is H (t) 1 a2 (t) = , = (1 + 2 ( t − t0 ))H0 H0 1 + 2 (t − t0 ) H0 a20 and presents a singularity at 2 H0 t0 − 1 tc = . 2H0 Now we consider a general α and that there is no pressure. Then by relation (5) we have that ρ0 = −αρ. Then from equations (2) (we recall that κ = 0) we obtain kρ0 Ḣ = , 2α and note that Rastall parameter λ does not appear in this equation. Now considering H(t0 ) = H0 and a(t0 ) = a0 we obtain a non–singular universe   k (t − t0 ) ρ0 + 2 H0 a (4 H0 α + kρ0 (t − t0 )) (t − t0 ) H(t) = , a (t) = a0 exp . 2a 4α For the particular case where ρ0 = 0 (so we have an empty space time, p = ρ0 = ρ = 0) the Hubble constant as time pass it remains constant (H(t) = H0 ) and there is an exponential growth in the scale factor a(t) = a0 exp(H0 (t − t0 )). 4. Darboux integrability and qualitative study of system (11) Now we present the proofs of Theorems 4 and 5. Proof of Theorem 4. First we consider that b 6= 0, 2. System (12) has the algebraic curves f1 = by + a = 0 and f2 = (2 − b)x2 + 2y + a = 0 invariant with cofactors K1 = bx and K2 = 2x respectively. Note that for b = 2 the two curves coincide. Therefore for b 6= 2 we have that bK2 − 2K1 = 0, and so system (12) admits the first integral (13) of Darboux type, see for more details Theorem 10 of Appendix 3.1. This proves statement (a). 8 The previous first integral is rational if b is a positive integer different from 0 and 2, and a polynomial first integral if b is a negative integer. In both cases all the orbits of system (12) are contained in algebraic curves.  For b = 2 we have f1 = f2 and system (12) has the exponential factor exp (2x2 − a)/(2y + a) with cofactor L1 = 2x. Hence the straight line f1 = 0 has multiplicity at least two, see [40]. In this case we have L1 − K1 = 0, and the first integral is given by (14) and is of Darboux type. This proves statement (b). For b = 0 and a 6= 0 system (12) admits the exponential factor exp(y) with cofactor L2 = ax. In this case, the straight line at infinity has multiplicity bigger than one, see [40]. So we have (a/2)K2 −L2 = 0, and system (12) has the first integral (15). Therefore statement (c) is proved. For b = a = 0 system (12) is the system ẋ = y + x2 , ẏ = 0, and clearly it has the first integral H(x, y) = y. Hence statement (d) is proved.  Remark 9. Note that the two invariant curves f1 = 0 and f2 = 0 of system (12) are actually the invariant curves g1 = 0 and g2 = 0 fof system (7). This can be easily checked taking into account the scaling of the variables (8), the relations (10) and the expression of system (11). 4.1. The finite singular points of systems (12). We recall that system (12) is the same as system (7) for the case where 1 − 4kλ 6= 0 and under the scaling of coordinates (8), see also Remark 3. System (12) has the finite singular point P0 = (0, 0) and it represents an empty flat model. p  For ab > 0 there are the two additional finite singular points P± = ± ab , − ab and represent spatially flat de-Sitter models. Note that the points P± belongs to the algebraic curves f1 = 0 and f2 = 0. For b 6= 2 the curve f2 = 0 is a parabola whereas for b = 2 coincides with the straight line f1 = 0. Note that a (1 − 4 kλ) ρ0 = . b 12λ2 k (1 + α) Since for system (12) we have that 1 − 4kλ 6= 0, we have that the only way that the points P± collide with the origin P0 is when a = 0 which corresponds to ρ0 = 0. We recall that in this case the equation of state becomes p = αρ. √ The Jacobian matrix at the origin has eigenvalues ± a. For a > 0 the origin is a hyperbolic saddle, whereas for a < 0 the point P0 it may be a focus or a center. By Theorem 4 system (12) with a < 0 is has a first integral defined in a neighborhood of the origin P0 , so P0 is a center. Therefore there is a neighborhood of P0 filled with periodic orbits. Since the singular point P0 represents an empty flat model which is a saddle if a > 0 and a center if a < 0, it follows statements (a) and (b) of Theorem 7. For a = b = 0 system (12) becomes ẋ = y + x2 , ẏ = 0 and after the change of time ds = (y + x2 )dt it writes as the linear system ẋ = 1, ẏ = 0. See Figure 25 for its phase portrait in the Poincaré disc. For a = 0 and b > 0 we apply Theorem 3.5 of [20] and the origin P0 is a union of one hyperbolic and one elliptic sector, we will denoted like HE. For a = 0 and b < 0 we apply Theorem 3.5 of [20] and we conclude that the origin P0 is a saddle. p p The Jacobian matrix of psystem (12) atpthe point P+ has eigenvalues 2 a/b and b a/b, and at P− has eigenvalues −2 a/b and −b a/b. Hence for a > 0 and b > 0 the singular point P+ is an unstable node, and P− is a stable node. So P+ represent an expanding flat de-Sitter model, whereas P− represents a contracting flat de-Sitter model. For a < 0 and b < 0 both singular points P± are saddles, so in this case de flat de-Sitter model is unstable. This proves statements (c), (d) and (e) of Theorem 7. 9 b P0 centre P0 saddle P+ unstable node P0 HE P- stable node P0 centre P0 saddle a P0 centre P+ saddle P0 saddle P0 saddle P- saddle Figure 3. The bifurcation diagram for the finite singular points of system (12). Hence the straight lines a = 0 and b = 0 are the bifurcation lines that characterize the local phase portraits at the finite singular points, see Figure 3. 4.2. The infinite singular points of system (12). In order to study the behaviour of the orbits of system (12) in a neighborhood of infinity we need to write the expression of system (12) into two charts, according to the Poincaré compactification, see chapter 5 of [20]. Hence at the first chart (U1 , F1 ) system (12) writes ż1 = −z1 2 z2 + az2 + bz1 − z1 , ż2 = −z2 (z1 z2 + 1) , (18) and for b 6= 1 the origin O1 is the only singular point of system (18). The Jacobian matrix of system (18) at the origin O1 has the eigenvalues −1 and b − 1. Hence, for b > 1 the origin is a hyperbolic saddle, whereas for b < 1 we have a stable node. For b = 1 then the infinity line z2 = 0 is filled of singular points and Theorem 11 of Appendix 3.2 applies. Now at the chart (U2 , F2 ) system (12) becomes ż1 = −az1 2 z2 − bz1 2 + z1 2 + z2 , ż2 = −z2 z1 (az2 + b) . The Jacobian matrix at the origin O2 is not identically zero and the two eigenvalues are both zero. Then from Theoem 3.5 of [20] we get: • If b = 0 and a < 0 the origin O2 is a saddle. • If b = 0 and a > 0 the origin O2 is a union of one hyperbolic and one elliptic sector, that we denote by HE. • If b ∈ (0, 1) the origin O2 is a saddle. • If b < 0 or b > 1 the origin O2 is a union of one hyperbolic and one elliptic sector. • If a = b = 0 we recall that system (12) becomes ẋ = x2 + y, ẏ = 0 and after the change of time ds = (x2 + y)dt it writes as the linear system ẋ = 1, ẏ = 0. The straight lines a = 0, b = 0 and b = 1 are the bifurcation curves for the infinite singular point of system (12), see Figure 4. Proof of Theorem 5. We use the Poincaré compactification to draw the phase portraits of system (12) in the Poincaré disc, see for more details chapter 5 of [20]. 10 b O1 saddle O1 saddle O1 saddle O2 HE O 2 HE O2 HE z2=0 full of singular points b=1 O1 stable node O1 stable node O2 HE O2 HE O1 stable node O2 HE O1 stable node O2saddle a O1 stable node O1 stable node O 2 HE O2 HE O1 stable node O 2 HE Figure 4. The bifurcation diagram for the infinite singular points of system (12). To system (12) we associate the polynomial vector field ∂ ∂ Xab = (y + x2 ) + x(a + by) . ∂x ∂y In Figure 2 we present all the phase portraits in the Poincaré disc of system (12) or equiva- lently of the compactified vector field p(Xab ). We have construct these phase portraits taking into account the local phase portraits at the finite and infinite singular points, the invariant algebraic curves, and the theorem of Markus–Neumann–Peixoto, see Theorem 13 of Appendix 3.3, and the Appendix 3.4 where in the different figures of that appendix we illustrate how we have obtained the phase portraits in the Poincaré disc for the different values of the parameters a and b. The key point of Theorem 13 for obtaining a phase portrait is to control its separatrices, where they born and where they die. For our system (12) the unique separatrices are the finite singular points, all the orbits at infinity and the separatrices of the hyperbolic sectors of the finite and infinite singular points. See for more details on the separatrices Appendix 3.3. Note that our system (12) has no limit cycles, because it is known that polynomial differential systems of degree two as our system (12) if they have a limit cycle then in the bounded region that they limited must be a focus, and our system (12) has no foci, for a proof of this result see [41]. Furthermore two phase portraits are not topological equivalent when it does not exist a homeomorphism to bring the separatrix configuration of one into the separatrix configuration of the other, see Theorem 13 in Appendix 3.3. In Figure 2 we have consider the bifurcation lines a = 0, b = 0, b = 1 and additionally the line b = 2. We recall that for b 6= 2 there are two invariant algebraic curves, the straight line f1 = 0 and the parabola f2 = 0, whereas for b = 2 we only have the straight line f1 = 0. Note that for b = 0 we only have the parabola f2 = 0. In the next table we put together the phase portraits of Figure 2 having the same number of separatrices S and the same number of canonical regions R, and we distinguish between them 11 the ones which are topological equivalent. S R 10 2 A12 = A13 6= A15 = A16 , 11 2 A10 11 3 A14 13 4 A5 = A6 = A7 , A5 6= A8 , A5 6= A9 , A5 6= A11 , A8 6= A9 , A8 6= A11 , A9 6= A11 17 4 A1 = A2 6= A4 17 5 A17 19 6 A3 So we obtain twelve non topological equivalent phase portraits. Additionally we must consider the phase portrats B0 , B1 , B2 and B3 having infinitely many separatrices. In summary systems (12) have 16 non–topological equivalent phase portraits.  5. Conclusions In this work we study the integrability and the dynamics of a Rastall gravity model filled with a perfect fluid that satisfies a linear equation of state. This model is given by a 2–dimensional polynomial differential system depending on four parameters. In general to prove the integrability or non integrability of a differential system depending on parameters is a very hard problem. In this work we prove that the model studied is Darboux integrable for all the values of the parameters. The existence and the expression of the first integral (see Theorem 4) shows that the equilibrium point localized at the origin of coordinates cannot be a focus and it must be a center. So there is a neighborhood of the origin full of periodic orbits, see section 4 and the phase portraits A12 , A13 , A14 , A15 , A16 , A17 and B3 in Figure 2. Moreover there is an invariant straight line that is involved in the expression of the first integral and is related with the zero energy condition see Remark 9 and also the phase portraits in the Appendix 6.4. In order to understand the dynamics of our model we study the stability not only of the finite singular points but also of the infinite ones. This is a big issue since we need to control the orbits that may escape to infinity or may come from infinity. This is the reason that we draw the global phase portraits in the Poincaré disc. In Figures 3 and 4 we present the bifurcation diagrams of the finite and infinite singular points respectively. Hence we can control the changes of the stability of the singular points when the parameters of the system vary. The physical and cosmological consequences of our dynamical study of the Rastall model here analyzed are given in Theorem 7 where we examine the stability of the flat de-Sitter models. Additionally, in Figure 2 we present all the non–topological equivalent global phase portraits of the model considered, see also Theorem 5 and Corollary 6. These phase portraits provide all the different possible dynamical behaviours of the model. Often the existence of several parameters, four parameters here, in a dynamical system makes very difficult its global qualitative study, and usually the numerical analysis is necessary in order to understand the dynamics of such models. There are many models in the literature relating Rastall gravity see for example Capone [38], Silva et al. [5], Oliveira et al. [42] and Cruz et al. [3] among others. We remark that we have used adequate transformations to reduce the number of parameters and we have used the polynomial structure of our model to achieve our results. Hence all the results in this work are not numeric, almost all of them are analytic, and some of them are algebraic. Note that the global phase portraits in Figure 2 are not numeric and provide all the possible dynamics of the Rastall model considered. An interesting future study could be the investigation of a similar model with a perfect fluid that satisfies a nonlinear equation of state. 12 We also check the existence of singular/non singular universe when Rastall parameter is varying with the time. We see that for a pressureless matter there is a non singular cosmology where for the radiation era there is a singular universe. 6. Appendices This section is composed by three appendices. 6.1. Appendix: Invariant curves, cofactors and Darboux first integrals. Consider the polynomial differential system ẋ = P (x, y), ẏ = Q(x, y). (19) Let f = f (x, y) be a polynomial in the variables x and y. The algebraic curve f (x, y) = 0 is an invariant algebraic curve of system (19) if ∂f ∂f P +Q = Kf, ∂x ∂y for some polynomial K = K(x, y) called the cofactor of the invariant algebraic curve f = 0. The polynomial structure of system (19) forces that the degree of the cofactor is at most one less than the degree of the polynomial system (19). Note, that the curve f = 0 is formed by the orbits of system (19), and consequently it is invariant under the flow of this system. There is another object called exponential factor which also provide cofactors and appear when invariant algebraic curves collide, i.e. when they have multiplicity larger than one. Let F (x, y) = exp(f (x, y)/g(x, y)) where f, g polynomials. Then F in an exponential factor of system (19) if ∂F ∂F P +Q = LF, ∂x ∂y for some polynomial L called the cofactor of the exponential factor F , see for more details [40]. A Darboux function is a function of the form, Y λ µj D= fi i Fj , where the fi = 0 are invariant algebraic curves of the system, and Fj ) are exponential factors. The following result started with Darboux [24], for the present version see for example [20]. Theorem 10. Suppose that a polynomial system (19) admits p irreducible invariant algebraic curves fi = 0 with cofactors Ki and q exponential factors Fj with cofactors Lj . Then there exist p P q P λi ’s and µj ’s in C not all zero such that λi Ki + µj Lj = 0, if and only if the function i=1 j=1 λ µ f1λ1 · · · fp p F1µ1 · · · Fq q (20) is a first integral of system (19). The first integrals of the form (20) are called the Darboux first integrals. 6.2. Appendix: Normally hyperbolic manifold. Let φt be a smooth flow on a manifold M and suppose C is a submanifold of M consisting entirely of equilibrium points for the flow. We say that C is normally hyperbolic if the tangent bundle to M over C splits into three subbundles T C, E s and E u invariant under dφt and satisfying (i) dφt contracts E s exponentially, (ii) dφt expands E u exponentially, (iii) T C = tangent boundle of C. 13 The following theorem holds, see [43]. Theorem 11. Let C be a normally hyperbolic submanifold of equilibrium Lpoints for φt .L Then there exist smooth stable and unstable manifolds tangent along C to E s T C and E u T C respectively. Furthermore, both C and the stable and unstable manifolds are permanent under small perturbation of the flow. 6.3. Appendix: Separatrix configuration. Our aim is to present the topological classifi- cation of the phase portraits of system (12). We use the definition of parallel flows given by Markus [44] and Neumann in [45]. Let φ be a C ω local flow on the two dimensional manifold R2 or R2 \ {0}. The flow (M, φ) is C k parallel if it is C ω –equivalent to one of the following ones: strip: (R2 , φ) with the flow φ defined by ẋ = 1, ẏ = 0. spiral: (R2 \ {0}, φ) with the flow φ defined by ṙ = r, θ̇ = 1. annular: (R2 \ {0}, φ) with the flow φ defined (in polar coordinates) by ṙ = 0, θ̇ = 1. We say that two polynomial vector fields Xa1 b1 and Xa2 b2 on R2 are topologically equivalent in the Poincaré disc D if there exists a homeomorphism of D preserving the infinity S1 and carrying orbits of the flow induced by p(Xa1 b1 ) into orbits of the flow induced by p(Xa2 b2 ). This homeomorphism should preserve or reverse simultaneously the sense of all orbits of the two compactified vector fields p(Xa1 b1 ) and p(Xa2 b2 ). The separatrices of a vector field p(Xab ) in the Poincaré disc D are (i) All the finite singular points of p(Xab ). (ii) All the separatrices of the hyperbolic sectors of the finite and infinite singular points of p(Xab ). (iii) All the orbits of p(Xab ) which are in the boundary S1 of the Poincaré disc. (iv) All the limit cycles of p(Xab ). Let Σ be the union of all separatrices of the flow (D, φ) defined by the compactified vector field p(Xab ) in the Poincaré disc D. Then Σ is a closed invariant subset of D. We say that every connected component of D \ Σ, with the restricted flow, is called a canonical region of φ. Theorem 12. Let φ be a C ω flow in the Poincaré disc with finitely many separatrices, and let Σ be the union of all its separatrices. Then the flow restricted to every canonical region is C ω parallel. For a proof of Theorem 12 see [46] and [45]. The separatrix configuration Σc of a flow (D, φ) is the union of all the separatrices Σ of the flow together with an orbit belonging to each canonical region. The separatrix configuration Σc of the flow (D, φ) is topologically equivalent to the separatrix configuration Σ̃c of the flow (D, φ̃) if there exists a homeomorphism from Σc to Σ̃c which transforms orbits of Σc into orbits of Σ̃c . Theorem 13. Let (D, φ) and (D, φ̃) be two compactified Poincaré flows with finitely many sep- aratrices coming from two polynomial vector fields (12). Then they are topologically equivalent if and only if their separatrix configurations are topologically equivalent. The reader can find the proof of Theorem 13 in [44, 45, 47]. According to Theorem 13 the topological classification of the phase portraits in the Poincaré disc of a family of polynomial differential system with finitely many separatrices can be done by considering their separatrix configurations. 6.4. Appendix: The local and the phase portraits of system (12). In this appendix and through the different figures we explain how the phase portraits in the Poincaré disc that appear in Figure 2 are obtained from the local phase portraits of the finite and infinite singular points, taking also into account their invariant algebraic curves and their separatrix configurations. 14 . . . ... . . ... . P- P0 P+ . . Figure 5. The local and global the phase portraits of system (12) for a > 0 and b > 2. In both pictures we have drawn the straight line and the parabola. Here S = 17 and R = 4. . . . ... . . ... . P- P0 P+ . . Figure 6. The local and global the phase portraits of system (12) for a > 0 and b = 2. The parabola does not exist and in both pictures we have drawn the straight line. Here S = 17 and R = 4. . . . ... . . ... . P- P0 P+ . . Figure 7. The local and global the phase portraits of system (12) for a > 0 and b ∈ (1, 2). In both pictures we have drawn the straight line and the parabola. We have S = 19 and R = 6. Acknowledgments J. Llibre is partially supported by the Agència de Gestió d’Ajuts Universitaris i de Re- cerca grant 2017SGR1617, and the H2020 European Research Council grant MSCA-RISE-2017- 777911. J. Llibre and C. Pantazi are also supported by the Ministerio de Ciencia, Innovación y Universidades, Agencia Estatal de Investigación grant MTM2016-77278-P (FEDER). C. Pan- tazi is additionally partially supported by the Catalan Grant 2017SGR1049 and the Spanish MINECO-FEDER Grant PGC2018-098676-B-I00/AEI/FEDER/UE. 15 . . . P- P0 P+ . . . Figure 8. The local and global the phase portraits of system (12) for a > 0 and b = 1. In both pictures we have drawn the straight line and the parabola. Here the infinity is filled by singular points. . . . ... . . ... . P- P0 P+ . . Figure 9. The local and global the phase portraits of system (12) for a > 0 and b ∈ (0, 1). In both pictures we have drawn the straight line and the parabola. Here S = 17 and R = 4. . . . . . . . . . . Figure 10. The local and global the phase portraits of system (12) for a > 0 and b = 0. In this case the straight line does not exist and in both pictures we have drawn the parabola. Here S = 13 and R = 4. References [1] F. Darabi, H. Moradpour, I. Licata, Y. Heydarzade, and C. Corda. Einstein and rastall theories of gravitation in comparison. Eur. Phys. J. C, 78:78–25, 2018. [2] J.C. Fabris O.F. Piattella C.E.M. Batista, M.H. Daouda and D. C. Rodrigues. Rastall cosmology and the λ cdm model. Phys. Rev. D, 85:084008, 2012. [3] M. Cruz, S. Lepe, and G. Morales-Navarrete. thermodynamics revision of rastall gravity. Class. Quant. Grav., 36:225007, 2019. 16 . . . . . . . . . . Figure 11. The local and global the phase portraits of system (12) for a > 0 and b < 0. In both pictures we have drawn the straight line and the parabola. Here S = 13 and R = 4 . . . . . . . . . . . Figure 12. The local and global the phase portraits of system (12) for a = 0 and b < 0. In both pictures we have drawn the straight line and the parabola. Note that the origin P0 belongs to the straight line. Here S = 13 and R = 4. . . . ... . . ... . P- P+ . . Figure 13. The local and global the phase portraits of system (12) for a < 0 and b < 0. In both pictures we have drawn the straight line and the parabola. We have S = 17 and R = 5. [4] N.D. Birrell and P.C.W. Davies. Quantum Fields in Curved Space. Cambridge University Press, 1982. [5] G.F. Silva, O.F. Piattella, J.C. Fabris, L. Casarini, and T.O. Barbosa. Bouncing solutions in rastall’s theory with a barotropic fluid. Gravit. Cosmol., 19:156–162, 2013. [6] G.W. Gibbons and S.W. Hawking. Cosmological event horizons, thermodynamics, and particle creation. Phys. Rev. D, 15:2738, 1977. [7] L.H. Ford. Gravitational particle creation and inflation. Phys. Rev. D, 35:2955, 1987. 17 . . . . . . . . . . Figure 14. The local and global the phase portraits of system (12) for a < 0 and b = 0. In this case the invariant straight line does not exist. In both pictures we have drawn the parabola. Here S = 10 and R = 2. . . . . . . . . . . Figure 15. The local and global the phase portraits of system (12) for a < 0 and b ∈ (0, 1). In both pictures we have drawn the straight line and the parabola. Here S = 10 and R = 2. . . . Figure 16. The local and global the phase portraits of system (12) for a < 0 and b = 1. In both pictures we have drawn the straight line and the parabola. The infinity is filled by singular points. [8] C.E.M. Batista, J.C. Fabris, O.F. Piattella, and A.M. Velasquez-Toribio. Observational constraints on rastalls cosmology. Eur. Phys. J. C, 73:2425, 2013. [9] W. Khyllep and J. Dutta. Linear growth index of matter perturbations in rastall gravity. Phys. Lett. B, 797:134796, 2019. [10] Y. Heydarzade, H. Moradpour, and F. Darabi. Black hole solutions in rastall theory. Can. J. Phys., 95:12531256, 2017. [11] H. Moradpour. Thermodynamics of flat flrw universe in rastall theory. Phys. Lett. B, 757:187, 2016. 18 . . . . . . . . . . Figure 17. The local and global the phase portraits of system (12) for a < 0 and b ∈ (1, 2). In both pictures we have drawn the straight line and the parabola. We have S = 11 and R = 3. . . . . . . . . . . Figure 18. The local and global the phase portraits of system (12) for a < 0 and b = 2. In this case the parabola does not exist. In both pictures we have drawn the straight line. Here S = 10 and R = 2. . . . . . . . . . . Figure 19. The local and global the phase portraits of system (12) for a < 0 and b > 2. In both pictures we have drawn the straight line and the parabola. Here S = 10 and R = 2. [12] H. Moradpour, Y. Heydarzadeand F. Darabi, and I. G. Salako. A generalization to the rastall theory and cosmic eras. Eur. Phys. J. C, 77:259, 2017. [13] L.L. Smalley. Rastall’s and related theories are conservative gravitational theories although physically inequivalent to general relativity. J. Phys. A, 16:2179–2185, 1983. [14] S. Nojiri and S.D. Odintsov. Gravity assisted dark energy dominance and cosmic acceler- ation. Phys. Lett. B, 599:137–142, 2004. [15] T. Koivisto. A note on covariant conservation of energymomentum in modified gravities. Class. Quantum Gravit., 23:4289, 2006. 19 . . . . . . . . . . Figure 20. The local and global the phase portraits of system (12) for a = 0 and b > 2. In both pictures we have drawn the straight line and the parabola. Here the origin P0 belongs to the straight line. We have S = 13 and R = 4. . . . . . . . . . . Figure 21. The local and global the phase portraits of system (12) for a = 0 and b = 2. Here the parabola does not exist and in both pictures we have drawn the straight line. Note that the origin P0 belongs to the straight line. Here S = 11 and R = 2. . . . . . . . . . . Figure 22. The local and global the phase portraits of system (12) for a = 0 and b ∈ (1, 2). In both pictures we have drawn the straight line and the parabola. The origin P0 belongs to the straight line. Here S = 13 and R = 4. [16] T. Harko and F.S.N. Lobo. Generalized curvature-matter couplings in modified gravity. Galaxies, 2:410–465, 2014. [17] M. Visser. Rastall gravity is equivalent to einstein gravity. Phys. Lett. B, pages 83–86, 2018. [18] P. Rastall. Generalization of the einstein theory. Phys. Rev. D, 6:3357–3359, 1972. [19] S. Hansraj, A. Banerjee, and P. Channuie. Impact of the rastall parameter on perfect fluid spheres. Annals of Phys., 400:320–345, 2019. 20 . . Figure 23. The local and global the phase portraits of system (12) for a = 0 and b = 1. In both pictures we have drawn the straight line and the parabola. Note that the origin P0 belongs to the straight line. Additionally the infinity is filled by singular points. . . . . . . . . . . Figure 24. The local and global the phase portraits of system (12) for a = 0 and b ∈ (0, 1). In both pictures we have drawn the straight line and the parabola. The origin P0 belongs to the straight line. Here S = 13 and R = 4. . . . . Figure 25. The local and global the phase portraits of system (12) for a = b = 0. In both pictures we have drawn the straight line and the parabola. The origin P0 belongs to the straight line. Additionally, in this case the parabola y + x2 = 0 is fulfilled by singular points. [20] F. Dumortier, J. Llibre, and J. C. Artés. Qualitative theory of planar polynomial systems. Springer, 2006. [21] P.D. Lax. Integrals of nonlinear equations of evolution and solitary waves. Commun. Pure Appl. Math., 21:467–490, 1968. [22] P.J. Olver. Applications of Lie groups to differential equations. Springer, 1986. [23] F. Cantrijn and W. Sarlet. Generalizations of noether’s theorem in classical mechanics. SIAM Rev., 23:467–494, 1981. 21 [24] G. Darboux. Mémoire sur les équations différentielles algébriques du premier ordre et du premier degré (mélanges). Bull. Sci. Math. 2ème Série, 2:60–96, 123–144, 151–200, 1878. [25] T.C. Bountis, A. Ramani, B. Grammaticos, and B. Dorizzi. On the complete and partial integrability of non–hamiltonian systems. Phys. A, 128:268–288, 1984. [26] J.J. Morales-Ruiz. Differential Galois Theory and Non-integrability of Hamiltonian Sys- tems. Progress in Mathematics 179. Birkhäuser, 1999. [27] D. Das, S. Dutta, and S. Chakraborty. Cosmological consequences in the framework of generalized rastall theory of gravity. Eur. Phys. J. C, 78:810, 2018. [28] E. Babichev, V. Dokuchaev, and Yu. Eroshenko. Dark energy cosmology with generalized linear equation of state. Class. Quant. Grav., 22:143–154, 2005. [29] P. Rastall. A theory of gravity. Can. J. Phys., 54:66–75, 1976. [30] J. Wainwright and G. F. R. Ellis. Dynamical Systems in Cosmology. Cambridge University Press, 1997. [31] G. Kofinas, G. Leon, and E. N. Saridakis. Class. Quantum Grav., 31:175011, 2014. [32] L. Heisenberg, R. Kimura, and K. Yamamoto. Cosmology of the proxy theory to massive gravity. Phys. Rev.D, 89:103008, 2014. [33] N. Goheer, R. Goswami, and P. K. S. Dunsby. Dynamics of f(r)–cosmologies containing einstein static models. Class. Quant. Grav., 26:105003, 2009. [34] A. Ganguly, R. Gannouji, R. Goswami, and S. Ray. Global structure of black holes via the dynamical system. Class. Quantum Grav., 32:105006, 2015. [35] M. Cruz, A. Ganguly, R. Gannouji, G. Leon, and E. N. Saridakis. Global structure of static spherically symmetric solutions surrounded by quintessence. Class. Quantum Grav., 34:125014, 2017. [36] H. Moradpour and I. G. Salako. Thermodynamic analysis of the static spherically sym- metric field equations in rastall theory. Advances in High Energy Physics, vol. 2016:5p., 2016. [37] F.F. Yuan and P. Huang. Emergent cosmic space in rastall theory. Class. Quantum Grav., 34:077001, 2017. [38] M Capone, V. F. Cardone, and M. L. Ruggiero. The possibility of an accelerating cosmology in rastall’s theory. J. Phys. Conf. Ser., 222:012012, 2010. [39] A.H.Ziaie, H.Moradpour, and S.Ghaffari. Gravitational collapse in rastall gravity. Physics Letters B, 793:276–280, 2019. [40] C. Christopher, J. Llibre, and J.V. Pereira. Multiplicity of invariant algebraic curves in polynomial vector fields. Pacific J. Math., 229:63–117, 2007. [41] W.A. Coppel. A survey of quadratic systems. J. Diff. Eqns., 2:293–304, 1966. [42] A.M. Oliveira, H.E.S. Velten, J.C. Fabris, and L. Casarini. Neutron stars in Rastall Gravity. Phys. Rev. D, 92:044020, 2015. [43] M. W. Hirsch, C. C. Pugh, and M. Shub. Invariant manifolds, Lecture Notes in Mathe- matics. Springer–Verlag, 1977. [44] L. Markus. Global structure of ordinary differential equations in the plane. Trans. Amer. Math Soc., 76:127–148, 1954. [45] D.A. Neumann. Classification of continuous flows on 2–manifolds. Proc. Amer. Math. Soc., 48:73–81, 1975. [46] J. Llibre W. Li, M. Nicolau, and X. Zhang. On the differentiability of first integrals of two dimensional flows. Proc. Amer. Math. Soc., 130:2079–2088, 2002. [47] M.M. Peixoto. Dynamical Systems. Proccedings of a Symposium held at the University of Bahia, 389–420. Acad. Press, 1973. 1 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Edifici C, 08193 Bel- laterra, Barcelona, Catalonia, Spain E-mail address: jllibre@mat.uab.cat 22 2 Departament de Matemàtiques, Universitat Politècnica de Catalunya, (EPSEB), Av. Doctor Marañón, 44–50, 08028 Barcelona, Spain E-mail address: chara.pantazi@upc.edu 23