Study of chaos in rotating galaxies using extended force-gradient symplectic methods

We take into account the dynamics of three types of models of rotating galaxies in polar coordinates in a rotating frame. Due to non-axisymmetric potential perturbations, the angular momentum varies with time, and the kinetic energy depends on the momenta and spatial coordinate. The existing explicit force-gradient symplectic integrators are not applicable to such Hamiltonian problems, but the recently extended force-gradient symplectic methods proposed in a previous work are. Numerical comparisons show that the extended force-gradient fourth-order symplectic method with symmetry is superior to the standard fourth-order symplectic method but inferior to the optimized extended force-gradient fourth-order symplectic method in accuracy. The optimized extended algorithm with symmetry is used to explore the dynamical features of regular and chaotic orbits in these rotating galaxy models. The gravity effects and the degree of chaos increase with an increase of the number of the radial terms in the series expansions of the potential. There are similar dynamical structures of regular and chaotical orbits in the three types of models for the same number of the radial terms in the series expansions, energy and initial conditions.


I. INTRODUCTION
Chaos in a dynamical system means that the final state of the dynamical system displays exponentially sensitive dependence on the initial state. Based on the importance of chaos, many studies have focused on the subject of chaos in the solar system [1,2] and galactic dynamics [3][4][5][6][7].
Regular or chaotic motions of particles in galactic dynamics may affect the fraction of mass. The authors of Refs. [3,[5][6][7] investigated mass components of non-rotating N-body models of elliptical galaxies in ordered and in chaotic motion. Voglis et al. [3] found that the fraction of mass in chaotic motion is about 24% of the total mass in one non-rotating triaxial equilibrium model with smooth centers, and 32% in another non-rotating triaxial equilibrium model with smooth centers. This shows that the spatial distribution of the mass in chaotic motion is in disagreement with that in ordered motion. Muzzio et al. [7] pointed out that the fraction of mass in chaotic motion of about 53% in models of non-rotating galaxies with smooth centers. On the other hand, the fraction of mass in ordered or chaotic motion was also studied in spiral galaxies. Voglis et al. [8] further showed that rotation leads to increasing the fraction of mass in chaotic motion (up to the level of ≈ 65%) and shifting the Lyapunov numbers to larger values in N-body models of rotating galaxies. In other words, the extent of chaos is substantially enhanced by the rotation, and the fractions of mass in chaotic motion in the rotating models are larger than those in the non-rotating models. The spatial distributions of the dynamical structures along the spiral arms at the ends of the bar in the barred-spiral galaxy correspond to those of particles with masses in regular and in chaotic orbits. In fact, the mass in chaotic motion can almost completely form spiral arms emanating from the neighborhood of the Lagrangian points L 1 and L 2 at the ends of the bar in a barred-spiral galaxy. The distribution of mass must be associated with gravitational forces of particles. The gravity is the main driving mechanism for the formation and the stability of spiral arms in galaxies. The authors of Refs. [9,10] developed the Moser analytic series representing the invariant manifolds near the unstable Lagrangian equilibrium points L 1 , L 2 in a rotating barred galaxy. In this way, these series can represent the spiral arms, which are density waves and are composed of chaotic orbits. Besides the analytical theory, other methods such as the specific finite time Lyapunov characteristic number, the smaller alignment index, the surface of section and the frequency analysis were used to classify the orbits in regular and chaotic cases in barred galaxies [11][12][13][14][15].
The onset of chaos in two-dimensional Hamiltonian systems of rotating galaxies in the disc plane in polar coordinates [10] is due to the inclusion of non-axisymmetric potential perturbations. These perturbations cause the angular momentum with respect to the angle to vary with time and therefore the Hamiltonian systems are not integrable. This non-integrability leads to chaos under some circumstances. Reliable numerical results are always required to detect the chaotical behavior. In some cases, extremely long integration times are also required. The adopted computational schemes for the long-term integration of the Hamiltonian systems become crucial to reach better stability and higher precisions. The proper choice of the integrators should naturally be symplectic schemes, which preserve the symplectic nature of Hamiltonian dynamics [16,17]. Symplectic methods are a class of geometric integration algorithms [18] and make the local error in the total energy not grow with time. There are standard symplectic methods [19,20] that require some evaluations of the force, and force-gradient symplectic integrators [16,[21][22][23][24] that require some evaluations of force gradient in addition to several evaluations of the force. Because of good long-term behavior, the standard symplectic methods have been used in the solar system [1,25] and black hole spacetimes [26][27][28][29][30][31]. They are also suitable for the two-dimensional Hamiltonian systems of rotating galaxies in polar coordinates. However, the force-gradient symplectic integrators are not applicable to such Hamiltonian problems because the kinetic energies of the Hamiltonian systems depend on the momenta and spatial coordinates. In fact, they are only adapted to the integrations of Hamiltonian systems, where the kinetic energies are quadratic functions of momenta and the potential energies are functions of coordinates. Energy errors of the force-gradient symplectic integrators are several orders of magnitude smaller than those of the standard symplectic algorithms of the same order, as was confirmed in the literature [22][23][24]. Recently, our group extended these force-gradient integrators to the explicitly integrable kinetic energies, which are not only quadratic functions of momenta but also depend on coordinates [32]. When the original force-gradient operator is adjusted appropriately, the adjusted operator lacks the concept of force gradient and belongs to the momentum operator like the operator corresponding to the potential. As a result, the existing explicit force-gradient symplectic integrators [21][22][23][24] are still available in the extended Hamiltonian systems. The extended force-gradient symplectic integrators do not alter their symmetry and time reversibility compared with the existing force-gradient algorithms. The authors of Ref. [32] used a modified Hénon-Heiles system and a spring pendulum as two toy models to test the numerical performance of the extended algorithms. They showed that the fourth-order Forest-Ruth standard symplectic method does not give true dynamical properties of order and chaos to the modified Hénon-Heiles system under some circumstances, whereas the fourth-order extended force-gradient symplectic methods do. The obtained results are because the Forest-Ruth method performs poorer accuracy than the extended force-gradient algorithms. In fact, the optimized fourth-order extended force-gradient symplectic methods have energy errors that are three orders of magnitude smaller than those of the Forest-Ruth method.
Note that the extended force-gradient symplectic methods proposed in Ref. [32] were shown to have good numerical performance in the simulations of the two toy models. Now, we wonder whether the extended algorithms still exhibit excellent performance in the application of real astronomical and astrophysical problems. To evaluate the performance of the extended force-gradient symplectic methods applicable to the three types of models of rotating barred galaxies in Refs. [8,10] is one of the main aims in the present paper. On the other hand, we are interested in studying a distribution of dynamical structures regarding regular and chaotic orbits along the radial direction for a given number of the radial terms of the potential in one of the models using the techniques of Poincaré surface of section and fast Lyapunov indicators (FLIs) [33]. In this way, we desire to know how the number of the radial terms of the potential affects chaos in one of the models. We also desire to know what differences in the dynamics exist among the three types of models with the same number of the radial terms of the potential, energy and initial conditions.
To implement the three aims, we organize this paper as follows. In Section 2, we introduce three types of models of rotating barred galaxies. In Section 3, we apply the extended fourth-order force-gradient symplectic methods to Model A4 and evaluate the performance of these algorithms. In section 4, we explore the dynamical structures in these models using the techniques of Poincaré surface of section and FLIs. Finally, the main results are summarized in Section 5.

II. MODELS OF ROTATING BARRED GALAXIES
The motion of a test particle of mass m p = 1 in the plane of a galaxy with a rotating bar is described in polar coordinates (r, φ) in the rotating frame by a two-dimensional Hamiltonian [8,10] The above notations are given here. Ω is a pattern speed of the rotating frame. p r represents a canonical momentum vs the coordinate r, and p φ is a canonical momentum of the coordinate φ corresponding to the angular momentum in the rest frame. Φ(r, φ) denotes the gravitational potential of the galaxy in the rotating frame and satisfies the Poisson equation where G is the constant of gravity and ρ stands for the density of matter.
A complete expression of the solution of the Poisson equation is long and complicated. Harsoula et al. [10] took a simple solution of Eq. (2) as the potential for the m = 2 mode of the galactic bar Φ(r, φ) = Φ 0 (r) + Φ 1 (r) cos 2φ + Φ 2 (r) sin 2φ, where Φ 0 (r), Φ 1 (r) and Φ 2 (r) are functions of r. If Φ 1 (r) = Φ 2 (r) = 0, then the potential Φ 0 (r) = Φ(r, φ) is axisymmetric and the angular momentum p φ is conserved. When Φ 1 (r) = 0 or Φ 2 (r) = 0, the potential Φ(r, φ) is non-axisymmetric and the angular momentum p φ is not conserved. Namely, the second and third terms of Eq. (3) act as an m = 2 mode of the non-axisymmetric potential perturbation. The three potential functions are given in Ref. [10] by where R is a numerical constant for the description of the size of the system. The size is regarded as the N-body code boundary corresponding to the solutions of the Poisson equation matched with the solutions of the Laplace equation.
The other notations such as A 00 are expanded in terms of the series of spherical Bessel functions j 0 and j 2 [34]: In the above series expansions of the potential, n represents the number of the radial terms; that is, each of the functions like A 00 in the potential has (n + 1) terms. In addition, ξ il = a il r/R. For l = 0, a i0 = (i + 1/2)π is the (i + 1)th root of the equation j −1 (a i0 ) = 0. For l = 2, a i2 is the solution of the equation tan (a i2 ) = a i2 . The solution a i2 should be restricted in the range (i − 1/2)π < a i2 < (i + 1/2)π and is solved by the Newton iteration method. The coefficients B i00 , B i20 , C i21 and B i22 in Eqs. (7)-(10) are obtained through the N-body code on the positions of the N-body particles and Eq. (2). They are different for three types of galactic models, which are called as models A, B and C in Ref. [10]. The properties of the three models were described in the papers [8][9][10]. Here are some of them. Firstly, the total angular momentum is larger in model B than in model A, but smaller than in model C. This leads to decreasing the size of a rotating central bar formed by density waves and increasing the density in the region of the bar along the sequence of the models A, B and C. Secondly, the pattern speed of the bar in each of the models becomes smaller, whereas the corotation radius gets larger at the end of a Hubble time. Thirdly, the degree of chaos is enhanced by rotation increasing the fraction of mass in chaotic motion. Fourthly, mass in chaotic motion almost completely dominates the formation of spiral arms. Fifthly, invariant manifolds of all unstable periodic orbits near and beyond corotation support both the outer edge of the bar and the spiral arms. Harsoula provided the coefficients in the three models to the first author of the present paper via private communication. The coefficients with n = 19 are listed in Tables 1-3. The 20 coefficients B i00 , the 20 coefficients B i20 and the 40 coefficients C i21 , B i22 (i = 0, · · · , 19) in each model are associated to monopole terms, quadrupole terms and triaxial terms, respectively. In practice, each group of the coefficients determine the potential Φ(r, φ). Different potentials Φ(r, φ) correspond to different galactic models. Hereafter, the three models for n = 19 respectively correspond to models A19, B19 and C19. Similarly, models A9, B9 and C9 are called when n = 9; models A4, B4 and C4 are also called when n = 4.
The unit systems are those of Refs. [8,10]. The half mass radius R h is used as a scaling unit of length, i.e. r scal = r/R h , where r stands for the real distance and r scal is the scaling radial distance. For convenience, the scaled radial distance is still written as r in the later discussions. The half mass radius is R h = 0.1006 for Model A, R h = 0.0926 for Model B and R h = 0.1167 for Model C. The size of the system is R = 0.85. Here, one unit of length is 8kpc. The unit of time is the half mass crossing time T hmct = [2R 3 h /(GM )] 1/2 = t Hub /300, where t Hub represents a Hubble time. The pattern speeds of the rotating frame in the three models A, B and C are taken as Ω A = 5886.65, Ω B = 6010.36 and Ω C = 6137.14, which correspond to 20 ∼ 25 km sec −1 kpc −1 in real units.
As is mentioned above, the second and third terms of Eq. (3) destroy the axial-symmetry of the system (1) such that the angular momentum p φ varies with time. Thus, no additional constants of motion but the Hamiltonian H of Eq. (1) exists. This implies the nonintegrability of the system (1) with two degrees of freedom in a four-dimensional phase space. A numerical method is a good tool to solve this nonintegrable problem.

III. A CHOICE OF NUMERICAL INTEGRATOR
A prior choice to an integrator for a long-term integration of the Hamiltonian (1) is naturally a symplectic method with the conservation of the Hamiltonian flow. Following this idea, we consider the application of symplectic integration to this Hamiltonian problem.

A. Generalized force-gradient symplectic integrators
For convenience, we take p = (p 1 , p 2 ) = (p r , p φ ) and q = (q 1 , q 2 ) = (r, φ). We rewrite the Hamiltonian (1) as where the kinetic energy T is a function of the momenta p and coordinate r The lie derivative operators of T and Φ are defined as where the symbols {, } represent the Poisson brackets. Obviously, T and Φ can be exactly, analytically solvable, and A, B correspond to their solvers. The solvers A and B can symmetrically compose a second-order symplectic leapfrog integrator where τ is a time step, and two commutators are In the second line of Eq. (15), the first term τ (A + B) corresponds to the numerical solution of the Hamiltonian (1), and the second and third terms correspond to local truncation errors of the numerical solution remaining at an order O(τ 3 ). Because of such local truncation errors, the algorithm M2 can give a second-order accuracy to the explicit numerical solution. In terms of the solvers A and B, a fourth-order symplectic method of Forest and Ruth [19] is established by where β = 1/(2 − 3 √ 2) and α = β/2. If r 2 in the second term of Eq. (12) is absent, A of Eq. (13) is a position operator which acts on only the position coordinates. C = 2BAB is similar to the momentum operator B, which acts on only the momenta. In fact, C is an exactly, analytically solvable operator corresponding to the force gradient of the gravitational potential. However, D is a momentum and position mixed operator and is not an exactly analytical solver. In view of these facts, a composition of B and C appears in a class of explicit symplectic integrators, called the force-gradient symplectic integration algorithms [16,[21][22][23][24]. Now, r 2 in the second term of Eq. (12) remains. Although A of Eq. (13) is a momentum and position mixed operator, it is still an exactly analytical solver in this case. Is C is an exactly analytical solver? We said yes as an answer to this question in our previous work [32]. In fact, the answer can be given through the operators acting on the momenta and positions in the second line of Eq. (16). It is easy to check that ∂qi∂pj ∂p k . Clearly, C is also a momentum operator: The momentum operator C of Eq. (21) is no longer the force gradient of the gravity potential, but can be regarded as an extension of the force-gradient factor in Refs. [16,[21][22][23][24]. Without doubt, it is an exactly analytical solver similar to the momentum operator B. Symmetric composition methods of the three operators A, B and C can yield various explicit symplectic algorithms. For instance, a five-stage fourth-order symplectic method including the operator C was proposed by Chin [22] in the form Omelyan et al. [24] also constructed a seven-stage fourth-order optimized symplectic algorithms with the inclusion of the operator C: where the time coefficients are The concept of optimization is the time coefficients minimizing the norm of the leading term of fifth-order truncation errors.
The construction mechanisms of the algorithms with the extended operator C was discussed in the previous work [32]. In fact, the extended force-gradient algorithms acting on the original system (1) are equivalent to the standard symplectic methods (like the method M4 without the extended operator C) acting on the modified Hamiltonian systems.

B. Numerical tests
Model A4 is used as a test model to evaluate the performance of the algorithms M4, N4 and N4P. For comparison, an eighth-and ninth-order Runge-Kutta-Fehlberg integrator [RKF8 (9)] with adaptive step sizes is also employed. The time step is τ = 0.0005 × T hmct = 0.0005 × t hub /300 = 0.0005/82.4/300 ≈ 2.0 × 10 −8 (sec·Mpc/km)=2.0 × 10 −5 (sec·kpc/km) ≈ 2.0 × 10 −5 (9.5 × 10 8 y) = 0.19 × 10 5 y. In short, τ ≈ 2.0 × 10 −5 (sec·kpc/km) is used in our codes. The authors of Ref. [8] pointed out that the time step is a good choice because the hysteresis between the bar of the potential and the bar of the real density of particles during the time is minimal and the cumulative effect of a numerical retarding torque in a Hubble time is small.
The energy of the system (1) is E = H = −7 × 10 6 . The initial conditions are r = 0.164/R h = 0.164/0.1006, φ = 1.89π, p r = 0. The initial momentum p φ > 0 is determined by E = H. Fig. 1(a) plots the relative energy errors |∆H/H 0 | = |(E t −E)/E| for these several algorithms, where E t is the numerical energy at time t. Here, each algorithm has 2.5 × 10 8 steps, which correspond to time t ≈ 5056.63 (sec·kpc/km). The three symplectic methods M4, N4 and N4P show no secular growth in the energy errors. This is an inherent advantage of these symplectic methods. Among the three symplectic integrators, the standard fourth-order symplectic method M4 without the extended operator C performs the poorest accuracy, whereas the fourth-order optimized symplectic algorithm N4P with the extended force-gradient operator C exhibits the best accuracy. These numerical results show that the inclusion of the extended force-gradient operator has an advantage over the exclusion of the extended force-gradient operator in the accuracy. The optimized method is also better than the corresponding non-optimized method. Unlike these symplectic methods, the non-symplectic integrator RKF8(9) makes the energy error grow with time. RKF8(9) is superior to N4P in the accuracy, but inferior to N4P in the computational efficiency, as is shown in Table 4. Taking the solutions of RKF8(9) as reference solutions, we obtain the relative position errors of the three symplectic algorithms M4, N4 and N4P in Fig. 1(b). As a result, the method M4 still has the largest position error, while the method N4P yields the smallest one. When many other values of the energy and initial conditions are also used, they do not affect the numerical performances of these integrators. In fact, the numerical performance of an integrator is independent of a choice of the parameters and initial conditions [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. In spite of this, it is necessary to choose bounded orbits as testing the integrator's performance. Thus, an appropriate choice of the parameters and initial conditions is still necessary.
Because of the fourth-order optimized symplectic algorithm N4P with the extended force-gradient operator C having the best performance in the stabilization of energy errors, it is used to study the dynamics of orbits.

IV. REGULAR AND CHAOTIC DYNAMICS
The authors of Ref. [10] applied the Moser theory of invariant manifolds to the chaotic spiral arms. The invariant manifolds starting at the Lagrangian points L 1 , L 2 , or unstable periodic orbits around L 1 and L 2 are described in terms of series. The convergence of the series is an analytical means for studying chaotic orbits with initial conditions in the neighborhood of the invariant manifolds of the unstable point L 1 or L 2 . A domain of the convergence of the series around every Lagrangian point is a Moser domain. The intersection of the orbits inside the Moser domain with an apocentric section produces the spiral structure of this domain. The chaotic orbits with initial conditions near, but in the exterior of the boundary of the Moser domain of convergence become chaotic attractors. Such chaotic orbits are also identified by the techniques of Poincaré sections and fast Lyapunov indicators (FLIs) [33]. The method of Poincaré sections, which shows intersections of the particles' trajectories with the surface of section in phase space, can clearly describe the phase space structure of a conservative 4-dimensional system. One or several isolated points on the Poincaré sections correspond periodic orbits. Kolmogorov-Arnold-Moser (KAM) tori on the Poincaré sections correspond regular quasi-periodic orbits. If there are many plotted points that are distributed randomly in an area, the motion is chaotic. In a word, the distribution of the points in the Poincaré map can show whether or not the motion is chaotic. For ordered and chaotic orbits, the length of a deviation vector increases in completely different time rates. The FLI uses the completely different time rates to distinguish between the ordered and chaotic cases. The Moser theory of invariant manifolds for allowing to study chaotic orbits requires that the series describing the Hamiltonian dynamics near an unstable equilibrium point or an unstable periodic orbit be convergent. It may not work well for the non-convergent Birkhoff normal form series around stable invariant points, or stable periodic orbits. The techniques of Poincaré sections and FLIs for finding chaos do not have this restriction. Clearly, they can detect chaotic orbits in larger regions compared with the Moser theory of invariant manifolds. Here, we employ the techniques of Poincaré sections and FLIs to investigate the dynamical behavior of the three types of bar spiral galaxy models A, B and C. The considered orbits are not restricted to those around the unstable Lagrangian points L 1 and L 2 . We are mainly interested in comparing the differences in the dynamical behavior among these types of models. The effect of the number (n + 1) of the radial terms in the series expansions of the potential on the dynamical behavior in each set of models is also considered.

A. Models A4, A9 and A19
The parameters and initial conditions are those of Fig. 1, but the initial values of r and p φ are different. Model A4 exhibits different phase space structures in two ranges of the initial separations r, as is shown through Poincaré sections at the plane φ = π/2 with p φ > 0 in Figs. 2 (a) and (b). The three orbits in Fig. 2(a) are KAM tori, which correspond to the regularity of the orbits. One orbit with the initial separation r = 4.430/0.1006 and another orbit with the initial separation r = 4.424/0.1006 in Fig. 2(b) seem to have no explicit difference from the phase space structures on the section, but they exhibit different regular and chaotic dynamical features, which can be identified clearly by means of the technique of FLI in Fig. 2(c). The FLI is that with two nearby orbits defined in Ref. [35]: where d(0) and d(t) are the distances between two nearby orbits at time 0 and t, respectively. The orbit with the initial separation r = 4.430/0.1006 corresponds to the FLI increasing in a power law with time log 10 t and is regular. However, the orbit with the initial separation r = 4.424/0.1006 has the FLI increasing in an exponential law with time and should be chaotic. After 5 × 10 4 integration steps, the FLIs with completely different increasing laws with time are sensitive to distinguish the two cases of order and chaos. The phase space structures for Model A9 in Fig. 3(a) are somewhat unlike those for Model A4 in Fig. 2. When the initial separation r is given in a small range of 1.45 < r < 1.7 for Model A4, the orbits are regular KAM tori (not plotted). However, the orbit with the initial separation r = 0.153/0.1006 in Model A9 has several hyperbolic points. Its chaoticity is confirmed by the FLI of Fig. 3(c). There are no chaotic orbits in a range of 10 < r < 12.3 for Model A4 in Fig. 2(a), but there are two chaotic orbits with the initial separations r = 1.210/0.1006 and r = 1.235/0.1006 in the range of 10 < r < 12.3 for Model A9 in Fig. 2 (b) and (c).
When Model A19 is considered in Fig. 4(a), an orbit with the initial separation r = 0.163/0.1006 is a regular KAM torus. Another orbit with the initial separation r = 0.154/0.1006 consists of many islands and therefore is still regular. However, the third orbit with the initial separation r = 0.160/0.1006 is chaotic because it has many discrete points which are randomly filled with an area. The degree of chaos seems to be strengthened in the small range of 1.45 < r < 1.7 from A4 to A9 and to A19. Chaos in the range of 10 < r < 12.5 seems to be stronger for A19 in Fig. 4(b) than for A9 in Fig. 3(b). When the initial separation r spans from 20 to 23 in Fig. 2(c), a lot of regular KAM torus orbits exist in the interior of the ordered orbit with the initial separation r = 2.230/0.1006. There is a small chaotic region between the two orbits of the initial separations r = 2.230/0.1006 and r = 2.290/0.1006. The regularity of the orbit with the initial separation r = 2.230/0.1006 and the chaoticity of the orbit with the initial separation r = 2.290/0.1006 can be identified clearly by means of the technique of FLIs in Fig. 4(d).
Seen from Figs. 2-4, the number and degree of chaotic orbits seem to increase as the number of the radial terms n gets large. This result is clearly shown in Fig. 5 that lists the dependence of FLI on the initial separation r in each of Models A4, A9 and A19. In other words, Fig. 5 shows that a distribution of the phase space structures of regular and chaotic orbits for each model depends on a range of the initial separation r. Each of the FLIs is obtained after 5 × 10 4 integration steps. 7 is a critical value of the FLI between the ordered and chaotic two cases. The FLIs not more than 7 indicate the regularity of bounded orbits, whereas those larger than 7 correspond to the chaoticity of bounded orbits. The result on an increase of n increasing the number and degree of chaotic orbits is based on n corresponding to the number of the radial terms in the series expansions of the potential. An increase of n means that of the number of the radial terms. Equivalently, the nonlinear gravitational interaction effects in the potential are enhanced.
B. Models B4, B9 and B19 Fig. 6 describes the phase space structures Model B4 in two ranges of the initial separations r. There is a weak chaotic orbit with the initial separation r = 0.750/0.0926 in the range of 7.2 < r < 8.2 in Fig 6(a). The two orbits in the range of 10.0 < r < 12.2 are regular KAM tori in Fig 6(b). However, chaos occurs for the initial separation r = 3.746/0.0926 in the range of 40.25 < r < 43 in Fig. 4(c). The dynamical features of the chaotic orbit in panel (a) and the two orbits in panel (c) are also shown in Fig. 6(d). The phase space structures of Model B4 in Fig. 6 (b) and (c) are similar to those of Model A4 in Fig. 2 (a) and (b).
The phase space structures of Model B9 in Fig. 7 (a) and (b) also look like those of Model A9 in Fig. 3 (a) and (b). The initial separations r = 0.157/0.0926 and r = 0.139/0.0926 correspond to two regular orbits consisting of three loops in Fig. 7(a). Chaos is present for r = 0.141/0.0926 in Fig. 7(a) and r = 1.128/0.0926 in Fig. 7(b), as the FLIs of Fig. 7(c) show.
Model B19 seem to have stronger chaos in the ranges of 1.425 < r < 1.7 and 10 < r < 12.25 in Fig. 8 (a) and (b) than Model B9 in Fig. 7 (a) and (b). For the initial separation r belonging to the range of 1.425 < r < 1.7, the orbital dynamical behavior for Model B19 in Fig. 8(a) is greatly different from that for Model A19 in Fig. 4(a). However, the phase space structures of Model B19 in Fig. 8 (b) and (c) also resemble those of Model A19 in Fig. 4 (b) and (c). The FLIs of Fig. 8(d) support the chaoticity of the two orbits for r = 2.368/0.0926 and r = 2.420/0.0926 in Fig.  8(c).
The relation between the FLI and the initial separation r in Fig. 9 clearly shows that the number and degree of chaotic orbits in Models B increase with n increasing.

C. Models C4, C9 and C19
As far as Models C are concerned, chaos seems to get stronger in the range of 3.5 < r < 5.1 from n = 4 in Fig.  10(a) to n = 9 in Fig. 10(c) and n = 19 in Fig. 10(e). The result is also suitable for the range of 10 < r < 12.2 in Fig. 10 (b), (d) and (f). The relation between the FLI and the initial separation r in Fig. (11) supports that chaos becomes stronger with n increasing.
Several points can be concluded from the above demonstrations. The three types of models A, B and C have the same expressions, but have differences in the pattern speeds Ω and the coefficients B i00 , B i20 , C i21 and B i22 . The three models have similar dynamical behaviors when they have same radial term number n, energy E and initial conditions. The dynamical structures of ordered and chaotic orbits in each model have different distributions along the radial direction. When the radial term number n increases, the gravity effects are enhanced. This leads to strengthening the degree of chaos. It is worth noting that the main structures of the spiral arms in the models A, B and C are formed at radial distances between 3 and 5 scaled radii, as is shown in Refs. [8] and [10]. Some ranges of the scaled radial distances, such as the range of the scaled radial distances from 40 to 43, seem to be very far from the main structures of the spiral arms. These results show that chaos can occur at the radial distances close to the main structures of the spiral arms and those far from the main structures of the spiral arms. In addition, the chaotic orbits we find are not restricted in the regions near the unstable equilibrium points L 1 and L 2 . Thus, the energies selected in this paper are unlike those of [8][9][10].

V. CONCLUSIONS
This paper mainly focuses on three classes of models of rotating galaxies in the polar coordinates in the rotating frame. They are models A, B and C, which have the same expressions but have different pattern speeds Ω and different coefficients B i00 , B i20 , C i21 and B i22 . When the potentials with large deviations from axial symmetry are included in these models, the angular momentum p φ is not a constant of motion. In this case, the kinetic energy is a function of the momenta and spatial coordinate.
The existing explicit force-gradient symplectic integrators [21][22][23][24] do not work for the present Hamiltonian problems. However, our extended force-gradient symplectic methods in the previous work [32] are still available. Numerical tests show that the fourth-order symmetric symplectic method without the extended force-gradient operator performs poorer accuracy than that with the extended force-gradient operator. In particular, the optimized extended forcegradient symplectic method is superior to the corresponding non-optimized one in accuracy.
The fourth-order optimized symplectic algorithm with the extended force-gradient operator is applied to survey the dynamical features of regular and chaotic orbits in these rotating galaxy models. It is shown through the techniques of Poincaré sections and fast Lyapunov indicators that an increase of the radial term number of the potential strengthens the gravity effects and the degree of chaos. The three types of models have similar dynamical structures for the same radial term number, energy and initial conditions. Data Availability Statement: Our paper is a theoretical work. All of the data are calculated and given in the paper.
Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. In accuracy, M4 is the poorest one among these symplectic integrators, while N4P is the best one. (b) The relative position error ∆r/r0 = (r − r0)/r0 between the radial separation r0 given by RKF8(9) and the radial separation r given by one of the symplectic methods. The position error for N4 is smaller than for M4 but larger than for N4P.   The structures are considered in the ranges of 1.5 < r < 1.7 (a) and 10 < r < 12.4 (b). The initial conditions in Fig. 3(b) are the same as those in Fig. 2(a). However, the orbits for r = 1.      Fig. 3 but Model B9 is considered. The phase space structures in panels (a) and (b) respectively look like those in Fig. 3 (a) and (b). The initial conditions in Fig. 7(b) are the same as those in Fig. 6(b). The orbit for r = 1.128/0.0926 is chaotic, but the orbit for r = 1.120/0.0926 is regular.   Fig. 4(a), but these structures in panels (b) and (c) respectively resemble those in Fig. 2 (b) and (c). The initial conditions in Fig. 8(b) are the same as those in Fig. 6