Chaos in $SU(2)$ Yang-Mills Chern-Simons Matrix Model

We study the effects of addition of Chern-Simons (CS) term in the minimal Yang Mills (YM) matrix model composed of two $2 \times 2$ matrices with $SU(2)$ gauge and $SO(2)$ global symmetry. We obtain the Hamiltonian of this system in appropriate coordinates and demonstrate that its dynamics is sensitive to the values of both the CS coupling, $\kappa$, and the conserved conjugate momentum, $p_\phi$, associated to the $SO(2)$ symmetry. We examine the behavior of the emerging chaotic dynamics by computing the Lyapunov exponents and plotting the Poincar\'{e} sections as these two parameters are varied and, in particular, find that the largest Lyapunov exponents evaluated within a range of values of $\kappa$ are above that is computed at $\kappa=0$, for $\kappa p_\phi<0$. We also give estimates of the critical exponents for the Lyapunov exponent as the system transits from the chatoic to non-chaotic phase with $p_\phi$ approaching to a critical value.


Introduction
Recently, there has been growing interest in exploring the structure of chaotic dynamics emerging from the matrix quantum mechanics [1][2][3][4][5][6][7][8][9][10][11][12], such as the BFSS and the BMN models [13][14][15][16][17][18][19] which appear in the DLCQ quantization of M-theory in the flat and the pp-wave background, respectively. These models are SU (N ) gauge theories, describing the dynamics of the N -coincident D0-branes, in the flat and spherical backgrounds. It is well-known that the gravity dual is obtained in the 't Hooft limit, i.e. at large N and strong Yang-Mills (YM) coupling and describes a phase in which D0-branes form a so called black-brane, i.e. a string theoretical black hole [18][19][20]. While the earlier investigations (and some recent as well [21][22][23][24][25][26] ) on the quantum mechanical behaviour of these models were performed in the Euclidean time formulation using both analytical perturbative and Monte-Carlo methods, in the past few years, there has been increasing interest on accessing the quantum dynamics using real-time formulations [10,11]. These studies are propelled by a result due Maldacena-Shenker-Stanford (MSS) [6], which states that under general circumstances, the Lyapunov exponent (which is a measure of chaos in both classical and quantum mechanical systems) for quantum chaos is bounded, that this bound is controlled by the temperature of the system, and given by λ L ≤ 2πT . It is conjectured that systems which are holographically dual to the black holes, are expected to be maximally chaotic. This is already demonstrated for the Sachdev-Ye-Kitaev (SYK) [27] model, and expected to be so for the BFSS model too. Numerical studies reported in [4] found that, for the BFSS model treated at the classical level, the largest Lyapunov exponent is given as λ L = 0.2924(3)(λ t Hoof t ) 1/4 . This is parametrically smaller than the MSS bound 2πT and violates it only temperatures below ≈ 0.015, while the quantum correction recently evaluated using Gaussian state approximation [11], indicate that the largest Lyapunov exponent vanishes below a non-zero temperature, and hence ensuring that the MSS bound is not violated.
It is important to note that not only the BFSS, BMN matrix models, but even their subsectors at small values of N appear as non-trivial many-body systems, and we lack a complete solution to these or even for the smallest Yang-Mills (YM) matrix model to date. The latter may be described as being composed of two 2 × 2 Hermitian matrices with SU (2) gauge and SO(2) global symmetries. It can be obtained by dimensionally reducing the YM theory from 2 + 1-to 0 + 1-dimensions. Classical dynamics of this system was recently investigated in [5] (see also the references [28,29] in this context) and it was shown that, using the SU (2) gauge and SO(2) rotations of the two matrices among themselves and a judicious choice of coordinates to fully implement the Gauss law constraint leads to a Hamiltonian with two degrees of freedom and their conjugate momenta. In addition, the angular momentum, p φ , associated to the rigid SO(2) symmetry appears as a conserved quantity via a term proportional to the square of p φ and strongly controls the structure of the effective potential and the ensuing dynamics. At p φ = 0, the model collapses to the usual x 2 y 2 potential, which is already known to lead to almost completely chaotic dynamics [30][31][32][33]. In [5], the response of the system to a range of different values of p φ is investigated and it is found that, at fixed energy, there is a value of p φ above which the chaos ceases to exist and the dynamics is essentially described by quasi-periodic motion. Therefore, the model is conjectured to have two phases, namely a chaotic phase corresponding to a toy model for a black hole, and a phase consisting of two D0-branes tied with fixed number of open strings stretching between them, with a force that depends on the number of excited strings. The latter can be roughly thought of as the "adiabatic invariant" for the quasi-periodic orbits, which appear as the Kolmogorov-Arnold-Moser (KAM) tori (see, for instance, [36]) in the Poincaré sections. For a given value of energy these two phases can coexist within a range of values of p φ , while the end of chaotic dynamics is argued to correspond to the end of the black hole phase. Quantum aspects of the 2 × 2 matrix model is addressed in [34], where the ground state energy is also estimated.
In order to gain further insight into the matrix model composed of 2 × 2 matrices with SU (2) gauge symmetry, in this paper, we set out to investigate the dynamics in the presence of the Chern-Simons (CS) term. It is possible to obtain the corresponding action starting from the SU (2) Yang-Mills Chern-Simons (YMCS) model in 2 + 1 dimensions and reducing it to 0 + 1 1 . In a manner similar to the one followed in [5], while paying attention to the differences in the procedure due to the CS term, which is first order in time derivative, we obtain the Hamiltonian of the system. The latter has the same degrees of freedom as the pure YM model, while the effective potential is governed not only by p φ , but also the CS coupling κ, which enters into the effective potential via κp φ and another term ∝ κ 2 r 2 . Varying κ at different values of p φ , we probe the impact on the chaotic dynamics. Our new findings are as follows. Firstly, we find that at p φ = 0, values of the largest (and only) Lyapunov exponent are above that evaluated at κ = 0, approximately within the range of values of 4π|κ| 4. This can be attributed to shrinking of the sharp edges of the effective potential contours (see figure 1a), but not sustained further for 4π|κ| 4 as the harmonic term ∝ κ 2 r 2 starts to dominate and chaotic dynamics gradually declines. The second and more interesting effect is due to the κp φ term, which alters the Lyapunov spectrum depending on its sign, in other words, the orientation of p φ matters. For instance, we find the values of the largest Lyapunov exponent for κp φ < 0 for a range of values of κ at fixed p φ are above that is evaluated at κ = 0. These results are presented and discussed in detail in section 3, where our findings obtained from the Lyapunov data are further corroborated via the use of Poincaré sections. As another important finding, we give estimates for the critical exponents for λ L and the value of the order parameter, p c φ , as the system transits from chaotic to non-chaotic phase. Rest of the paper is organized as follows. Section 2 gives the developments leading to the Hamiltonian of the model. Most of the details of the calculations in this section are relegated to the appendices for completeness. We summarize our results and briefly state our conclusions in section 4.

