Bifurcation of Limit Cycles and Center in 3D Cubic Systems with Z 3 -Equivariant Symmetry

: This paper focuses on investigating the bifurcation of limit cycles and centers within a speciﬁc class of three-dimensional cubic systems possessing Z 3 -equivariant symmetry. By calculating the singular point values of the systems, we obtain a necessary condition for a singular point to be a center. Subsequently, the Darboux integral method is employed to demonstrate that this condition is also sufﬁcient. Additionally, we demonstrate that the system can bifurcate 15 small amplitude limit cycles with a distribution pattern of 5 − 5 − 5 originating from the singular points after proper perturbation. This ﬁnding represents a novel contribution to the understanding of the number of limit cycles present in three-dimensional cubic systems with Z 3


Introduction
The study of limit cycles in differential systems comes from practice, and many applied disciplines, such as physics, mechanics, biology, economics, ecology, electronics, and finance, may exhibit limited cycle phenomena.Limit cycles play a very important role in the study of nonlinear system dynamics.The limit cycle phenomenon was first proposed by H. Poincaré [1], and further development was promoted by the famous Hilbert's 16th problem [2].The latter part of the sixteenth problem is as follows: what is the maximum number of limit cycles (Hilbert number: H(n)) that can be generated by a planar differential system of degree n?For n = 2, the best result was obtained by Shi [3] and Chen [4] more than forty years ago, H(2) ≥ 4. For n = 3, scholars have conducted a lot of research work, and the best result for cubic polynomial systems is H(3) ≥ 13 [5,6].For more information of n = 4, 5, 6, . .., see [7][8][9].
The limit cycles of a differential autonomous system are always perturbed by singular points, but it is difficult to obtain a higher number of limit cycles in terms of a single singular point.Considering this situation, it is a good direction to study from multiple singular points.In the case of two singular points, Du, Liu and Huang [10] studied a class of Kolmogorov system that is symmetrical about the line y = x, and showed that under certain conditions, each singularity can bifurcate into five limit cycles so there are ten limit cycles in total.Du, Huang and Zhang [11] proved that a class of systems symmetric with regard to the origin can generate 12 limit cycles of 6 − 6 distribution by Hopf bifurcation.For three singular points, Yu and Han [12] used the perturbation technique to calculate the focus values and proved that a class of the Z 2 -equivariant system has at most 12 limit cycles.For four singular points, Wang, Liu and Du [13] considered a class of Z 3 -equivariant systems and obtained the conclusion that H(4) ≥ 16.For more singular points, Yao and Yu [14] applied the normal form and the technique of solving coupled multivariate polynomial equations to prove that the maximal number of small limit cycles is 25 in Z 5 -equivariant planar vector fields, so H(5) ≥ 25.Li, Chan and Chung [15] gave a concrete numerical example of at least 24 limit cycles for a perturbed Hamiltonian system with Z 6 symmetry, and conjectured that H(2k + 1) ≥ (2k + 1) 2 − 1 for the perturbed Hamiltonian systems.Li and Zhang [16] used the detection function method to show that H (7) ≥ 49 by considering a Z 8 -equivariant vector field.Through analyzing a Z 10 symmetric system, Wang and Yu [17] obtained the result H(9) ≥ 80.Recently, the fractional-order dynamic model has been widely studied.It is considered to be more valuable than the integer-order dynamic model to describe many objective natural phenomena in the real world because of its good memory characteristics and genetic function in many materials and evolutionary processes.Huang et al. [18,19] and Xu et al. [20][21][22] focused on the understanding of the bifurcation phenomena and dynamical behavior in fractional-order neural networks with different delays.They provide insights into the impact of fractional order and multiple delays on the stability, equilibrium points, and bifurcation patterns in neural network dynamics.
Surprisingly, unlike the planar systems, where the number of limit cycles is finite [23], a simple three-dimensional system may have an infinite number of small amplitude limit cycles [24].At present Huang, Gu and Wang [25] introduced some results on the bifurcation of limit cycles for three-dimensional smooth systems, such as Lotka-Volterra systems, Chen systems, Lorenz systems, Lü systems, and three-dimensional piecewise smooth systems.Sanchez and Torregrosa [26] studied several classes of three-dimensional polynomial differential systems by using the parallel computing method, and determined that there are at most 11, 31, 54, and 92 limit cycles at the origin for three-dimensional quadratic, cubic, quartic, and quintic systems, respectively.For multiple singular points, Gu et al. [27] proved the existence of seven limit cycles in a class of three-dimensional cubic Kolmogorov systems.Lu et al. [28] proved the existence of eight limit cycles in a class of Lorenz systems, which are Z 2 symmetric and quadratic three-dimensional systems.Guo, Yu and Chen [29] applied the theory of center manifold and normal forms to prove the existence of 12 limit cycles with 4 − 4 − 4 distribution in three-dimensional vector fields with Z 3 symmetry.Then, they analyzed a 3D quadratic system with two symmetric singular points by using the same approach in [30] and showed that there are twelve limit cycles with 6 − 6 distributed around the singular points.
Based on the above research results, numerous issues remain to be investigated in the dynamics of three-dimensional differential systems.In this study, we analyze the limit cycle bifurcation of a class of Z 3 -equivariant symmetric cubic systems by using symmetry to find the focal points with the same topological structure, and the form of the systems is as follows : Regarding the investigations into Hopf bifurcation of ordinary differential models, several classical methods are available.These include the Poincaré normal form, Liapunov-Schmidt method [31] and the dimension reduction method of directly calculating the center manifold, see, for example, [32].Additionally, various other research approaches have been proposed, such as the averaging theory [33,34], the technique of the inverse Jacobi multiplier [35,36], the simple normal form method [37] and the formal first integral method [38].Recently, the authors of [39] presented a theorem on the bound of the cyclicity in terms of the Bautin ideal to determine the maximum number of limit cycles within the center manifold generated from a center.
In our research, we employ the linear and straightforward algorithm proposed in a previous work [40] to compute the quantities of singular points without the need for dimensionality reduction.This algorithm is linear in nature and avoids complex integration operations.As a result, the calculations can be easily performed using symbolic computation systems, such as Mathematica.Furthermore, the weight polynomial and characteristic curve methods introduced in [41] are utilized to search for invariant algebraic surfaces, namely, all the Darboux polynomials under certain conditions.Then, the center conditions are derived, and the maximum number of limit cycles via a Hopf bifurcation is determined.Finally, we demonstrate that the system (1) can bifurcate 15 small amplitude limit cycles with a 5 − 5 − 5 distribution from the singular points after proper perturbation.The result is a new lower bound on the number of limit cycles in three-dimensional cubic systems with Z 3 -equivariant symmetry.
This paper is structured as follows.Section 2 introduces fundamental theories and methods of differential Hopf bifurcation systems, which are essential for a comprehensive discussion of system (1).In Section 3, we present a recursive formula for singular point values, which is linear and can simplify complex integral operations.As a result, we can conveniently use computer symbolic operation systems, such as Mathematica or Maple, to calculate the singular point values of system (1).In Section 4, we establish a necessary condition for the singular point to be a center, and then we use the Darboux integral method to demonstrate that this condition is also sufficient.Finally, in Section 5, we show that system (1) can bifurcate 15 small amplitude limit cycles with a 5 − 5 − 5 distribution from the singular points after proper perturbation.

