mathematics Article On the Zero-Hopf Bifurcation of the Lotka–Volterra Systems in R3 Maoan Han 1 , Jaume Llibre 2, * and Yun Tian 1 1 Department of Mathematics, Shanghai Normal University, 100 Guilin Road, Xuhui District, Shanghai 200234, China; mahan@shnu.edu.cn (M.H.); ytian22@shnu.edu.cn (Y.T.) 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain * Correspondence: jllibre@mat.uab.cat  Received: 5 May 2020; Accepted: 5 July 2020; Published: 12 July 2020  Abstract: Here we study 3-dimensional Lotka–Volterra systems. It is known that some of these differential systems can have at least four periodic orbits bifurcating from one of their equilibrium points. Here we prove that there are some of these differential systems exhibiting at least six periodic orbits bifurcating from one of their equilibrium points. We remark that these systems with such six periodic orbits are non-competitive Lotka–Volterra systems. The proof is done using the algorithm that we provide for computing the periodic solutions that bifurcate from a zero-Hopf equilibrium based in the averaging theory of third order. This algorithm can be applied to any differential system having a zero-Hopf equilibrium. Keywords: Lotka–Volterra differential systems; periodic orbit; Hopf bifurcation; averaging theory MSC: Primary 34C07; 34C08; 37G15 1. Introduction and Statement of Results An equilibrium point of a 3-dimensional autonomous differential system having a pair of purely imaginary eigenvalues and a zero eigenvalue is a zero-Hopf equilibrium. A 2-parameter unfolding of a 3-dimensional autonomous differential system with a zero-Hopf equilibrium is a zero-Hopf bifurcation. More precisely, when the two parameters of the unfolding are zero we have an isolated zero-Hopf equilibrium, and the dynamics of the unfolding is complex and sometimes chaotic in a small neighborhood of this isolated equilibrium when we vary the two parameters in a small neighborhood of the origin, see for more details [1–8] and references quoted there. A Lotka–Volterra system in R3 with coordinates ( x1 , x2 , x3 ) is a quadratic polynomial differential system of the form 3 dxi = xi ri − ∑ aij x j , i = 1, 2, 3,  (1) dt j =1 where the dot denotes derivative with respect to the independent variable t, usually called the time, and the ri ’s and the aij ’s are parameters. Many natural phenomena can be modeled by the Lotka–Volterra systems, starting in biology with the time evolution of conflicting species that now continuing being studied intensively see [9–20], later on problems of plasma physics [21], or problems in hydrodynamics [22], · · · . It is known that Lotka–Volterra systems can exhibit zero-Hopf equlibria, see for instance [23]. Then a natural question is if we perturbed a Lotka–Volterra system (1) having a zero-Hopf equilibrium point inside the class of all Lotka–Volterra systems (1) how many periodic orbits can bifurcate from such an equilibrium? Mathematics 2020, 8, 1137; doi:10.3390/math8071137 www.mdpi.com/journal/mathematics Mathematics 2020, 8, 1137 2 of 14 Note that the unfolding of a Lotka–Volterra system (1) with a zero-Hopf equilibrium needs at least a 3-parameter family. Arnold [24] in 1973 proposed to investigate bifurcations of 3-parameter families with a zero–Hopf equilibrium. As far as we know the number of periodic orbits which can bifurcate from a zero-Hopf equilibrium point when this is perturbed inside the class of all Lotka–Volterra systems (1) only has been studied partially in the paper [23] using averaging theory of second order. There the authors provided explicit conditions for the existence of one or two periodic orbits bifurcating from one of these equilibria. Here we shall use the averaging theory of third order for studying the number of periodic orbits which can bifurcate from a zero-Hopf equilibrium point of a Lotka–Volterra system (1). Previous results in this direction are the following. The bifurcation of periodic orbits in a Hopf equilibrium of a Lotka–Volterra system (1) have been studied by many authors. Thus in the papers [25–27] the authors proved that two periodic orbits can bifurcate from a Hopf equilibrium of system (1). While in [28–30] it is shown that three periodic orbits can bifurcate from a Hopf equilibrium. Recently in [31] it is proved that four periodic orbits can bifurcate from a Hopf equilibrium of system (1). All these previous results on the number of periodic orbits bifurcating from a Hopf equilibrium are when system (1) has all its coefficients aij and ri positive (i.e., for the so called competitive Lotka–Volterra systems), and under this assumption in [28] it is conjectured that at least five periodic orbits can bifurcate from a such Hopf equilibrium, but this conjecture remains open. Our main objective is to show that the previous conjecture cannot be extended to the general Lotka–Volterra systems (1), because we shall prove that there are non-competitive Lotka–Volterra systems having at least six periodic orbits bifurcating from a zero-Hopf equilibrium. In short, until now it is known that there are Lotka–Volterra systems (1) having at least four periodic orbits bifurcating from one of their equilibrium points. Our main result is the following one. Theorem 1. There are non-competitive Lotka–Volterra systems (1) having at least six periodic orbits bifurcating from a zero-Hopf equilibrium. The proof of Theorem 1 is given in Section 3. In Section 2 we describe the algorithm used for computing the periodic solutions bifurcating from a zero–Hopf equilibrium. 2. The Algorithm for Computing the Periodic Solutions Bifurcating from a Zero-Hopf Equilibrium Assume that we have a differential system ẋi = f i ( x1 , . . . , xn , λ) for i = 1, . . . , n, (2) defined in an open subset of Rn , and that the origin of coordinates is a zero–Hopf equilibrium for this system, i.e., the eigenvalues of the linear part of the system at the origin are ±ωi and 0 with multiplicity n − 2. Here λ = (λ1 , . . . , λm ) denotes the parameters of the system. (1) Since we want to apply the averaging theory of order three (see the Appendix A) for studying the periodic solutions bifurcating from the zero–Hopf equilibrium at the origin and the averaging theory uses a small parameter ε, we write the parameters of the system in the form λk = λk0 + ελk1 + ε2 λk2 + ε3 λk3 , for k = 1, . . . , m. (2) Due to the fact that the zero–Hopf bifurcation will take place in a neighborhood of the origin, where it is localized the zero–Hopf equilibrium, we blow up this neighborhood doing the scaling of variables ( x1 , . . . , xn ) = (εX1 , . . . , εXn ), Mathematics 2020, 8, 1137 3 of 14 using again the small parameter ε, and we obtain a differential system of the form Ẋi = Fi ( X1 , . . . , Xn , λ, ε) for i = 1, . . . , n. (3) (3) In order to simplify the future computations and also for applying the averaging theory described in the Appendix A we need that the right hand part of the differential system starts with order ε, for these two reasons we shall pass the linear part of the differential system (3) to its real Jordan normal form doing a convenient linear change of variables ( X1 , . . . , Xn ) → (u1 , . . . , un ). Thus the differential system (3) in the new variables (u1 , . . . , un ) writes u̇1 = −ωu2 + εg1 (u1 , . . . , un , λ, ε), u̇2 = ωu1 + εg2 (u1 , . . . , un , λ, ε), u̇3 = εg3 (u1 , . . . , un , λ, ε), (4) ··· ··· u̇n = εgn (u1 , . . . , un , λ, ε). (4) In order to apply the averaging theory to a differential system the right hand part of that differential system must be a periodic function in the independent variable of the system, see again the Appendix A. For this reason we first pass the differential system (4) to the generalized cylindrical coordinates (r, θ, u3 , . . . , un ) where u1 = r cos θ and u2 = r sin θ, and system (4) becomes ṙ = εG1 (r, θ, u3 , . . . , un , λ, ε), θ̇ = ω + εG2 (r, θ, u3 , . . . , un , λ, ε), u̇3 = εG3 (r, θ, u3 , . . . , un , λ, ε), (5) ··· ··· u̇n = εGn (r, θ, u3 , . . . , un , λ, ε). Now this differential system has its right hand part periodic in the variable θ, because this variable appears only through the functions cos θ and sin θ. Since the cylindrical coordinates are not well defined at r = 0, we are studying only the periodic orbits which does not intersect the set r = 0. (5) Now we take the variable θ as the new independent variable, and system (5) in this new independent variable writes r0 = εH1 (r, θ, u3 , . . . , un , λ, ε), u30 = εH3 (r, θ, u3 , . . . , un , λ, ε), (6) ··· ··· u0n = εHn (r, θ, u3 , . . . , un , λ, ε), where the prime denotes derivative with respect to the variable θ. Note that the differential system (6) is already written in the normal form (A1) for applying to it the averaging theory described in the Appendix A. We also must take care to look for the periodic solutions into the region where θ̇ does not vanish. (6) We apply the averaging theory of third order and according with it we may get s periodic solutions (r k (θ, ε), u3k (θ, ε), . . . , ukn (θ, ε)) of system (6) for k = 1, . . . , s such that (r k (0, ε), u3k (0, ε), . . . , ukn (0, ε)) = (r k∗ , u3k∗ , . . . , ukn∗ ) + O(ε), (7) Mathematics 2020, 8, 1137 4 of 14 where the values (r k∗ , u3k∗ , . . . , ukn∗ ) are the zeros of the first averaging function which is not identically zero, see again the Appendix A. (7) Now we go back through the changes of variables until the initial differential system (2), and we shall see how look the periodic solutions (r k (θ, ε), u3k (θ, ε), . . ., ukn (θ, ε)) in the initial differential system. First the periodic solutions (7) in the differential system (5) are (r k (t, ε), θ k (t, ε), u3k (t, ε), . . . , ukn (t, ε)) verifying (r k (0, ε), θ k (0, ε), u3k (0, ε), . . . , ukn (t, ε)) = (r k∗ , 0, u3k∗ , . . . , ukn∗ ) + O(ε), note that θ k (t, ε) = ωt + O(ε). Now the periodic solution (r k (t, ε), θ k (t, ε), u3k (t, ε), . . . , ukn (t, ε)) of system (5) in system (4) becomes (u1k (t, ε), u2k (t, ε), u3k (t, ε), . . . , ukn (t, ε)) satisfying (u1k (0, ε), u2k (0, ε), u3k (0, ε), . . . , ukn (0, ε)) = (r k∗ , 0, u3k∗ , . . . , ukn∗ ) + O(ε), because u1 = r cos θ and u2 = r sin θ. Furthermore, the periodic solution (u1k (t, ε), . . . , ukn (t, ε)) of system (4) in system (3) writes ( X1k (t, ε), . . . , Xnk (t, ε)) satisfying ( X1k (0, ε), . . . , Xnk (0, ε)) = ( X1k∗ , . . . , Xnk∗ ) + O(ε), here the values of X jk∗ for k = 1, . . . , n depend on the linear change of variables that pass the linear part of system (3) at the origin into its real Jordan normal form. Finally we pass the periodic solutions ( X1k (t, ε), . . . , Xnk (t, ε)) of system (3) to system (2) and we get the periodic solutions ( x1k (t, ε), . . . , xnk (t, ε)) satisfying ( x1k (0, ε), . . . , xnk (0, ε)) = (εX1k∗ , . . . , εXnk∗ ) + O(ε2 ). So all these periodic solutions when ε → 0 tend to the origin of coordinates, consequently are periodic solutions bifurcating from the zero-Hopf equilibrium localized at the origin. Here we have described this algorithm for the averaging theory up to third order, but the algorithm is the same for the averaging theory of an arbitrary order. See the averaging theory at any order in the papers [32–34]. 3. Proof of Theorem 1 If system (1) has a zero-Hopf equilibrium ( a, b, c) with non-zero components without loss of generality we can consider this equilibrium at the point (1, 1, 1) doing the rescaling ( x, y, z) → ( x/a, y/b, z/c). Then every Lotka–Volterra system (1) having the equilibrium (1, 1, 1) can be written as  ẋ = x a11 ( x − 1) + a12 (y − 1) + a13 (z − 1) ,  ẏ = y a21 ( x − 1) + a22 (y − 1) + a23 (z − 1) , (8)  ż = z a31 ( x − 1) + a32 (y − 1) + a33 (z − 1) , where now we denote the coordinates of R3 by ( x, y, z). Since we shall use the averaging theory of third order for studying the periodic orbits of this system we take the coefficients aij as follows aij = aij0 + εaij1 + ε2 aij2 + ε3 aij3 , with i and j varying in {1, 2, 3}, being ε a small parameter. Note that in the differential system (8) there are 37 parameters. This big number of parameters produce that the computations for studying the number of periodic orbits which can bifurcate from the equilibrium (1, 1, 1) are Mathematics 2020, 8, 1137 5 of 14 tedious and huge. All the computations of this paper has been done with the help of the algebraic manipulator mathematica. First we translate the equilibrium (1, 1, 1) to the origin of coordinates and system (8) becomes ẋ = (1 + x ) a110 x + a120 y + a130 z + ε( a111 x + a121 y + a131 z)+ ε2 ( a112 x + a122 y + a132 z) + ε3 ( a113 x + a123 y + a133 z) ,  ẏ = (1 + y) a210 x + a220 y + a230 z + ε( a211 x + a221 y + a231 z)+ (9) ε2 ( a212 x + a222 y + a232 z) + ε3 ( a213 x + a223 y + a233 z) ,  ż = (1 + z) a310 x + a320 y + a330 z + ε( a311 x + a321 y + a331 z)+ ε2 ( a312 x + a322 y + a332 z) + ε3 ( a313 x + a323 y + a333 z) .  In order that the linear part of system (9) at the origin has eigenvalues 0 and ±ωi with ω 6= 0 we choose the conditions a110 = a120 = a130 = a210 = 0, a320 = −( a2220 + ω 2 )/a230 and a330 = − a220 , (10) with a230 ω 6= 0. So the origin of system (9) is a zero-Hopf equilibrium, and consequently system (8) has a zero-Hopf equilibrium at the point (1, 1, 1). We remark that there are other conditions on the parameters which also provide that the point (1, 1, 1) be a zero-Hopf equilibrium. In what follows we shall study the periodic orbits bifurcating from the zero-Hopf equilibrium (0, 0, 0) of system (9) under conditions (10). Note that according with conditions (10) not all coefficients of the Lotka–Volterra system are positive, so we are working with non-competitive Lotka–Volterra systems. As we shall see the amount of computations for studying this Hopf-bifurcation following our algorithm are huge due to the big number of parameters in system (9). In order to study the periodic orbits bifurcating from the zero-Hopf equilibrium at the origin of the differential system (9) using the averaging theory of third order (see the Appendix A), we need to introduce a small parameter and take a new independent variable in which the differential system be periodic. The small parameter for the averaging theory will be the parameter ε, and we do the rescaling ( x, y, z) = (εX, εY, εZ ). Then system (9) in the new variables ( X, Y, Z ) writes Ẋ = ε( a111 X + a121 Y + a131 Z ) + ε2 ( a112 X + a111 X 2 + a122 Y + a121 XY + a132 Z + a131 XZ ) + ε3 ( a113 X + a112 X 2 + a123 Y + a122 XY + a133 Z + a132 XZ ) + O(ε4 ), Ẏ = a220 Y + a230 Z + ε( a211 X + a221 Y + a220 Y 2 + a231 Z + a230 YZ )+ ε2 ( a212 X + a222 Y + a211 XY + a221 Y 2 + a232 Z + a231 YZ )+ ε3 ( a213 X + a223 Y + a212 XY + a222 Y 2 + a233 Z + a232 YZ ) + O(ε4 ), (11) Ż = ( a230 a310 X − a2220 Y − a220 a230 Z − Yω 2 )/a230 + ε( a230 a311 X + a230 a321 Y + a230 a331 Z + a230 a310 XZ − a2220 YZ − a220 a230 Z2 − YZω 2 )/a230 + ε2 ( a312 X + a322 Y + a332 Z + a311 XZ + a321 YZ + a331 Z2 ) + ε3 ( a313 X + a323 Y + a333 Z + a312 XZ + a322 YZ + a332 Z2 ) + O(ε4 ). Mathematics 2020, 8, 1137 6 of 14 In order to simplify the computations of the averaging theory we shall write the linear part of the differential system (11) into its real Jordan normal form doing the linear change of variables ( X, Y, Z ) → (u, v, w) given by X= w, a230 a310 w a ωv − a220 a230 u Y= + 230 2 , ω2 a220 + ω 2 Z= − a220 a230 a310 w + a230 ω 2 u. Now the differential system (11) in the new variables (u, v, w) becomes ε u̇ = −ωv + 2 ( a131 a3220 a310 − a121 a2220 a230 a310 + 4 ω ( a220 + ω 2 ) a131 a220 a310 ω 2 − a220 a230 a321 ω 2 + a2220 a331 ω 2 + a331 ω 4 )ω 2 u+ a230 ( a121 a220 a310 + a321 ω 2 )ω 3 v − ( a2220 + ω 2 ) ( a131 a2220 a2310 − a121 a220 a230 a2310 − a111 a220 a310 ω 2 − a230 a310 a321 ω 2 + a220 a310 a331 ω 2 − a311 ω 4 )w − ω 5 uv + a220 a310 ω 3 vw + O(ε2 ),  ε v̇ = ωu + (− a3220 a221 a230 + a4220 a231 − a230 ω ( a2220 + ω 2 ) 3 a131 a2220 a230 a310 + a121 a220 a2230 a310 − a2220 a2230 a321 + a3220 a230 a331 − a220 a221 a230 ω 2 + 2a2220 a231 ω 2 − a131 a230 a310 ω 2 + a220 a230 a331 ω 2 + a231 ω 4 )ω 2 u + a230 ( a2220 a221 − a121 a230 a310 + a220 a230 a321 + a221 ω 2 )ω 3 v − ( a2220 + ω 2 )(− a2220 a221 a230 a310 + a3220 a231 a310 − a131 a220 a230 a2310 + a121 a2230 a2310 − a220 a2230 a310 a321 + a2220 a230 a310 a331 − a211 a2220 ω 2 + a111 a230 a310 ω 2 − a221 a230 a310 ω 2 + a220 a231 a310 ω 2 − (12) a220 a230 a311 ω 2 − a211 ω 4 )w − a220 a2230 ω 4 u2 − a230 ( a3220 + a2220 a230 + a220 ω 2 − a230 ω 2 )ω 3 uv + a2230 a310 ( a2220 + ω 2 )ω 2 uw+ a220 a2230 v2 ω 4 + a220 a230 ( a220 + a230 ) a310 ( a2220 + ω 2 )ωvw + O(ε2 ),  ε ẇ = 2 ( a131 a2220 − a121 a220 a230 + a131 ω 2 )ω 2 u+ ω ( a220 + ω 2 ) 2 a121 a230 ω 3 v − ( a2220 + ω 2 )( a131 a220 a310 − a121 a230 a310 − a111 ω 2 )w  +O(ε2 ). In the computations of the previous differential system we have obtained the expressions of u̇, v̇ and ẇ until terms of O(ε4 ), but here we only present them until terms of order O(ε2 ), otherwise the expression of system (12) would need several pages. Using an algebraic manipulator as mathematica or mapple it is relatively easy to repeat our computations. Now we write the differential system (12) in cylindrical coordinates (r, θ, w) where u = r cos θ and v = r sin θ, and taking θ as the new independent variable of the differential system defined we get the new differential system r0 = εF11 (θ, r, w) + ε2 F21 (θ, r, w) + ε3 F31 (θ, r, w) + O(ε4 ), (13) w0 = εF12 (θ, r, w) + ε2 F22 (θ, r, w) + ε3 F32 (θ, r, w) + O(ε4 ), Mathematics 2020, 8, 1137 7 of 14 defined in in r > 0, where the prime denotes derivative with respect to the variable θ. Here we only provide the explicit expressions of F11 = F11 (θ, r, w) and F12 = F12 (θ, r, w) which are the shorter ones, but our next computations will use the expressions of F21 , F22 , F31 and F32 . Thus we have 1 F11 = ( a230 ( a131 a3220 a310 − a121 a2220 a230 a310 + a230 ω ( a2220 + ω 2 ) 5 a131 a220 a310 ω 2 − a220 a230 a321 ω 2 + a2220 a331 ω 2 + a331 ω 4 ) cos2 θ + ( a4220 a231 − a3220 a221 a230 − a131 a2220 a230 a310 + 2a121 a220 a2230 a310 − a2220 a2230 a321 + a3220 a230 a331 − a220 a221 a230 ω 2 + 2a2220 a231 ω 2 − a131 a230 a310 ω 2 + a2230 a321 ω 2 + a220 a230 a331 ω 2 + a231 ω 4 )ω cos θ sin θ + a230 ( a2220 a221 − a121 a230 a310 + a220 a230 a321 + a221 ω 2 )ω 2 sin2 θ )ω 2 r − ( a2220 + ω 2 )( a230 ( a131 a2220 a2310 − a121 a220 a230 a2310 − a111 a220 a310 ω 2 − a230 a310 a321 ω 2 + a220 a310 a331 ω 2 − a311 ω 4 ) cos θ + ( a3220 a231 a310 − a2220 a221 a230 a310 − a131 a220 a230 a2310 + a121 a2230 a2310 − a220 a2230 a310 a321 + a2220 a230 a310 a331 − a211 a2220 ω 2 + a111 a230 a310 ω 2 − a221 a230 a310 ω 2 + a220 a231 a310 ω 2 − a220 a230 a311 ω 2 − a211 ω 4 )ω sin θ ])w − ( a230 ω ( a2220 + a220 a230 + ω 2 ) cos2 θ sin θ + a230 ( a3220 + a2220 a230 + a220 ω 2 − a230 ω 2 ) cos θ sin2 θ − a220 a2230 ω sin3 θ )ω 4 r2 + ( a230 ( a220 + a230 ) a310 ω 3 ( a2220 + ω 2 ) cos θ sin θ + a220 a230 ( a220 + a230 ) a310 ω 2 ( a2220 + ω 2 ) sin2 θ )rw ,  1 F12 = (( a131 a2220 − a121 a220 a230 + a131 ω 2 ) cos θ + ω 3 ( a2220 + ω 2 ) a121 a230 ω sin θ )ω 2 r − ( a2220 + ω 2 )( a131 a220 a310 − a121 a230 a310 − a111 ω 2 )w ,  We note that the differential system (13) is written in the normal form (A1) for applying the averaging theory of third order described in the Appendix A, where the variables t and x of the Appendix A are now θ and (r, w) respectively. Computing the averaged function of first order f 1 (r, w) = ( f 11 (r, w), f 12 (r, w)) defined in the Appendix A we get f 11 (r, w) = Ar, f 12 (r, w) = Bw, where ( a131 a220 − a121 a230 ) a310 + ( a221 + a331 )ω 2 + ( a220 + a230 ) a310 a220 w A= , 2ω 3 ( a a230 − a131 a220 ) a310 + a111 ω 2 B = 121 . ω3 We look for the zeros (r ∗ , w∗ ) of f 1 (r, w) with r > 0, and since the unique zero of the function f 1 (r, w) is (0, 0), or a continuum of zeros if the coefficient A or B is zero, the averaged function of first order does not give any information on the periodic solutions of system (13), see the Appendix A. Therefore we force that the averaged function of first order be identically zero and we shall use the averaged functions of higher order to obtain information on the periodic solutions of the differential system (13). Since the coefficient of rw in the function f 11 (r, w) is ( a220 + a230 ) a310 a220 we need to consider the following three cases in order that the averaged function of first order be identically zero: Case 1: a220 = − a230 , a331 = ( a121 a230 a310 + a131 a230 a310 − a221 ω 2 )/ω 2 , a111 = (− a121 a230 a310 − a131 a230 a310 )/ω 2 . Case 2: a310 = 0, a331 = − a221 , a111 = 0. Mathematics 2020, 8, 1137 8 of 14 Case 3: a220 = 0, a331 = ( a121 a230 a310 − a221 ω 2 )/ω 2 , a111 = −( a121 a230 a310 )/ω 2 . Case 1. Since the averaged function of first order f 1 (r, w) is identically zero, we compute the averaged function of second order f 2 (r, w) = ( f 21 (r, w), f 22 (r, w)) and we obtain f 21 (r, w) = (Cw + D )r, f 22 (r, w) = Ew, where C= −(− a121 a2230 a2310 − a131 a2230 a2310 + a121 a230 a310 ω 2 + a131 a230 a310 ω 2 + a221 a230 a310 ω 2 + a230 a231 a310 ω 2 + a211 ω 4 )/(2ω ), D= −( a2121 a4230 a2310 + 2a121 a131 a4230 a2310 + a2131 a4230 a2310 − 2a121 a221 a3230 a310 ω 2 −2a131 a221 a3230 a310 ω 2 − a121 a3230 a231 a310 ω 2 − a131 a3230 a231 a310 ω 2 + a121 a3230 a310 a321 ω 2 + a131 a3230 a310 a321 ω 2 − a121 a211 a2230 ω 4 − a131 a211 a2230 ω 4 − a131 a221 a230 a310 ω 4 + a122 a2230 a310 ω 4 + a132 a2230 a310 ω 4 − a131 a230 a231 a310 ω 4 + a121 a2230 a311 ω 4 + a131 a2230 a311 ω 4 − a131 a211 ω 6 − a222 a230 ω 6 − a230 a332 ω 6 )/(2a230 ω 7 ), E= ( a2121 a4230 a2310 + 2a121 a131 a4230 a2310 + a2131 a4230 a2310 − 2a121 a221 a3230 a310 ω 2 −2a131 a221 a3230 a310 ω 2 − a121 a3230 a231 a310 ω 2 − a131 a3230 a231 a310 ω 2 + a121 a3230 a310 a321 ω 2 + a131 a3230 a310 a321 ω 2 − a121 a211 a2230 ω 4 − a131 a211 a2230 ω 4 − a131 a221 a230 a310 ω 4 + a122 a2230 a310 ω 4 + a132 a2230 a310 ω 4 − a131 a230 a231 a310 ω 4 + a121 a2230 a311 ω 4 + a131 a2230 a311 ω 4 − a131 a211 ω 6 + a112 a230 ω 6 )/( a230 ω 7 ). Again the unique zero of the averaged function of second order f 2 (r, w) is the (0, 0) or a continuum of solutions in case that convenient coefficients C, D or E are zero. Therefore the averaging theory of second order does not provide any information on the periodic solutions of the differential system (13). Consequently we impose that the averaged function of second order f 2 (r, w) be identically zero, and we obtain that a211 = a121 a2230 a2310 + a131 a2230 a2310 − a121 a230 a310 ω 2 − a131 a230 a310 ω 2 − a221 a230 a310 ω 2 − a230 a231 a310 ω 2 /ω 4 ,  a332 = a2121 a2230 a310 + 2a121 a131 a2230 a310 + a2131 a2230 a310 − a121 a221 a2230 a310 − a131 a221 a2230 a310 − a121 a131 a230 a2310 − a2131 a230 a2310 + a121 a2230 a310 a321 + a131 a2230 a310 a321 + a121 a131 a310 ω 2 + a2131 a310 ω 2 + a122 a230 a310 ω 2 + a132 a230 a310 ω 2 + a121 a230 a311 ω 2 + a131 a230 a311 ω 2 − a222 ω 4 /ω 4 ,  a112 = − a2121 a2230 a310 − 2a121 a131 a2230 a310 − a2131 a2230 a310 + a121 a221 a2230 a310 + a131 a221 a2230 a310 + a121 a131 a230 a2310 + a2131 a230 a2310 − a121 a2230 a310 a321 − a131 a2230 a310 a321 − a121 a131 a310 ω 2 − a2131 a310 ω 2 − a122 a230 a310 ω 2 − a132 a230 a310 ω 2 − a121 a230 a311 ω 2 − a131 a230 a311 ω 2 /ω 4 .  We compute the averaged function of third order f 3 (r, w) = ( f 31 (r, w), f 32 (r, w)) and we get a0 r4 + a1 r3 + a2 r2 w + a3 r2 + a4 rw + a5 w2 + a6 r + a7 w f 31 (r, w) = , 384a230 ( a2230 + ω 2 )ω 13 r b0 r3 + b1 r2 w + b2 r2 + b3 rw + b4 r + b5 w f 32 (r, w) = . 24a230 ( a2230 + ω 2 )2 ω 9 We do not provide the explicit expressions of the coefficients a j and b j because we shall need approximately twenty pages for writing them. Mathematics 2020, 8, 1137 9 of 14 Now we shall study the zeros of the function f 3 (r, w). Since the variable w appears linearly in the equation f 32 (r, w) = 0, we isolate it and we get w = W (r ). Substituting w = W (r ) into the equation f 31 (r, w) = 0, we obtain an equation in the variable r of the form n (r ) c r 2 + c3 r 3 + c4 r 4 + c5 r 5 + c6 r 6 + c7 r 7 + c8 r 8 = 2 = 0. (14) d (r ) ( d0 + d1 r + d2 r 2 )2 The coefficients c j and d j are polynomials in some of the coefficients of the differential system (8), more precisely in the coefficients a113 , a121 , a122 , a123 , a131 , a132 , a133 , a212 , a221 , a222 , a223 , a230 , a231 , a232 , a310 , a311 , a312 , a321 , a322 , a333 , ω. We have computed the rank of the Jacobian matrix of the function (c2 , c3 , c4 , c5 , c6 , c7 , c8 ) with respect to the 21 previous coefficients, it is the rank of a 7 × 23 matrix, and we get that this rank is 7. Therefore the seven coefficients of the polynomial n(r ) are independent, and consequently we can choose them in such a way that the polynomial n(r ) has six positive real roots. Moreover, we also can choose those coefficients in such a way that the resultant of the polynomials n(r ) and d(r ) is not zero, and consequently both polynomials do not have a common root. So Equation (14) can have six positive solutions, r ∗j for j = 1, 2, 3, 4, 5, 6. In short, we have that (r ∗j , W (r ∗j )) for j = 1, 2, 3, 4, 5, 6 are six zeros of the third averaged function f 3 (r, w). These zeros can be chosen simple, i.e., the Jacobian of the function f 3 (r, w) evaluated in such zeros is not zero. Consequently by the averaging theory (see the Appendix A) the differential system (13) has six periodic solutions (r j (θ, ε), w j (θ, ε)) such that (r j (0, ε), w j (0, ε)) → (r ∗j , W (r ∗j )) when ε → 0. Going back to the differential system (12) we obtain for this system six periodic solutions (u j (t, ε), v j (t, ε), w j (t, ε)) such that (u j (0, ε), v j (0, ε), w j (0, ε)) → (r ∗j , 0, W (r ∗j )), when ε → 0. These periodic solutions provide six periodic solutions ( X j (t, ε), Yj (t, ε), Zj (t, ε)) for the differential system (11) such that X j (0, ε) → W (r ∗j ), a230 a310 W (r ∗j ) a220 a230 r ∗j Yj (0, ε) → − , ω2 a2220 + ω 2 Zj (0, ε) → a230 ω 2 r ∗j − a220 a230 a310 W (r ∗j ), when ε → 0. Finally going back to the differential system (8) we obtain six periodic solutions ( x j (t, ε), y j (t, ε), z j (t, ε)) such that x j (0, ε) = 1 + εW (r ∗j ) + O(ε2 ), a230 a310 W (r ∗j a220 a230 r ∗j  y j (0, ε) = 1 + ε − + O ( ε2 ), (15) ω2 a2220 + ω 2 z j (0, ε) = 1 + ε a230 ω 2 r ∗j − a220 a230 a310 W (r ∗j ) + O(ε2 ),  when ε → 0. Clearly from (15) these six periodic solutions ( x j (t, ε), y j (t, ε), z j (t, ε)) tend to the equilibrium point (1, 1, 1) of the differential system (8) when ε → 0. Hence they bifurcate from that zero-Hopf equilibrium at ε = 0. This completes the proof of Theorem 1. Case 2. Again, since the averaged function of first order f 1 (r, w) is identically zero, we compute the averaged function of second order f 2 (r, w) = ( f 21 (r, w), f 22 (r, w)) and we obtain f 21 (r, w) = (Cw + D )r, f 22 (r, w) = Ew, Mathematics 2020, 8, 1137 10 of 14 where a220 ( a211 a2220 + a211 a220 a230 + a220 a230 a311 + a2230 a311 + a211 ω 2 ) C= , 2a230 ω 3 1 D= a121 a211 a220 a230 − a131 a211 a2220 − a131 a220 a230 a311 + 2a230 ω 3 a121 a2230 a311 − a131 a211 ω 2 − a222 a230 ω 2 − a230 a332 ω 2 ,  1 E= a121 a211 a220 a230 − a131 a211 a2220 − a131 a220 a230 a311 + a230 ω 3 a121 a2230 a311 − a131 a211 ω 2 + a112 a230 ω 2 .  As in the previous case the unique zero of the averaged function of second order f 2 (r, w) is the (0, 0) or a continuum of solutions in case that convenient coefficients C, D or E are zero. Consequently we impose that the averaged function of second order f 2 (r, w) be identically zero, but since the coefficient of rw in the function f 21 (r, w) is a product of two factors we have to consider two subcases. Subcase 2.1: a220 = 0. Then, in order that the averaged function of second order f 2 (r, w) be identically zero we take a121 a2230 a311 − a131 a211 ω 2 − a222 a230 ω 2 a332 = , a230 ω 2 a131 a211 ω 2 − a121 a2230 a311 a112 = . a230 ω 2 We compute the averaged function of third order f 3 (r, w) = ( f 31 (r, w), f 32 (r, w)) and we get a0 r3 + a1 r2 + a2 rw + a3 w2 + a4 r + a5 w f 31 (r, w) = , 384a2230 ω 5 r (16) b0 r3 + b1 r2 + b2 rw + b3 r + b4 w f 32 (r, w) = . 24a2230 ω 5 Here the expressions of the coefficients a j ’s and b j ’s are relatively short, but we do not need them explicitly. We shall study the zeros of the function f 3 (r, w). Since the variable w appears linearly in the equation f 32 (r, w) = 0, we isolate it and we get w = W (r ). Substituting w = W (r ) into the equation f 31 (r, w) = 0, we obtain an equation in the variable r of the form c2 r 2 + c3 r 3 + c4 r 4 + c5 r 5 + c6 r 6 = 0. ( d0 + d1 r + d2 r 2 )2 So at most we have four positive solutions for the variable r, and consequently at most four zeros for the averaged function of third order f 3 (r, w). In any case less than the six obtained in Case 1. Subcase 2.2: a211 a2220 + a211 a220 a230 + a220 a230 a311 + a2230 a311 + a211 ω 2 = 0. Then in order that the averaged function of second order f 2 (r, w) be identically zero we take a211 a2220 + a211 a220 a230 + a211 ω 2 a311 = − , a230 ( a220 + a230 ) a a + a131 a211 + a220 a222 + a222 a230 a332 = − 121 211 , a220 + a230 a121 a211 + a131 a211 a112 = . a220 + a230 We compute the averaged function of third order f 3 (r, w) = ( f 31 (r, w), f 32 (r, w)) and we get again the expressions given in (16), of course the coefficients a j ’s and b j ’s are now different. Repeating the arguments of the previous subcase we obtain at most four zeros for the averaged function of third order f 3 (r, w). Mathematics 2020, 8, 1137 11 of 14 Case 3. Again since the averaged function of first order f 1 (r, w) is identically zero, we compute the averaged function of second order f 2 (r, w) = ( f 21 (r, w), f 22 (r, w)) and we obtain f 21 (r, w) = (Cw + D )r, f 22 (r, w) = Ew, where a310 ( a121 − a221 ) a230 C= − , 2 ω3 1 D= − a121 a3230 a310 a321 − a131 a221 a230 a310 ω 2 + a122 a2230 a310 ω 2 + 2a230 ω 5 a121 a2230 a311 ω 2 − a131 a211 ω 4 − a222 a230 ω 4 − a230 a332 ω 4 ,  1 E= a121 a3230 a310 a321 − a131 a221 a230 a310 ω 2 + a122 a2230 a310 ω 2 + a230 ω 5 a121 a2230 a311 ω 2 − a131 a211 ω 4 + a112 a230 ω 4 .  As in the previous case the unique zero of the averaged function of second order f 2 (r, w) is the (0, 0) or a continuum of solutions in case that convenient coefficients C, D or E are zero. Consequently we impose that the averaged function of second order f 2 (r, w) be identically zero, but since the coefficient of rw in the function f 21 (r, w) is a product of two factors which can be zero, namely a310 ( a121 − a221 ), we have two consider two subcases. Subcase 3.1: a310 = 0. Then in order that the averaged function of second order f 2 (r, w) be identically zero we take a121 a2230 a311 − a131 a211 ω 2 − a222 a230 ω 2 a332 = , a230 ω 2 − a121 a2230 a311 + a131 a211 ω 2 a112 = . a230 ω 2 We compute the averaged function of third order f 3 (r, w) = ( f 31 (r, w), f 32 (r, w)) and we get again the expression given in (16), consequently at most four solutions. Subcase 3.2: a221 = a121 . Then in order that the averaged function of second order f 2 (r, w) be identically zero we take 1 a332 = a121 a3230 a310 a321 − a121 a131 a230 a310 ω 2 + a122 a2230 a310 ω 2 + a230 ω 2 a121 a2230 a311 ω 2 − a131 a211 ω 4 − a222 a230 ω 4 ,  1 a112 = a121 a131 a230 a310 ω 2 − a121 a3230 a310 a321 − a122 a2230 a310 ω 2 − a230 ω 4 a121 a2230 a311 ω 2 + a131 a211 ω 4 .  We compute the averaged function of third order f 3 (r, w) = ( f 31 (r, w), f 32 (r, w)) and we get a0 r4 w + a1 r4 + a2 r2 w2 + a3 r3 + a4 r2 w + a5 w3 + a6 r2 + a7 rw + a8 w2 f 31 (r, w) = , 384a2230 ω 9 r b0 r3 + b1 r2 w + b2 r2 + b3 rw + b4 w2 + b5 r + b6 w f 32 (r, w) = . 24a2230 ω 7 Here the explicit expressions of the coefficients a j ’s and b j ’s only should need approximately three pages for writing them. However, unfortunately in this case we do not know how to control the zeros (r ∗ , w∗ ) of the function f 3 (r, w) with r ∗ > 0. We think that in this subcase it is possible that more than six simple zeros can be obtained, but for the moment this is an open problem. Mathematics 2020, 8, 1137 12 of 14 Author Contributions: All the authors have contributed in a similar way to this paper in methodology, software, validation, formal analysis, investigation, resources, data curation, writing original draft preparation, writing review and editing, visualization, and supervision. All authors have read and agreed to the published version of the manuscript. Funding: The first and the third authors are partially supported by National Natural Sciences Foundation 11431008, 11771296 and 11501370. The second author is partially supported by the Ministerio de Ciencia, Innovación y Universidades, Agencia Estatal de Investigación grants MTM2016-77278-P (FEDER) and MDM-2014-0445, the Agència de Gestió d’Ajuts Universitaris i de Recerca grant 2017SGR1617, and the H2020 European Research Council grant MSCA-RISE-2017-777911. Acknowledgments: We thank to the reviewers their comments and suggestions. Conflicts of Interest: The authors declare no conflict of interest. Appendix A. The Averaging Theory of First, Second and Third Order The averaging theory of third order for studying periodic orbits was developed [35] and in [34] at any order. It can be summarized as follows. Consider the differential system ẋ = εF1 (t, x ) + ε2 F2 (t, x ) + ε3 F3 (t, x ) + ε4 R(t, x, ε), (A1) where F1 , F2 , F3 : R × D → R, R : R × D × (−ε f , ε f ) → R are continuous functions, T–periodic in the first variable, and D is an open subset of Rn . Assume that the following hypotheses (i) and (ii) hold. (i) F1 (t, ·) ∈ C2 ( D ), F2 (t, ·) ∈ C1 ( D ) for all t ∈ R, F1 , F2 , F3 , R, Dx2 F1 , Dx F2 are locally Lipschitz with respect to x, and R is twice differentiable with respect to ε. We define Fk0 : D → R for k = 1, 2, 3 as 1 T Z f1 (x) = F (s, x )ds, T 0 1 1 T Z f2 (x) = [ Dx F1 (s, x ) · y1 (s, x ) + F2 (s, x )] ds, T 0 1 T h1 ∂2 F 1 ∂F1 Z f3 (x) = y1 (s, x ) T 21 (s, x )y1 (s, x ) + (s, x )y2 (s, x ) T 0 2 ∂x 2 ∂x ∂F i + 2 (s, x )y1 (s, x ) + F3 (s, x ) ds, ∂x where Z s y1 (s, x ) = F1 (t, x )dt, Z0 s  Z t  ∂F1 y2 (s, x ) = (t, x ) F1 (r, x )dr + F2 (t, x ) dt. 0 ∂x 0 (ii) For an open and bounded set V ⊂ D and for each ε ∈ (−ε f , ε f ) \ {0}, there exists a ∈ V such that f 1 ( a) + ε f 2 ( a) + ε2 f 3 ( a) = 0 and d B ( f 1 + ε f 2 + ε2 f 3 , V, aε ) 6= 0 (i.e., the Brouwer degree of the function f 1 + ε f 2 + ε2 f 3 at the point a is not zero). Then for |ε| > 0 sufficiently small there exists a T–periodic solution x (t, ε) of system (A1) such that x (0, ε) → a when ε → 0. A sufficient condition in order that d B ( f 1 + ε f 2 + ε2 f 3 , V, aε ) 6= 0 is that the Jacobian of the function f 1 + ε f 2 + ε2 f 3 at a is not zero, see for details [36]. The averaging theory of first order takes place when f 1 is not identically zero. Therefore the zeros of f 1 + ε f 2 + ε2 f 3 are mainly the zeros of f 1 for ε sufficiently small. The averaging theory of second order takes place when f 1 is identically zero and f 2 is not identically zero. Then the zeros of f 1 + ε f 2 + ε2 f 3 are mainly the zeros of f 2 for ε sufficiently small. Mathematics 2020, 8, 1137 13 of 14 Finally the averaging theory of third order takes place when f 1 and f 2 are identically zero and f 3 is not identically zero. Therefore the zeros of f 1 + ε f 2 + ε2 f 3 are mainly the zeros of f 3 for ε sufficiently small. References 1. Baldomá, I.; Castejón, O.; Seara, M.T. Exponentially small heteroclinic breakdown in the generic Hopf-Zero singularity. J. Dyn. Diff. Equ. 2013, 25, 335–392. [CrossRef] 2. Broer, H.W.; Vegter, G. Subordinate Silnikov bifurcations near some singularities of vector fields having low codimension. Ergodic Theory Dyn. Syst. 1984, 4, 509–525. [CrossRef] 3. Champneys, A.R.; Kirk, V. The entwined wiggling of homoclinic curves emerging from saddle-node/Hopf instabilities. Physica D 2004, 195, 77–105. [CrossRef] 4. Dumortier, F.; Ibá nez, S.; Kokubu, H.; Simó, C. About the unfolding of a Hopf-zero singularity. Discrete Contin. Dyn. Syst. 2013, 33, 4435–4471. [CrossRef] 5. Guckenheimer, J. On a codimension two bifurcation. Lect. Notes Math. 1980, 898, 99–142. 6. Guckenheimer, J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields; Springer: Berlin/Heidelberg, Germany, 1983. 7. Kuznetsov, Y.A. Elements of Applied Bifurcation Theory, 3rd ed.; Springer: Berlin/Heidelberg, Germany, 2004. 8. Scheurle, J.; Marsden, J. Bifurcation to quasi-periodic tori in the interaction of steady state and Hopf bifurcations. SIAM J. Math. Anal. 1984, 15, 1055–1074. [CrossRef] 9. Hirsch, M.W. Systems of differential equations that are competitive or cooperative. I: Limit sets. SIAM J. Math. Anal. 1982, 13, 167–169. [CrossRef] 10. Hirsch, M.W. Systems of differential equations that are competitive or cooperative. II: Convergence almost everywhere. SIAM J. Math. Anal. 1985, 16, 423–439. [CrossRef] 11. Hirsch, M.W. Stability and convergence in stronly monotone dynamical systems. J. Die Reine Angew. Mathematik 1988, 383, 1–53. 12. Hirsch, M.W. Systems of differential equations that are competitive or cooperative. III: Competing species. Nonlinearity 1988, 1, 117–124. [CrossRef] 13. Hirsch, M.W. Systems of differential equations that are competitive or cooperative. IV: Structural stability in 3-dimensional systems. SIAM J. Math. Anal. 1990, 21, 1225–1234. [CrossRef] 14. Hirsch, M.W. Systems of differential equations that are competitive or cooperative. V: Convergence in 3-dimensional systems. J. Differ. Equ. 1989, 80, 94–106. [CrossRef] 15. Hirsch, M.W. Systems of differential equations that are competitive or cooperative. VI: A local Cr closing lemma for 3-dimensional systems. Ergod. Theory Dyn. Syst. 1990, 11, 443–454. [CrossRef] 16. Lotka, A.J. Analytical note on certain rhythmic relations in organic systems. Proc. Natl. Acad. Sci. USA 1920, 6, 410–415. [CrossRef] 17. May, R.M. Stability and Complexity in Model Ecosystems; Princeton University Press: Princeton, NJ, USA, 1974. 18. Smale, S. On the differential equations of species in competition. J. Math. Biol. 1976, 3, 5–7. [CrossRef] 19. Volterra, V. Lecons sur la Théorie Mathématique de la Lutte pour la vie; Gauthier Villars: Paris, France, 1931. 20. Zeeman, M.L. Hopf bifurcations in competitive three-dimensional Lotka–Volterra systems. Dyn. Stab. Sys. 1993, 8, 189–217. [CrossRef] 21. Laval, G.; Pellat, R. Plasma Physics. Proceedings of Summer School of Theoretical Physics; Gordon and Breach: New York, NY, USA, 1975. 22. Busse, F.H. Transition to Turbulence via the Statistical Limit Cycle Route Synergetics; Springer: Berlin/Heidelberg, Germany, 1978; p. 39. 23. Llibre, J.; Xiao, D. Limit cycles bifurcating from a non-isolated zero-Hopf equilibrium of 3-dimensional differential systems. Proc. Am. Math. Soc. 2014, 142, 2047–2062. [CrossRef] 24. Arnold, V.I. Arnold’s Problems; Springer: Berlin/Heidelberg, Germany, 2000. 25. Hofbauer, J.; So, J.W.H. Multiple limit cycles for three dimensional Lotka–Volterra equations. Appl. Math. Lett. 1994, 7, 65–70. [CrossRef] 26. Lu, Z.; Luo, Y. Two limit cycles in three-dimensional Lotka–Volterra systems. Comput. Math. Appl. 2002, 44, 51–66. [CrossRef] Mathematics 2020, 8, 1137 14 of 14 27. Xiao, D.; Li, W. Limit cycles for the competitive three dimensional Lotka–Volterra system. Differ. Equ. 2000, 164, 1–15. [CrossRef] 28. Gyllenberg, M.; Yan, P. On the number of limit cycles for the three-dimensional Lotka–Volterra systems. Discrete Contin. Dyn. Syst. Ser. S 2009, 11, 347–352. [CrossRef] 29. Gyllenberg, M.; Yan, P.; Wang, Y. A 3D competitive LotkaVolterra system with three limit cycles: A falsification of a conjecture by Hofbauer and So. Appl. Math. Lett. 2006, 19, 1–7. [CrossRef] 30. Lu, Z.; Luo, Y. Three limit cycles for a three-dimensional Lotka–Volterra competitive system with a heteroclinic cycle. Comput. Math. Appl. 2003, 46, 231–238. [CrossRef] 31. Yu, P.; Han, M.; Xiao, D. Four small limit cycles around a Hopf singular point in 3-dimensional competitive Lotka–Volterra systems. J. Math. Anal. Appl. 2016, 436, 521–555. [CrossRef] 32. Giné, J.; Grau, M.; Llibre, J. Averaging theory at any order for computing periodic orbits. Physica D 2013, 250, 58–65. [CrossRef] 33. Giné, J.; Llibre, J.; Wu, K.S.; Zhang, X. Averaging methods of arbitrary order, periodic solutions and integrability. J. Differ. Equ. 2016, 260, 4130–4156. [CrossRef] 34. Llibre, J.; Novaes, D.D.; Teixeira, M.A. Higher order averaging theory for finding periodic solutions via Brouwer degree. Nonlinearity 2014, 27, 563–583. [CrossRef] 35. Buica, A.; Llibre, J. Averaging methods for finding periodic orbits via Brouwer degree. Bull. Sci. Math. 2004, 128, 7–22. [CrossRef] 36. Lloyd, N.G. Degree Theory; Cambridge University Press: Cambridge, UK, 1978. c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).