SU (2) Matrix Model with the Chern-Simons Term
The action of the model may be given as

2)
and In these expressions X 1 , X 2 are 2×2 traceless Hermitian matrices whose entries are functions of time only. They transform under the adjoint representation of SU (2): are the covariant derivatives, and A 0 is a gauge field which transforms accordingly under the local SU (2) gauge group. S is invariant under the local SU (2) gauge symmetry as well as under a global SO(2); i.e. the "rigid" rotations of the X i 's among themselves. In (2.3), κ is the CS coupling constant. Note that, due the gauge invariance of S CS term in 0 + 1-dimensions, κ is not level quantized 2 . Let us also note that we have implicitly set the Yang-Mills coupling, an overall factor 1 g 2 in S Y M , to unity. YM coupling can easily be restored back in the action by performing the scalings t → g 3 κ. Pure CS model limit is obtained by letting g → ∞ and it is discussed in detail in the appendices C and D.
It is convenient to work in the A 0 = 0 gauge. In the presence of the CS term, Gauss law constraint takes the form We may express the matrices X i as where 1 √ 2 is a normalization factor and σ α are the usual Pauli matrices. For future notational convenience, it is also useful to arrange components of X i into column vectors (2.6) Substituting (2.5) into the action (2.1) yields the Lagrangian