Preliminary Theory and Method
For the Z q -equivariant system, we only present the results we need in the next section in this section.For more details, the reader can refer to [42,43] Lemma 1 (see [43]).If a system is a Z q -equivariant system and it has a real singular point except for the origin, then the system must have q real singular points, which are q-symmetric about the origin.
To discuss the limit cycle of a system, it is usually necessary to calculate all the focal values at the singular point.The consequent difficulty arises that the calculation of all the focus values of the singular point of the system is too complicated for us to calculate all the focus values.However, according to Hilbert's finite basis principle, there must be a finite number of focus values that can express all the focus values of the singular point of the system (the finite number of focus values in front is called the focus basis).Wang [44] proved that the first non-zero focal values v 2m+1 at the origin of system (2) are algebraically equivalent to the first non-zero singular values µ m at the origin of its complex system (4) (and the Lyapunov constant V 2m at the origin of the center manifold of system (2)).Therefore, the study of the focal values of the differential system can be simplified to the investigation of its singular point values.The focus basis issue, namely the highest-order fine focus problem, can be effectively solved in this process, and then the maximum number of limit cycle branches is also analyzed.In this paper, we apply the formal series method proposed by Wang [40] to study system (1).
Next, we show (for details see [40]) some preparations for the calculation of the singular point values of the system.Consider the following three-dimensional system: A kjl x k y j u l = X(x, y, u), where x, y, u, t, d, A kjl , B kjl , d kjl ∈ R(k, j, l ∈ N) and d > 0.
The linear matrix of system (2) has a pair of purely imaginary eigenvalues ±i and a negative eigenvalue −d at the origin of the singular point.Reviewing the related concepts and results, it is obvious that Hopf bifurcation occurs in system (2).The qualitative analysis of the high-dimensional systems is generally utilized to reduce their dimensions to plane systems by using the center manifold theorem, and then the limit cycles, centers, and isochronous centers are studied by means of the plane bifurcation theory.For system (2), there exists an approximation to the center manifold u = u(x, y), which can be represented as the polynomial series in x and y formally as follows: where h.o.t.denotes a higher-order term with an order greater than or equal to 3.
By transformation the system (2) can be transformed into the following complex system: in which z, w, T, a kjl , b kjl , dkjl ∈ C(k, j, l ∈ N).It is worth noting that the coefficients a kjl and b kjl of system (4) satisfy a conjugate relationship, namely, b kjl = a kjl , k ≥ 0, j ≥ 0, l ≥ 0, k + j + l ≥ 2. System (2) and system (4) are referred to as being concomitant.To make it easier to write dkjl and Ũ, we still use d kjl , U.
A singular point of a system becomes a center if and only if all its focal values are zero.So we have the following results.
Definition 1 (see [45]).For system (2), if v 1 = 1, then the origin is called the rough focus (strong focus); if v 1 = 1, and v then the origin is called the fine focus (weak focus) of order k, and the quantity of v 2k+1 is called the k-th focal value at the origin (k = 1, 2, . ..); if v 1 = 1, and for any positive integer k, v 2k+1 = 0, then the origin is called a center.Lemma 3 (see [44]).Let µ k be the k-th singular point value at the origin of system (4), and v 2k+1 be the k-th focal value at the origin of system (2 The maximum number of limit cycles bifurcated from the fine focus of the system is closely related to the highest order of the fine focus.So solving the highest-order fine focus question of the system is equivalent to analyzing the issue of the maximum number of limit cycle bifurcations as follows: Lemma 4 (see [46]).If the origin of system (2) is a fine focus of order k, then the origin of system (2) can bifurcate at most k small-amplitude limit cycles under a suitable perturbation.
Lemma 5 (see [45]).Suppose that the focal values v i of system (2) associated with the origin are determined by k independent parameters λ 1 , λ 2 , . . ., λ k , i.e., v i then system (2) has exactly k small-amplitude limit cycles bifurcating from the origin under a suitable perturbation on λ * .

Singular Point Values
For the Z q -equivariant symmetric system, it can be seen from Lemma 1 that its singular points appear as multiples of q.Moreover, these q singular points have the same topological structure, so their limit cycles are also generated in multiples of q, and the properties around the corresponding limit cycles are similar.Therefore, when we study the centerfocus problem and the limit cycle bifurcation of the Z q -equivariant symmetric system, we only need to select one of the q singular points with the same topological structure for research.We can easily determine that the points (1, 0, 0), (− 1 2 , 2 , 0) are three singular points of the system (1), which are symmetric about the Z-axis.So we only need to choose one of them, and at the same time for the convenience of calculation, we only calculate the singular point values of the singular point (1, 0, 0).
To calculate the values of a single non-origin singular point, the following steps are typically involved: First, the singular point is moved to the origin by a transformation.Second, the real system is changed into a complex system by complex transformation.Third, calculate the singular point value in the complex system.
Firstly, let and system (1) can be rewritten in the following form: Then, let and system ( 7) can be rewritten in the following form: ), Then, by using Lemma 2, we can calculate the singular point values of the singular point (0, 0, 0) of the complex system (8).With the aid of ycomputer algebra software Mathematica 12, the recurrence formula can be obtained as follows: Lemma 6.For system (8), the singular point values µ m (m = 1, 2, . ..) at the origin are determined by the following recursive formulas: where c k,k,0 = 0, k = 2, 3, . .., and if k < 0 or j < 0 or l < 0, we let c k,j,l = 0, else Next, using the computer algebra system Mathematica, the recursive formula is programmed, and the singular point values of system (8) corresponding to the Hopf critical point located at the origin are calculated, and the following theorem is obtained.Theorem 1.For the flow on the center manifold of system (8), the first five singular point values at the origin are given as follows: where Case 1 : if A 020 = 0, then B 030 = 0 and Case 2 : if A 020 = 0, then iA 020 H 3 , iA 020 H 4 , where H 4 and H 5 are given in Appendix A. For every µ k that was calculated, we already imposed that Because µ 3 , µ 4 and µ 5 are so complicated, we can simplify them appropriately.According to Theorem 1, the simplification of µ 3 , µ 4 , and µ 5 is equivalent to the simplification of H 3 , H 4 and H 5 .To simplify the calculation, we can first consider whether u 3 can be simplified.Obviously, with the program command Factor, PolynomialReduce, we can easily find the greatest common factor (G 3 ) of F 1 and H 3 .Then µ 3 = 1 124416 iA 020 H 3 can be reduced to µ 3 = − 1 331776 iA 020 F 2 by factorization.The same method is adopted for the simplification of µ 4 and µ 5 , so we can obtain the following theorem.Theorem 2. Through factorization and polynomial reduction, the first five order singular point values of system (8) can be reduced to the following form: Case 1 : if A 020 = 0, then B 030 = 0 and Case 2 : if A 020 = 0, then iA 020 F 3 , where F 3 and F 4 are given in Appendix A. For every µ k that was calculated, we already imposed that

