This is a preprint of: “Global qualitative dynamics of the Brusselator system”, Jaume Llibre, Clàudia Valls, Math. Comput. Simulation, vol. 170, 107–114, 2020. DOI: [10.1016/j.matcom.2019.01.017] GLOBAL QUALITATIVE DYNAMICS OF THE BRUSSELATOR SYSTEM JAUME LLIBRE1 AND CLÀUDIA VALLS2 Abstract. We study the dynamics of the Brusselator model by de- scribing the global flow of this differential system in the Poincaré disc, and we prove that its unique equilibrium point is a global attractor. 1. Introduction and statement of the results We consider the differential system ẋ = a + x2 y − (b + c)x, (1) ẏ = bx − x2 y, called the Brusselator model where x and y are the dimensionless concen- tration of some species and the parameters a, b, c are all positive. Such differential system appears in several branches of the sciences, mainly in chemistry because it exhibits kinetics of model trimolecular irreversible re- actions (see [3, 4] for details). The first integrals of system (1) in function of its parameters were studied in [5]. We will determine the qualitative behavior of all the solutions of system (1), by describing its phase portraits in the Poincaré disc in function of its parameters. See the appendix for the definitions and the basic results that we use, and in particular for the definition of the Poincaré disc. We say that two vector fields on the Poincaré disc are topologically equiva- lent if there exists a homeomorphism of the Poincaré disc which sends orbits to orbits preserving or reversing their orientation. Our main result is the following one. Theorem 1. The global phase portrait of system (1) in the Poincaré disc is topologically equivalent to the one described in Figure 1. Note that its unique finite singular point is a global attractor. The proof of Theorem 1 is given in section 2, the more delicate step in this proof is to show the non–existence of periodic solutions. We have also 2010 Mathematics Subject Classification. 34C05, 34C07, 34C08. Key words and phrases. Brusselator system, nilpotent saddle, phase portrait, Poincaré compactification, global attractor. 1 http://www.gsd.uab.cat 2 J. LLIBRE AND C. VALLS Figure 1. This phase portrait in the Poincaré disc has 16 sepa- ratrices and 3 canonical regions. included an appendix with all the definitions and basic results related with the Poincaré compactification. 2. Proof of Theorem 1 By the change of coordinates and reparametrization of time of the form x → aX/b, y → b2 Y /a, t → τ /a, we can write system (1) as ẋ = (2) ẏ = x − x2 y, where k = c/b > 0. The dot denotes derivative with respect to the new variable τ and we have renamed the new variables (X, Y ) as (x, y). 2.1. The finite singular points. The unique finite singular point is (1/k, k). The Jacobian matrix M of system (2) at this singular point is   1 − k 1/k2 M= . −1 −1/k2 The eigenvalues of M are p −1 + k2 − k3 ± (1 − k2 + k3 )2 − 4k3 . 2k2 We study the zeroes of 1 − k2 + k3 = 0 and we get two complex solutions and the following real one 1 1 2 1/3 1  1 √ 1/3 k̄ = − √ − (25 − 3 69) < 0. 3 3 25 − 3 69 3 2 Since we are interested for the values of k > 0, and then we have that −1 + k2 − k3 < 0, the equilibrium point is always stable. Moreover, the determinant of M is 1/k > 0, and so the finite singular point is either a BRUSSELATOR SYSTEM 3 hyperbolic node or a hyperbolic focus. We compute the zeroes of (1 − k2 + k3 )2 − 4k3 = 0, and we get four complex zeroes and two real ones, which are 1 5 2 1/3 1  1 √ 1/3 k1 = − √ + (11 + 3 69) = 0.5698.. > 0 3 3 11 + 3 69 3 2 and √ 1 1  47 3 93 1/3 1  1 √ 1/3 k2 = + − + (47 + 3 93) = 2.1478.. > k1 . 3 3 2 2 3 2 So if k ∈ (0, k1 ] ∪ [k2 , +∞) the singular point (1/k, k) is a hyperbolic stable node, and if k ∈ (k1 , k2 ) it is a hyperbolic stable focus. Now we shall show that this system has no periodic solutions, the proof will use the Dulac–Bendixson criterion, i.e.: Theorem 2. If there exists a C 1 function B(x, y) in a simply connected region R such that ∂(BP )/∂x + ∂(BQ)/∂y has constant sign and is not zero in a simply connected domain R, then the C 1 differential system ẋ = P (x, y), ẏ = Q(x, y) does not have periodic orbits in R. For a proof of this criterion see Theorem 7.12 of [2]. Proposition 3. The polynomial differential system (2) has no periodic so- lutions. Proof. First note that for k > 0 the singular point (1/k, k) belongs to the positive quadrant x > 0 and y > 0. We have that y (3) ẋ|{x=1/(1+k),y>0} = > 0 and ẏ|{x>0,y=0} = x > 0. (1 + k)2 If there exists a periodic solution, since it must surrounds the unique singular point (1/k, k) (see for more details Theorem 1.31 of [2]), from (4) the periodic solution must be contained in the simply connected region 1 R = {x > , y > 0}. 1+k For any positive integer n consider the function 1 B(x, y) = . (x − 1/(1 + k))2 We apply the Dulac–Bendixson criterion to the differential system (2) in the region R with the function B(x, y), note that B(x, y) > 0 in R. Then for P (x, y) = 1 − (1 + k)x + x2 y and Q(x, y) = x − x2 y, we get that ∂(BP ) ∂(BQ) g(x, y) = + ∂x ∂y   (k + 1)2 k2 x − k x3 − 2x + 1 − x3 + x2 − 2xy + x − 1 = . ((k + 1)x − 1)3 4 J. LLIBRE AND C. VALLS Note that in the region R the curve g(x, y) = 0 is equivalent to the curve  ((k + 1)x − 1) k + 1 − x2 (4) f (x, y) = y − = 0, 2x √ defined for x > k + 1. We claim that the curve f (x, y) = 0 in the region R is transversal with respect to the flow of the differential system (2), i.e. ∂f ∂f P+ Q 6= 0 in R. ∂x ∂y f =0 Due to this claim a periodic orbit cannot intersect the curve f (x, y) = 0 in the region R, therefore the function g(x, y) does not vanish on a periodic orbit, and by the Dulac–Bendixson criterion the differential system (2) has no periodic solutions, as we want to prove. Now we prove the claim. We have that ∂f ∂f h(x) P+ Q =− , ∂x ∂y f =0 4x2 where h(x) = 2(k + 1)2 x7 − 5(k + 1)x6 + (−2k3 − 6k2 − 6k + 1)x5 + 8(k + 1)2 x4 −4(2k + 3)x3 + (k3 + 3k2 + 3k + 3)x2 − 3(k + 1)2 x + 2(k + 1). Using the Sturm sequence associated to the polynomial√h(x) it follows that the polynomial h(x) has no real roots in the interval ( k + 1, ∞). Conse- quently the claim is proved. See [9] for the definition and the properties of a Sturm sequence, we note that the algebraic manipulator mapple allows to compute the Sturm sequence associated to a polynomial. This completes the proof of the proposition.  We have all the finite information. In the next subsection we study the infinite singular points. 2.2. The infinite singular points. On the local chart U1 system (2) be- comes u̇ = −u − u2 + v 2 + (1 + k)uv 2 − uv 3 , v̇ = v(−u + (1 + k)v 2 − v 3 ). The singular points (u, v) of this system on v = 0 must satisfy −u(1+u) = 0. So we have two infinite singular points in the local chart U1 : q1 = (u, v) = (0, 0) and q2 = (u, v) = (−1, 0). Computing the Jacobian matrix M1 at the point q1 we get   −1 0 M1 = . 0 0 BRUSSELATOR SYSTEM 5 So, q1 is semi-hyperbolic. Using Theorem 2.19 in [2] we get that it is a saddle. On the other hand, the Jacobian matrix M2 at the point q2 is   1 0 M1 = , 0 1 and so the point q2 is a hyperbolic unstable node. On the local chart U2 system (2) becomes u̇ = u2 + u3 − (1 + k)uv 2 + v 3 − u2 v 2 , (5) v̇ = uv(u − v 2 ). So the origin of U2 is an equilibrium point whose linear part is identically zero, and by the Poincaré–Hopf Theorem (see for instance chapter 6 of [2]) its topological index must be 0. We apply the blow-up techniques to study it, see for instance [1]. We first do a vertical blow-up. For this, we consider the new variables (u, w) where w = v/u. In these new variables system (5) writes u̇ = u2 (1 + u − u(1 + k)w2 − u2 w2 + uw3 ), ẇ = −uw(1 − (1 + k)uw2 + uw3 ). We eliminate the common factor u by making a rescaling of time and we get u̇ = u(1 + u − u(1 + k)w2 − u2 w2 + uw3 ), (6) ẇ = −w(1 − (1 + k)uw2 + uw3 ). System (6) has a unique singular point on u = 0 which is (u, w) = (0, 0). Computing the eigenvalues of the Jacobian matrix of system (6) at the origin we get that they are 1 and −1. Hence, the origin is a saddle. Going back through the changes of variables passing we do not have all the information for determining the local phase portrait at the origin of the local chart U2 . So we must do a horizontal blow-up. Hence, we consider the new variables (w, v) where w = u/v. In these new variables system (5) can be written as ẇ = v(v + w2 − (1 + k)wv), v̇ = wv 3 (w − v). We eliminate the common factor v by making a rescaling of time and we get ẇ = v + w2 − (1 + k)wv, (7) v̇ = wv 2 (w − v). System (7) has a unique singular point on v = 0 which is the origin. Com- puting the eigenvalues of the Jacobian matrix at (w, v) = (0, 0) we get that the origin is a nilpotent singular point. In order to determine its local be- havior we need to do another blow-up. We thus consider another vertical 6 J. LLIBRE AND C. VALLS V U Figure 2. The local phase portrait at the origin of the local chart U2 . blow-up of the form (w1 , v) where w1 = w/v. In these new variables (7) writes ẇ1 = 1 − w1 v(1 + k) + w12 v + w12 v 3 − w13 v 3 , v̇ = −w1 v 4 (1 − w1 ). Note that this system has no singular points. Going back through the changes of variables from (w1 , v) until (u, v), and taking also into account the information on the horizontal blow-up we get that the origin of the local chart U2 in the variables (u, v) is the degenerate saddle node of Figure 2.2. 2.3. The phase portrait in the Poincaré disc. Taking into account the local phase portraits at the finite and infinite singular points, together with the behavior of the vector field associated to system (2) on the axes of coordinates, and that in view of Proposition 3 system (2) has no limit cycles surrounding the finite singular point, the phase portrait in the Poincaré disc of system (2) is the one given in Figure 1. From this global phase portrait it follows that the equilibrium point (1/k, k) is a global attractor. This concludes the proof of the theorem. Appendix: Basic results and the Poincaré disc In this section we present some basic results and notations which are necessary for stating and proving our results. 2.4. Singular points of differential systems in R2 . Consider a differ- ential system ẋ = P (x, y), ẏ = Q(x, y) in R2 . A singular point or an equilibrium point of the differential system (2) is a point (x0 , y0 ) ∈ R2 such that P (x0 , y0 ) = Q(x0 , y0 ) = 0. BRUSSELATOR SYSTEM 7 We recall that a hyperbolic singular point is a singular point (x0 , y0 ) such that the eigenvalues of the Jacobian matrix  ∂P ∂P  (x , y ) (x0 , y0 )  ∂x 0 0 ∂y  (8)  ∂Q ∂Q  (x0 , y0 ) (x0 , y0 ) ∂x ∂y have nonzero real part. A semi-hyperbolic singular point is a singular point such that one of the eigenvalues of the matrix (8) is zero. A nilpotent singular point is a singular point such that the eigenvalues of the matrix (8) are both zero but the matrix is not identically zero. Finally, a linearly zero singular point is a singular point for which the matrix (8) is identically zero. For details on the local phase portraits of the hyperbolic, semi-hyperbolic and nilpotent singular points see [2]. On the other hand, the local phase portraits of the linearly zero singular points must be studied using the change of variables called blow-ups, see for instance [1]. 2.5. Poincaré compactification. Let pi (x1 , x2 ) be a real polynomial in the variables x1 , x2 , and let X = (p1 , p2 ) be a polynomial vector field of degree d in Rn , being d the maximum of the degrees of the polynomials pi for i = 1, 2 (see all the details in Chapter 5 of [2]). P The Poincaré sphere is defined as S2 = {y = (y1 , y2 , y3 ) ∈ R3 : 3i=1 yi2 = 1} and its tangent space at the point y ∈ S2 is denoted by Ty S2 . We identify the space R2 where is defined the polynomial vector field X with the tangent space T(0,0,1) S2 . We define the central projection f : T(0,0,1) S2 −→ S2 as follows: to each point q ∈ T(0,0,1) S2 the central projection f associates to q the two intersection points of the sphere S2 with the straight line which connects the points q with the origin of coordinates. Note that the infinity of R2 ≡ T(0,0,1) S2 corresponds to the equator S1 = {y ∈ S2 : y3 = 0} of S2 . This central projection f provides two copies Df ◦ X of the polynomial vector field X in S2 , one in the northern hemisphere and the other in the southern. Denote by X ′ the vector field defined by these two copies of X into the sphere S2 minus its equator S1 . We can extend the vector field X ′ on S2 \S1 to an analytic vector field p(X) on S2 defining p(X) = y3d+1 X ′ . Let π : R3 −→ R2 the projection π(y1 , y2 , y3 ) = (y1 , y2 ). Then the projection by π of the closed northern hemisphere of S2 is the Poincaré disc D. Note that the interior of D is diffeomorphic to R2 and its boundary S1 is the infinity of R2 , and π(p(X)) is the extension of the polynomial vector field X to the Poincaré disc D. For computing the analytic expression of p(X) we consider the sphere S2 as a differentiable manifold. We consider Ui = {y ∈ S2 : yi > 0} and Vi = {y ∈ S2 : yi < 0} for i = 1, . . . , 3. The corresponding coordinates maps 8 J. LLIBRE AND C. VALLS are given by Fi : Ui → R2 and Gi : Vi → R2 by 1 Fi (y) = Gi (y) = (yj1 , yj2 ) yi with 1 ≤ j1 < j2 ≤ 3 and jk 6= i for k = 1, 2. After a rescaling in the independent variable in the local chart (U1 , F1 ) the expression of p(X) is z2d (−z1 p1 + p2 , −z2 p1 ), where pi = pi (1/z2 , z1 /z2 ) for i = 1, 2. In a similar manner we can deduce the expression of p(X) in the local chart (U2 , F2 ). This is z2d (−z1 p2 + p1 , −z2 p2 ), pi = pi (z1 /z2 , 1/z2 ), for i = 1, 2. In the chart (Vi , Gi ) the expression of p(X) is the same than in the chart (Ui , Fi ) multiplied by (−1)d for i = 1, 2. 2.6. Topological equivalent polynomial vector fields. We recall that two polynomial vector fields X and Y on R2 are topologically equivalent if there is a homeomorphism on the Poincaré disc D preserving the infinity S1 and carrying orbits of the flow of π(p(X)) into trajectories of the flow of π(p(Y )), either preserving or reversing the sense of all the orbits. A separatrix of the Poincaré compactification π(p(X)) is a trajectory which is either an equilibrium point, or a limit cycle, or an orbit on the boundary of a hyperbolic sector at an equilibrium point, finite or infinity, or any orbit contained at the infinity S1 . The closed set (see Neumann [7]) formed by all separatrices of π(p(X)) is denoted by ΣX . A canonical region of π(p(X)) is an open connected component of D \ ΣX . The union of ΣX plus one orbit chosen from each canonical region is the separatrix configuration of π(p(X)), denoted by Σ′X . Two separatrix configurations Σ′X and Σ′Y are equivalent if there is a homeomorphism in B preserving the infinity S1 carrying orbits of Σ′X into orbits of Σ′Y , either preserving or reversing the sense of all orbits. Markus [6], Neumann [7] and Peixoto [8] characterized the topologically equivalence between two Poincaré compactified vector fields by showing that two Poincaré compactified polynomial vector fields π(p(X)) and π(p(Y )) with finitely many separatrices are topologically equivalent if and only if their separatrix configurations Σ′X and Σ′Y are equivalent. Acknowledgements The first author is partially supported by a FEDER-MINECO grant MTM2016-77278-P, a MINECO grant MTM2013-40998-P, and an AGAUR BRUSSELATOR SYSTEM 9 grant number 2014SGR-568. The second author is partially supported by FCT/Portugal through UID/MAT/04459/2013. References [1] M.J. Álvarez, A. Ferragut and X. Jarque, A survey on the blow up technique, Int. J. Bifurcation and Chaos 21 (2011), 3103–3118. [2] F. Dumortier, J. Llibre and J.C. Artés, Qualitative Theory of Planar Differential Systems, Springer Verlag, New York, (2006). [3] P. Gray and S.K. Scott, Chemical Oscillations and Instabilities, Non-linear Chem- ical Kinetics (Reihe: International Series of Monographs on Chemistry, Vol. 21). Von P. Gray and S. K. Scott. Clarendon Press, Oxford 1990. [4] A.I. Lavrovaa, E.B. Postnikovb and Y.M. Romanovsky, Brusselator: an ab- stract chemical reaction?, Phys. Usp. 52 (2009), 1239?-1244. [5] J. Llibre and C. Valls, Liouvillian and analytic first integrals for the Brusselator system, J. of Nonlinear Mathematical Physics 19 (2012), 1250016–pp. 8. [6] L. Markus, Global structure of ordinary differential equations in the plane: Trans. Amer. Math Soc. 76 (1954), 127–148. [7] D. A. Neumann, Classification of continuous flows on 2–manifolds, Proc. Amer. Math. Soc. 48 (1975), 73–81. [8] M.M. Peixoto, Dynamical Systems. Proccedings of a Symposium held at the Uni- versity of Bahia, 389–420, Acad. Press, New York, 1973. [9] J. Stoer and R. Bulirsh, Introduction to Numerical Analysis, Texts in Applied Mathematics, Springer–Verlag, 2002. 1 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain E-mail address: jllibre@mat.uab.cat 2 Departamento de Matemática, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais 1049–001, Lisboa, Portugal E-mail address: cvalls@math.ist.utl.pt