7)
while the constraint (2.4) takes the form The canonical conjugate momenta are easily obtained from the Lagrangian (2.7) as which clearly show that the kinematical and conjugate momenta are no longer the same in the presence of the CS term; a fact which is widely known in the literature (see for instance, [39]). Defining

10)
Gauss law constraint in (2.8) can be expressed as the condition of the vanishing of the SU (2) angular momentum: In order to obtain the corresponding Hamiltonian, we need to observe that the Lagrangian involves a term which is first order in time derivatives. Let us note that the generic form of such a Lagrangian can be given as where g ab is the metric associated to the generalized coordinates q a , f a are functions of the generalized coordinates, i.e. f a ≡ f a (q b ) and V ≡ V (q a ) is the potential. Corresponding Hamiltonian can be shown to take the form (see the appendix C for details) Adapting (2.13) to (2.7), in the Cartesian coordinates we obviously have g ab as the Euclidean flat metric δ ij , we may write f i = −κε ij x j (i, j : 1, 2) and observe that V = ( Putting all these together, we find that the Hamiltonian corresponding to (2.7) takes the form 3 with the equations of motion easily evaluated to bė (2.15) 3 In the pure CS limit the Hamiltonian becomes zero as explained in the appendices C and D.
Using (2.15), the time derivative of L 1 may be expressed aṡ

16)
A similar result holds for˙ L 2 . Although the second term in (2.16) remains aligned with L 1 as it does in the pure YM matrix model, this is not manifest for the first term. Nevertheless, the subsequent analysis will show, upon implementing the Gauss law in appropriate coordinates, that the dynamics remain planar. Taking advantage of the local SU (2) ≈ SO(3) and the global SO(2) rotations, we may introduce the coordinates (α, β, γ, r, θ, φ). Following [5], we may consider the 3 × 2 matrix M whose columns are the vectors x 1 and x 2 , i.e. M = ( x 1 , x 2 ) and express M as where R(α, β, γ) is a SO(3) Euler matrix using z − x − z active rotation with the angles (α, β, γ), respectively. Its explicit form is given in the appendix for quick reference.
:= (r cos θ, r sin θ, 0) may be thought as a configuration of the two D0-branes oriented coplanarly with a relative angle θ obtained via a SU (2) gauge choice. The latter is not preserved in general by the global SO(2) rotations on x 1 and x 2 , which can be taken to act on the right of M 0 , nor it is preserved by the SU (2) ≈ SO(3) gauge rotations, which acts from the left on M 0 . Thus, taking these facts together, (2.17) is a convenient way to introduce new coordinates for the present dynamical system. The advantage of this choice of the coordinates is that, the Gauss law constraint in (2.11) can be fully solved and manifestly imposed on the Hamiltonian expressed in terms of the new variables, as we will demonstrate in what follows. Let us also remark that, this is essentially the same approach followed in [5] except that, we no longer restrict the gauge SU (2) ≈ SO(3) rotations to an SO(2) subgroup in advance, since, it is not readily seen thaṫ L i (i = 1, 2) remain aligned with L i .
Let us note in advance that the components of angular momentum L can be expressed in terms of the conjugate momenta (p α , p β , p γ ) corresponding to the Euler angles (α, β, γ) 4 as [42] which immediately implies that the Gauss law constraint L = 0 is equivalent to We will make use of (2.19) to fully impose the Gauss law in what follows. 4 As this is not frequently encountered in the literature, we provide a quick derivation in the appendix E.
The metric in the new coordinates (r, θ, φ, α, β, γ) is straightforwardly obtained from the expression We give the components of g ij and its inverse g ij in appendix D and also provide there the details of the evaluation of the Hamiltonian in the new coordinates using the generic form in (2.13) together with the inverse metric g ij . Employing these facts and imposing the Gauss law constraint (2.19), we find Since this Hamiltonian is cyclic in φ, as in the pure YM case [5], p φ is a constant of motion and taking advantage of this fact, we have defined the effective potential, V ef f , in the second line of (2.21). A number of remarks regarding this Hamiltonian are now in order 5 . Firstly, we observe that the terms involving the CS coupling κ are new and therefore we are now in a position to examine the chaotic dynamics emerging from (2.21) as κ and the angular momentum p φ assume a range of different values. Also note the presence of the r term in V ef f . In [5] this term is motivated by the fact that for sin θ ≈ θ, the motion can be considered to be adiabatic in θ with an effective frequency ω θ,ef f ≈ r. With taken as a small parameter, the term r can then be considered as the quantum mechanical correction to the energy, which lifts the flat direction of the pure YM model, that is, the case corresponding to the commuting matrices. In the present case, dependence of V ef f on θ is the same as the pure YM model, leading to the same interpretation for this term. Though, the interesting new fact is that, for κ = 0, V ef f already develops a minimum even at = 0. This minimum is at θ = 0, and the real positive root of the quartic equation κ 2 r 4 + r 3 − p 2 φ = 0. For = 0, we obtain r 2 = p φ κ , which yields E > 2κp φ for κp φ > 0 and simply E > 0 for κp φ < 0. Let us also note that for κ = 0, r ∝ p 2/3 φ −1/3 , and for a typical value of = 0.1, V ef f ≈ 0.32 at p φ = 1 [5], while for κ = 0, this minimum shifts upward for κ > 0 and downward for κ < 0. For instance, we have V ef f ≈ 0.53 and 0.21 at 4πκ = 2 and 4πκ = −2, respectively; this is illustrated in figure 1b. In general, the positive shift of the V ef f with increasing values of κp φ > 0 reinforces the harmonic term in the potential and they together act to decrease the Lyapunov exponent, while κp φ < 0 gives a window of negative values (−5 < 4πκ < 0), in which we observe a slight increase in the Lyapunov spectrum as will be made manifestly in the next section.
It is also useful to have the contour plots of the V ef f at p φ = 0, 1, 2 for various values of κ as will refer to them in the next section. These are given in figure 1a, 1c and 1d. Sharp edges in these potential contours near θ ≈ 0 correspond to the flat direction of the pure YM potential. In the present case, the CS term helps to lift this, as the harmonic term in V ef f assists to shrink the sharp edges for all values of p φ and also acts to pull the contours toward closed loops for p φ = 0. For p φ > 0, the latter happens faster for κ > 0 as opposed to κ < 0 and vice versa for p φ < 0.