Center Condition
In this section, we obtain the necessary and sufficient condition for the singular point (1, 0, 0) of system (1) to be a center by using the singular point values and the Darboux integral method.
From Theorem 2, we obtain the following result.
Proof.First, we prove its sufficiency, and substitute B 030 = − 3 2 (−3A 020 + 3A 010 A 020 + A 020 A 030 ), A 020 = 0 into the expressions of the first five singular point values in Theorem 2. It is clear that µ 1 = µ 2 = µ 3 = µ 4 = µ 5 = 0 is established.The sufficiency of Theorem 3 is proved.Next, we prove that the condition in Theorem 3 is also necessary.Taking the first singular point value µ 1 = 0, then we obtain B 030 = − 3 2 (−3A 020 + 3A 010 A 020 + A 020 A 030 ).Next, taking the second singular point value µ 2 = 0, then we obtain A 020 = 0 or A 020 = 0, F 1 = 0, so we need to have a case-by-case discussion.Case A : If A 020 = 0, that is B 030 = 0, then under this condition, Case 1 of theorem 2 holds, so this situation discussion is completed.Case B : If A 020 = 0, then Case 2 of Theorem 2 holds.Next, we need to determine whether µ 2 = µ 3 = µ 4 = µ 5 = 0 is true, that is, to discuss whether or not the polynomials F 1 , F 2 , F 3 and F 4 have common zeros.For this purpose, we apply the command "GroebnerBasis" in Mathematica to calculate the Gröbner basis of the < F 1 , F 2 , F 3 , F 4 > about the independent variables A 010 , A 020 and A 030 .We obtain This indicates that the polynomials F 1 , F 2 , F 3 and F 4 have no common zeros.Therefore, in Case 2, the first five singular values of system (8) (i.e., (7)) cannot be zero at the same time.In conclusion, the proof of Theorem 3 is complete.
It is necessary that the system's first finite focal values are equal to zero for a singular point to become a center.In order to obtain the necessary and sufficient center condition, its sufficiency needs to be proved.The common methods of proof include using the symmetry principle, finding the first integral means, determining the integral factor method, etc.It is not difficult to see that Theorem 3 gives the necessary condition for the singular point (0, 0, 0) to be a center of system (7).
Theorem 4. For system (7), if and only if B 030 = − 3 2 (−3A 020 + 3A 010 A 020 + A 020 A 030 ) and A 020 = 0, the surface Γ : u 1 = 0 is an invariant algebraic surface, which defines a global center manifold, and on this center manifold, the origin of system (7) is a center.
Proof.The necessity of Theorem 4 is obvious.Now we use the Darboux integrable theory [30,47] to demonstrate the sufficiency of this condition.
When condition(B 030 = − 3 2 (−3A 020 + 3A 010 A 020 + A 020 A 030 ) and A 020 = 0) in Theorem 4 holds, system (7) can be rewritten in the form It is known from [47] that for the following system where the degree of P(x, y, u),Q(x, y, u),R(x, y, u) does not exceed n.The polynomial equation F(x, y, u) = 0 is an invariant algebraic surface of system (10) if and only if there exists a polynomial K(x, y, u), which satisfies the following condition: where the polynomial F is called the Dardoux polynomial of system (10).F(x, y, u) is generally defined as where m is an integer not less than 1, F i (x, y, u) is a homogeneous polynomial of degree i, and F i (x, y, u) = 0.The polynomial K is the cofactor of F(x, y, u) = 0, and its degree does not exceed n − 1.
For system (9), in order to prove the existence of the Dardoux polynomial F(x 1 , y 1 , u 1 ) by searching from m = 1, for m = 1, let Since system (9) is of degree 3, the corresponding cofactor K(x 1 , y 1 , u 1 ) has a degree of at most two, and we assume that After we substitute ( 13) and ( 12) into (11) and obtain an algebraic equation, then comparing the coefficients of the same powers of the equation, we can obtain the following algebraic equations: Solving the algebraic equations ( 14), we can obtain c 1 = 0, c 2 = 0 and c 3 = 1, which results in F(x 1 , y 1 , u 1 ) = u 1 so that the surface Γ : F(x 1 , y 1 , u 1 ) = u 1 = 0 is truly an invariant algebraic surface of system (9).
Next, we will demonstrate that the surface Γ : F(x 1 , y 1 , u 1 ) = u 1 = 0 is actually a global center manifold of system (9).First, we calculate the normal vector of the surface F(x 1 , y 1 , u 1 ) = u 1 at the origin, and obtain ∇F(0, 0, 0) = (0, 0, 1).Since the eigenvalues of the linear matrix of system (9) at the origin are two purely imaginary eigenvalues ±i and a negative eigenvalue −d, it is obvious that the tangent space of the center manifold at the origin is spanned by these vectors: In addition, one can verify that ∇F(0, 0, 0) • e 1 = 0, ∇F(0, 0, 0) • e 2 = 0, which means that the surface F(x 1 , y 1 , u 1 ) = u 1 = 0 is truly a global center manifold of system (9).Finally, it is necessary to prove that the origin of system ( 9) is a center on the center manifold F(x 1 , y 1 , u 1 ) = u 1 = 0.When u 1 = 0, system (9) is simplified to the following planar system: It is easy to see that system (15) is a symmetric system about the X-axis, and according to the principle of symmetry, the origin of the system is a center (i.e., the singular point (0, 0, 0) of system ( 7) is a center).This completes the proof of Theorem 4.
Next, we use the system (9) to simulate the case of center at (0, 0, 0).All trajectories starting from the initial points on the invariant surface will remain on the surface.For example, when the initial points are chosen as (x 1 , y 1 , u 1 ) = (0.05, 0.02, 0), the periodic orbits are located on the invariant surface as shown in Figure 1.
However, if the initial points are chosen from outside the invariant surface, then all trajectories first converge to the invariant surface u 1 = 0, and once they reach the invariant surface, they become periodic orbits on the invariant surface, as shown in Figure 1, where one initial point is chosen below the invariant surface, and another initial point is chosen above the invariant surface given by (x 1 , y 1 , u 1 ) = (−0.1,−0.1, −0.01), (0.07, 0.03, 0.01).9), The blue line is the trajectory of the initial point (x 1 , y 1 , u 1 ) = (0.05, 0.02, 0), the yellow line is the trajectory of the initial point (x 1 , y 1 , u 1 ) = (0.07, 0.03, 0.01) and the red line is the trajectory of the initial point (x 1 , y 1 , u 1 )=(−0.1,−0.1,−0.01).