Analysis of the Chaotic Dynamics
We now explore the chaotic structure of the system governed by (2.21) by studying the Lyapunov spectrum and the Poincaré sections.

Lyapunov Spectrum
Setting the energy E = 1, = 0.1, and letting p φ assume the values 0, 1, 2, which is convenient for ease in comparison with the pure YM matrix model results in [5], we obtain the Largest Lyapunov Exponents (LLE), λ L , as the CS coupling takes on a range of values, in which the typical behavior of the LLE's are captured. Our results are obtained after averaging over 120 randomly selected initial conditions 6 in each case. They are presented in the figures 2a, 2b, 2c and we will elaborate on them shortly.
Chaotic structure of the pure YM model is explored in [5] and it is found that the system is fully chaotic at p φ = 0 and essentially becomes non-chaotic with increasing values of p φ . At the intermediate values 0 < p φ < 2, for example at p φ = 1, there are regions in the phase space, in which quasi-periodic motion is present as signaled by KAM tori appearing in the Poincaré section plots given in [5], while the rest of the phase space is filled with chaotic motion.
In figure 2a, profile of the Lyapunov spectrum of the model values of κ in the interval 4π|κ| ≤ 15 and at p φ = 0 is presented. The plot is essentially symmetric w.r.t. the κ = 0-axis as may be expected from (2.21), which is even under κ ↔ −κ for p φ = 0 and although LLE values tend to decrease in an almost monotonic manner for 4π|κ| > 4, they are essentially non-vanishing for 4π|κ| ≤ 10, which makes us conclude that the model is chaotic and behaves similar to the pure YM case within this range of the CS coupling. The rather mild increase in the LLE values observed in this plot in the narrow range 4π|κ| < 4 can be explained as follows. As |κ| increases, sharp edged regions in the contour plot of the effective potential V ef f , as illustrated in figure 1a, become less pronounced, and consequently, compared to κ = 0, the system spends relatively less time in these regions where the dynamics is adiabatic in θ and therefore no appreciable contribution to chaos is expected [5]. In turn, this results in a mild increase in the LLE spectrum within the indicated range of κ values. Nevertheless, for 4π|κ| > 4, harmonic term starts to become significant and the chaotic dynamics is gradually lost.
For p φ = 0, the κp φ term in V ef f impacts the Lyapunov spectrum asymmetrically depending on its sign, as it causes a fixed negative or a positive shift on the latter. At p φ = 1, for instance, which is illustrated in figure 2b, we immediately observe that λ L values with in the range of values −5 < 4πκ < 0 are above what is computed at κ = 0. This can be attributed to the downward shift in V ef f due to κp φ < 0, which clearly also lowers the minimum of V ef f as we have already discussed toward the end of previous section. The increase in λ L can not be sustained for 4πκ < −5, since then the harmonic term ∝ κ 2 r 2 becomes sufficiently strong even at short distances to dominate V ef f and initiates the decline of the chaotic dynamics. For κ > 0, this term acts to strengthen the harmonic terms and the chaotic motion becomes sharply suppressed before 4πκ ≈ 5. At p φ = 2, we still observe a mild increase in the Lyapunov exponents roughly in the range −4 < 4πκ < 0, but the maximum value of λ L now appears to be ≈ 0.03, an order of magnitude less than that is found for p φ = 0 and p φ = 1, and not significant enough to conclude that any dense chaotic dynamics remain for p φ ≥ 2.

Poincaré sections
All of the conclusions of the previous subsection regarding the chaotic dynamics of the present dynamical system are well supported by the Poincaré sections. We have obtained the latter at the θ = 0 intersections of the phase space and projected on to the p θ , p r plane. Figures 3, 4, 5 and6 show the Poincaré sections on the first quadrant of the p θ , p r plane.
From figure 3, we see that chaotic dynamics appears to fill the phase space at p φ = 0, for a large range of values of κ, which is approximately 4π|κ| 10, while the periodic motion starts to compete and take over after this range of κ values as can be observed from figure 3c. At p φ = 1 and 4πκ = 1, from figure 4, we observe that the phase space is still dominated by chaos, while a few KAM tori indicating quasi-periodic motion are visible. As κ continues to increase, more KAM tori start to occur, and system swiftly becomes non-chaotic for 4πκ 4 and gets dominated by quasi-periodic orbits. However, for κ < 0 as illustrated in figure  5, the system appears to remain densely chaotic with only a few KAM tori appearing until around 4πκ ≈ −5, while the quasi-periodic motion starts to spread for 4πκ −7.5 and start to take over only after 4πκ −10. Let us also note that, some KAM tori appear to intersect, especially as seen in the Poincaré sections at larger values of |κ|, for instance, in the figures (4d) and (5f). This is due to possible different values of the r coordinate appearing in the evolution of the system starting with distinct initial conditions being projected to the same point on the p θ ,p r plane. For p φ = 2, we see that there is very little chaos remaining in the phase space regardless of the value of κ and quasi periodic motion dominates the phase space. This can be seen from the Poincaré sections in figure 6. There is no chaos for κ > 0, and although some randomly spread points appear for negative κ values, for small |κ|, KAM tori quickly dominate the phase space and quasi periodic motion is all that is left.