Bifurcation of Limit Cycles
In this section, we will discuss the maximum number of limit cycles for system (1) branching from the singular point (1, 0, 0).
From Lemma 4, we first need to consider the highest order of the fine focus of the origin of system (7), i.e., find the value of k such that v 3 = • • • = v 2k = 0, but v 2k+1 = 0 holds.From Theorem 1 and Lemma 3, we obtain the first nine focal values of system (7) at the origin: v 2m+1 = iπµ m (m = 1, 2, . . ., 5).
According to the discussion of Theorem 3 and Theorem 4, we can know that the value of k is 5, so the highest degree of the fine focus at the origin of system (7) is 5.Then, can the order of the fine focus at the origin of system (7) reach 5? That is, we discuss whether system (7) has such a set of parameter conditions that v Here is what we will study next.
From the proof of Theorem 3, as long as the condition (B 030 = − 3 2 (−3A 020 + 3A 010 A 020 + A 020 A 030 )) is met, we can simply find the values of the coefficients that make v 3 = 0 true.So, we now only need to look for the values of the coefficients such that v 5 = v 7 = v 9 = 0, v 11 = 0.According to the proof of Theorem 3, we know that the polynomials F 1 , F 2 , F 3 and F 4 have no common root, but F 1 , F 2 and F 3 should have.Then, finding the coefficient values such that v 5 = v 7 = v 9 = 0, v 11 = 0 holds is equivalent to solving the equations F 1 = 0, F 2 = 0, F 3 = 0 for finding the solutions of A 010 , A 020 , A 030 .With the help of Mathematica, we can obtain that this equation has four groups of real solutions, and take one of them as follows (only showing the first 50 digits): A 010 = 4.5629298768000231029699982427738544924493861105275 . . ., A 020 = 2.3050977788352747166914549291286457145886259511924 . . ., A 030 = −10.163059853074365253530894833667331792823007039006 . . ., (16) under which This result implies that there is a common solution such that This means that we found the values of the coefficients that make v 3 = v 5 = v 7 = v 9 = 0, v 11 = 0 to be true.Substitute the value of ( 16) in the condition B 030 = − 3 2 (−3A 020 + 3A 010 A 020 + A 020 A 030 ), obtaining Hence, there exists a set of parameter conditions of system (7) such that v 3 = v 5 = v 7 = v 9 = 0, but v 11 = 0, that is, the origin of system ( 7) is a fine focus of the fifth order, i.e., the singular point (1, 0, 0) of system ( 1) is a fine focus of fifth order.
To sum up, we have the following theorem.
Theorem 6.For the flow on the center manifold of system (7), the highest order of the fine focus at the origin is five, and the origin becomes a fine focus of the fifth order if and only if the following conditions are satisfied: where the polynomials F 1 , F 2 and F 3 have the same form as in Theorem 2. So for system (1), if and only if condition (18) is satisfied, the highest order of the fine focus is five, and the singular point (1, 0, 0) is a fine focus of fifth order.
By applying Theorem 6 and Lemma 4 to system (1), we can observe that, under suitable perturbations, the singular point (1, 0, 0) can generate at most five small amplitude limit cycles.Nonetheless, it remains to be determined whether the system indeed has five small amplitude limit cycles.The following theorem provides an affirmative answer to this question.Theorem 7.For system (1), when the coefficients satisfy the condition of Theorem 6, the singular point (1, 0, 0) is the fifth-order fine focus.Moreover, system (1) can bifurcate five limit cycles at the singular point (1, 0, 0) by suitable perturbations of the parameters.

Conclusions
In this paper, we utilize the formal series method to compute the values of singular points and investigate the bifurcation of limit cycles and centers in a specific class of three-dimensional cubic systems with Z 3 -equivariant symmetry.We not only establish a necessary and sufficient condition for the singular point of system (1) to transform into a center but also demonstrate that the system can give rise to 15 small amplitude limit cycles with a 5 − 5 − 5 distribution from the singular points following appropriate perturbation.It is worth noting that the obtained result of fifteen limit cycles represents a novel lower bound on the number of limit cycles in three-dimensional cubic systems with Z 3 -equivariant symmetry.
However, for more general three-dimensional symmetric systems, the problem remains open.Further investigations are required to explore additional dynamic properties, such as isochronous centers and periodic orbits.The challenges lie in calculating the quantities of singular points and determining the center conditions.Future research efforts may focus on enhancing computational tools to overcome these difficulties and achieve further improvements.