Transition from Chaotic to Non-Chaotic Phase
In order to investigate the transition of the system from the chaotic phase, i.e. black-brane phase, to the non-chaotic, integrable phase dominated by quasi periodic motion, it is useful to examine, the change of λ L with p φ treated as the order parameter, while keeping κ fixed. The fitting curves presented in figure 7 helps to illustrate the situation with sufficient clarity. In particular, from this figure, we see that for −7 4πκ < 0, λ L decreases linearly with increasing p φ , and the transition between the two-phases occur at p c φ ≈ 2 with a critical exponent of 1 for λ L . Therefore, we conclude that approximately within this range of κ, the transition from the chaotic to non-chaotic phase has the same characteristic as found for the pure YM model in [5]. For 4πκ < −7, λ L value is already below 0.1 (see figure 2b) at p φ = 1 and tends to decrease faster with increasing p φ ; at 4πκ = −10 we estimate that what little remains of the chaotic phase approaches p c φ ≈ 2 with a critical exponent ≈ 3/2. For κ > 0, on the other hand, not only the approach to non-chaotic phase appears to be faster, but also it tends to occur at smaller values of p φ at larger κ. For instance, we estimate that at 4πκ = 1, p c φ ≈ 1.75, while at 4πκ = 2, p c φ ≈ 1.55, with critical exponents ≈ 3/2 and ≈ 5/2, respectively.
(a) 4πκ = 1  Estimates for p c φ and the critical exponents for λ L are obtained from the best fitting curves to the data.

Conclusions and Outlook
In this paper, we have studied the chaotic structure of the minimal Yang Mills Chern Simons matrix model. Using the gauge and global symmetries, and with a suitable choice of the coordinates, the Hamiltonian of the system is obtained in a form in which the Gauss law constraint is fully solved and manifestly imposed. We have studied the chaotic dynamics of the model, and in particular, probed the changes in the Lyapunov exponent as the values of both the CS coupling, κ, and the conserved conjugate momentum, p φ , are varied. We have found that, even for p φ = 0, there is a range of CS coupling values, approximately given as 4π|κ| 4 within which the Lyapunov exponent is larger in value compared to that evaluated at κ = 0. We have also seen that κp φ term in the effective potential alters the Lyapunov spectrum depending on its sign. We have found that the largest Lyapunov exponents evaluated within a range of values of κ are above that is computed at κ = 0, for κp φ < 0. These results are discussed in detail in section 3, where we also presented estimates for the critical exponents for λ L and the value of the order parameter, p c φ , as the system transits from chaotic to non-chaotic phase.
Let us finally note that out of time order correlators (OTOCs) approach recently applied in [35] to a system involving a x 2 y 2 term in the potential to probe quantum chaos may also be suitable for the model treated in this paper and we hope to report on any developments along these directions elsewhere in the near future.

A. Chern-Simons Action in 0 + 1-Dimensions
In 2 + 1-dimensions the CS action for the Hermitian SU (2) gauge fields A µ may be given as [39], [40] Invariance of e iS CS (2+1) under large gauge transformations requires that k is an integer 7 . Dimensional reduction from 2 + 1 to 0 + 1-dimension is facilitated by requiring that all spatial dependence of the gauge fields A µ (x, t) are dropped and only their time dependence is retained. This implies that in S CS (2+1) all the spatial derivatives collapse to zero. Introducing the notation A µ ≡ (A 0 , X i ), i : 1, 2, we immediately find where V 2 := d 2 x denotes the result of the 2-dimensional volume integral. Thus the coupling of the CS action in 0+1 dimensions take the form κ := kV 2 4π and due to the volume factor V 2 it differs from the coupling of the CS action in 2+1-dimensions. In particular, due to the V 2 factor, which can take arbitrary real values, κ is not an integer multiple of 1 4π . We also see that the CS coupling κ obtained in this manner is consistent with the fact that S CS (0+1) is manifestly gauge invariant (as it is trivially observed form the second line of (A.2) upon using the cyclicity of the trace and U † U = U U † = 1). These facts clearly indicate that κ is not level quantized.
In this paper we work with 4π|κ| 10 as this conveniently gives the relevant range of κ values to explore the dynamics at E = 1 and |p φ | 2. The latter are the values of energy and p φ used in [5], which we use since it gives us an ease in comparison. In our study, we are only considering the classical theory but writing out the explicitly in the action there will not be any reason to keep 4πκ within the values of O(1). 7 Let us also note that for gauge fields valued in the Lie algebra of the Lie group G, CS action on a 2D + 1 dimensional space-time manifold M changes under G transformation by a term of the form Ω ∝ M T r[(g −1 dg) 2D+1 ] (up to a constant factor) where g ∈ G. In general Ω does not vanish [40]. For M = S 2D+1 , Ω is determined via the (2D + 1) th homotopy group of G, π 2D+1 (G), and the requirement of gauge invariance of e iS CS leads to the level quantization of the CS coupling. If M is the 2D + 1 Minkowski space, Ω is still determined by π 2D+1 (G) if we demand that g → 1 as time goes to ±∞ and at spatial boundaries [39,40]. For D = 0 and G ≡ SU (2), we clearly have Ω = 0 since g −1 dg is valued in the Lie algebra of SU (2) and hence its trace is vanishing. Thus, this general consideration is applicable to the S CS in 0 + 1-dimensions and indicates that S CS (0+1) is indeed gauge invariant. In particular, we also have that π 1 (SU (2)) ≡ Π 1 (S 3 ) = 0 therefore Ω vanishes in 0 + 1-dimensions implying once again that S CS (0+1) is gauge invariant.

B.1. Calculation of Lyapunov Exponents
Lyapunov exponents are useful to determine the sensitivity of a system to given initial conditions. More precisely, they measure the exponential growth in perturbations and therefore give a reliable way to establish the presence of chaos in a dynamical system [36][37][38]. For a Hamiltonian system if we denote the perturbations in the phase space coordinates g(t) ≡ (g 1 (t), g 2 (t) , · · · , g 2N (t)) by δg(t), then we may conclude that the system is chaotic if, at large t, δg(t) deviates exponentially from its initial value at t = t 0 : ||δg(t)|| = e λ(t−t 0 ) ||δg(t 0 )||. Here λ > 0 are called the positive Lyapunov exponents and there are 2N of them for a phase space of dimension 2N . Let us also note that this description is in parallel with the statement that even slightly different initial conditions give trajectories in the phase space, which are exponentially diverging from each other and hence lead to chaos. In a dynamical system presence of at least one positive Lyapunov exponent is sufficient to conclude the presence of chaotic motion. In Hamiltonian systems, due to the symplectic structure of the phase space, Lyapunov exponents appear in λ i and −λ i pairs and a pair of the Lyapunov exponents vanishes as there is no exponential growth in perturbations along the direction of the trajectory specified by the initial condition. and sum of all the Lyapunov exponents is zero as a consequence of Liouville's theorem. These facts are well-known and their details may be found in many of the excellent books on chaos [36][37][38].
We follow the appendix in [25] to describe the method to compute all the Lyapunov exponents.With U (t) denoting a time evolution operator we may write and Lyapunov exponents are defined by Dividing the time into n equal steps such that t = n∆t Lyapunov exponents can be expressed as We may consider that (h 1 0 , h 2 0 , · · · , h 2N 0 ) span an orthonormal basis for the set of vectors tangent to the phase space trajectory g(0) at t = 0. After a time ∆t, time evolved vectors can be written as k i 1 = U (∆t)h i 0 , where (k 1 1 , k 2 1 , · · · , k 2N 1 ) spans a basis of tangents vectors to the trajectory g(∆t). However, this basis of vectors need not be orthogonal. Using the Gram-Schmidt orthogonalization process we can obtain the orthogonal set from (k 1 1 , k 2 1 , · · · , k 2N 1 ), which we may denote as (h 1 1 ,h 2 1 , · · · ,h 2N 1 ). The expansion rate of the vectorh i 1 can be determined as since all h i 0 is normalized to 1 already. An orthonormal basis after time ∆t is therefore given This procedure defines the time evolution after one step of ∆t . After n∆t steps we may write the Lyapunov exponents as The set {λ 1 , ..., λ 2N } is called the Lyapunov spectrum and as a consequence of this construction λ 1 is the largest Lyapunov exponent, since h 1 i leads toward the direction in the phase space which is most sensitive to the initial conditions and therefore the expansion rate r 1 i has the largest value in this region of the phase space.
A Matlab code solving the Hamilton's equations of motion and evaluating the Lyapunov exponents according to the procedure outlined above is used for our numerical calculations. Let us note that since the phase space is 4-dimensional there are only four Lyapunov exponents, two of which are zero and the remaining two may be denoted as λ L and −λ L in view of our earlier remarks in this section.
We pick the initial conditions that are used both in the evaluation of the Lyapunov spectrum and the Poincaré sections as follows. Both p r and θ are initially taken to be equal to zero. Using the Hamiltonian (2.21), p θ can be expressed as where we have already set E = 1 and = 0.1. We determine the intervals of r values, which make the argument of the square root in (B.8) positive and restrict to the one in which r > 0. Initial value of r is chosen randomly from this interval and the initial value of p θ is then determined from (B.8).
We run the Matlab code evaluating the largest positive Lyapunov exponent for 120 randomly selected initial conditions according to this procedure at each value of κ and take their average to obtain the each data point. The error bars are obtained by computing the mean square variances. In the simulation, we take a time step of 0.5 and run the code from time 0 to 3000. The λ L obtained in this manner is recorded at several values of p φ and κ.

C. Generic form of the Hamiltonian and the Pure Chern-Simons Limit
Restoring the YM coupling we may write the YMCS action in the form Let us note that the dimensional analysis shows that g has the units [Length] − 3 2 , κ has the units [Length] 2 , while the fields X i and A 0 are of dimension [Length] −1 . In the paper we 4 3 κ to work in the units where g is set to unity.
With the YM coupling included explicitly using (2.7) we can write while the corresponding conjugate momenta and the Hamiltonian are given as We may obtain the pure CS limit by letting g → ∞. In this limit, we see from (C.3) that p 1 → −κ x 2 and p 2 → κ x 1 . Substituting this in the Hamiltonian (C.4) we find Therefore we see that in the pure CS limit the Hamiltonian vanishes and hence no dynamics or chaos remains. This result is consistent with the general considerations arising from the Chern-Simons theories, as it is known that the pure CS theory has no dynamics, but nontrivial dynamics emerges from coupling to dynamical matter fields, or by considering the CS action on a manifold with boundaries [39]. More generally, for a system with generalized coordinates q i and velocitiesq i , Lagrangian involving first order time derivatives have the generic form where g ij is the metric, f i is some function of the generalized coordinates i.e. f i ≡ f i (q j ) and V is a potential V ≡ V (q i ). Canonical momenta are evaluated as In terms of p i ,q i can be solved using the inverse metric in the forṁ The Hamiltonian then takes the form given in (2.13). In order to discuss the limit in which the Lagrangian consists of only first order time derivatives, we may first set g ij → 1 We may, therefore, write the generic form of the Hamiltonian as while the conjugate momenta take the form p i = 1 g 2 g ijqj + f i . Thus, in the limit with g → ∞ we have, p i → f i and upon substituting this in (C.10) we immediately see that H = 0.
It is somewhat more subtle to see the pure CS limit from the Hamiltonian given in the angular coordinates in (2.21). We will provide a concrete discussion of this limit in the next appendix.