Chaos from Massive Deformations of Yang-Mills Matrix Models

We focus on an $SU(N)$ Yang-Mills gauge theory in $0+1$-dimensions with the same matrix content as the bosonic part of the BFSS matrix model, but with mass deformation terms breaking the global $SO(9)$ symmetry of the latter to $SO(5) \times SO(3) \times {\mathbb Z}_2$. Introducing an ansatz configuration involving fuzzy four and two spheres with collective time dependence, we examine the chaotic dynamics in a family of effective Lagrangians obtained by tracing over the aforementioned ansatz configurations at the matrix levels $N = \frac{1}{6}(n+1)(n+2)(n+3)$, for $n=1,2,\cdots\,,7$. Through numerical work, we determine the Lyapunov spectrum and analyze how the largest Lyapunov exponents change as a function of the energy. Poincar\'e sections of these effective models are also given at selected energies that best reflects the dynamics prior to, at and after the onset of chaos.


Introduction
Recently, there has been intense interest in exploring, modeling and understanding the emergence of chaos from matrix models of Yang-Mills(YM) theories [1][2][3][4][5][6][7][8][9][10][11][12]. Attention has been focused on examining the many aspects and features of chaos in the matrix models of primary interest within the context of M-theory and string theories, namely in the BFSS and BMN models [13][14][15][16][17]. These models are supersymmetric SU (N ) gauge theories in 0 + 1dimensions, whose bosonic part consist of nine N × N matrices and commonly referred to as matrix quantum mechanics in the literature. BFSS model is associated to the type II-A string theory and appears as the DLCQ of the M-theory in the flat background, while the BMN model, being a deformation of the BFSS involving quadratic and qubic terms preserving the maximal amount of supersymmetry, is the DLCQ of M-theory on the pp-wave background. They describe the dynamics of N -coincident D0-branes in flat and spherical backgrounds, respectively; the latter being due to the fact that fuzzy 2-spheres are vacuum configurations in the BMN model. In the gravity dual, i.e. at large N and strong coupling/low temperature limit, it is known that these D0-branes form a black-brane phase, i.e. a string theoretic black hole [16][17][18].
From these considerations it is apparent that a strong motivation in studying possible chaotic dynamics possessed by these matrix models comes from possibilities in acquiring an in depth perspective and perhaps even discovering novel physical phenomena relating to the dual black-brane phases, such as their thermalization, evaporation processes. In some of the recent works, large N and high temperature limit of the BFSS and BMN models are probed using numerical methods and results indicating fast termalization are encountered [3] and interpreted to indicate fast scrambling [1]. The latter is a conjecture on black holes, proposed by Sekino and Susskind, which may be concisely stated as the fact that chaotic dynamics sets in faster in black holes than in any other physical system and the rate at which this happens is proportional to logarithm of the number of degrees of freedom, i.e. to the logarithm of the Bekenstein-Hawking entropy of the black hole. It may be remarked that although this large N , high temperature limit, is distinguished from the limit in which the gravity duals are obtained, as noted in [5], in numerical studies conducted so far [2],no indication of an occurrence of a phase transition between low and high temperature limits of these matrix models is reported, which makes it plausible to expect that features like fast scrambling of black holes in the gravity dual could survive at the high temperature limit too. In fact, the latter is the natural limit in which classical dynamics provides a good approximation to the matrix quantum mechanics, 1 which is free from fermions as the latter could alter or impact the dynamics of the bosonic matrices only at low temperatures. Analysis performed in [5] for the BFSS model by calculating the Lyapunov spectrum demonstrated a classical version of the fast scrambling holds with scrambling time being proportional to log N 2 .
Authors of [9] considered different ansatz configurations in the BMN model at small matrix sizes of N = 2, 3 to probe chaos. In particular, they have focused on configurations of matrices which are essentially fuzzy two spheres at angular momentum 1/2 and 1, but with a collective dependence to distinct real time dependent functions such that the Gauss law constraint is readily satisfied. The effective models obtained in this manner exhibit chaos as demonstrated in [9]. From a general perspective, appearance of chaos in BFSS and BMN models is not surprising, but also consistent with the early investigations on chaotic motion in Yang-Mills theories in 80's [19][20][21] and the initial work performed on the BFSS model in [22]. Let us also note that, BFSS as well as BMN models and their smaller subsectors at small values of N are highly non-trivial many-body system whose complete solution evades us to date. Studies made so far provide some qualitative implications on the black hole phases in such models; for instance, in an SU (2) YM matrix model with only two 2 × 2 matrices inspected in [6], the edge of the chaotic region is argued to correspond to the end of the black hole phase. In a recent paper [12] coauthored by one of the present authors, a Yang Mills five matrix models with a mass term is examined in detail and by considering equivariant fluctuations around its fuzzy four sphere [17,23,24] solutions an effective action is obtained. The chaotic dynamics in this model is revealed through detailed numerical analysis. Since non-trivial D0-brane dynamics leading to chaos require non-commutativity [8], fuzzy spheres which are described by matrices with certain non-vanishing commutation properties could serve as good candidates for background geometries to probe chaotic behaviour in YM matrix models. Dynamics of D0 − D2-and D0 − D4-brane systems with spherical D2-and D4-branes, i.e. on the fuzzy two-and four-spheres, are studied in the literature from various perspectives [23,[25][26][27][28] in the past, while studies oriented toward exploring chaos in this context appears to be limited in number.
We therefore think that it will be quite useful to further investigate the role played by fuzzy two-and four-sphere configurations in the development of chaotic dynamics in Yang Mills matrix models, in a setting which does not confine us to only small matrix sizes but rather allows us expand from small to reasonably larger size matrices and that could, in turn, yield more comprehensive data through which we may be able to comment further on the dynamics of D0-branes in the background of fuzzy spheres. In order to serve these purposes, we think that it may be useful to consider a matrix model in which mass deformation terms allow for both fuzzy two and four spheres. In particular, we take a Yang Mills matrix model with the same content as that of the BFSS, but with mass deformation terms breaking the global SO(9) symmetry of the latter to SO(5) × SO(3) × Z 2 . Introducing a configuration involving fuzzy four and two spheres with collective time dependence as an ansatz, we study the chaotic dynamics in a family of effective Lagrangians which we obtain by tracing over these configurations at the matrix levels N = 1 6 (n+1)(n+2)(n+3), for n = 1, 2, · · · , 7, corresponding to the sizes of fuzzy four spheres. Before presenting a full numerical analysis of the chaotic dynamics, we first identify the fixed points of the phase space for the model at each matrix level and perform a linear stability analysis of these points. The latter allows us to identify both the stable and unstable fixed points, an it also suggests that the possible chaotic motion in these models could start at energies around and above the lowest unstable fixed point energy within a given model. This expectation is corroborated through the numerical work we perform, in which we obtain the Lyapunov spectrum and analyze how the largest Lyapunov exponents change as a function of the energy. We also give the Poincaré sections of these effective models at selected energies that best describes the dynamics prior to, at and after, the onset of chaos. The fuzzy two spheres in these configurations have maximal spin N −1 2 at a given matrix level N . We briefly inspect the ansatz, where the other extreme, namely a direct sum of minimal spin 1 2 fuzzy spheres are present. In this latter case, the models develop chaotic dynamics only within a narrow interval of energies. We speculate on how these two results may be intertwined from a broader perspective before closing the paper.

Yang-Mills Matrix Models with Double Mass Deformation
In this section we introduce a Yang-Mills (YM) matrix models with two distinct mass deformation terms, which may be contemplated as a double mass deformation of the bosonic part of the BFSS model. The latter can be described by the action T r stands for the normalized trace, T r1 N = 1. The covariant derivatives are given by

3)
where A t is the U (N ) gauge field transforming as (2.7) is an immediate consequence of the algebraic equations of motion for the gauge field A t and amounts to having physical states as SU (N )-singlets in the quantized theory. A gauge invariant double mass deformation of (2.1) may be specified as where the indices a and i take on the values a = 1, .., 5 and i = 6, 7, 8, respectively. In (2.8) the terms proportional to µ 2 1 and µ 2 2 are the quadratic deformations, which respect the U (N ) gauge symmetry, but altogether break the SO(9) down to SO(5) × SO(3) × Z 2 . The discrete Z 2 factor is present, simply because, we have not introduced any mass terms for the matrix X 9 : in fact,in what follows, we are going to consider the sector in which X 9 is set equal to the zero matrix, or equally well, taken as the zero matrix as a part of the ansatz configuration which will be introduced shortly. Since, we are going to be essentially concerned with the classical dynamics of (2.8), we scale the coupling constant to unit value from now on, as it only determines the overall scale of energy classically.
In the A t = 0, gauge the equations of motion for X I take the form while the Gauss law constraint remains unchanged in the form as given in (2.7). The massive deformation of the BFSS model, which preserves maximal amount of supersymmetry is already known to be the BMN matrix model [14], which possesses fuzzy twospheres and their direct sums as possible vacuum configurations. However, in the present work, our focus is directed toward exploring the emerging chaotic dynamics from the Yang-Mills matrix models which could allow for not only fuzzy two-sphere configurations but also higher dimensional fuzzy spheres, and in particular a fuzzy four-sphere. It is possible to conceive deformations of S BF SS including two separate mass terms, which break the SO(9) symmetry down to several different product subgroups. The underlying motivation for introducing the specific massive deformation in (2.8) comes from the fact that in two distinct limiting cases the equations of motion can be solved either with fuzzy two-sphere our fuzzy four-sphere configurations. To be more precise, we have for X i = 0 = X 9 , (2.9b) and (2.9c) satisfied identically, while (2.9a) takes the form which is satisfied by the fuzzy four sphere configurations X a ≡ Y a for µ 2 1 = −16. Y a are N × N matrices carrying the (0, n) UIR of SO(5) with N = 1 6 (n + 1)(n + 2)(n + 3) and satisfying the defining properties given in (A.11). Whereas, in the other extreme, one may set X a = 0 = X 9 , with the only remaining non-trivial equation of motion which is solved by fuzzy two sphere configurations X i ≡ Z i or their direct sum for µ 2 2 = −2. In this case, X i are N × N matrices carrying the spin j = N −1 2 UIR of SO(3) ≈ SU (2). In view of these observations we consider an ansätz configurations involving fuzzy twoand four-spheres with collective time dependence, which fulfill the Gauss law constraint given in (2.7). Tracing over the fuzzy two-and four-sphere configurations, we aim to obtain reduced models with only four phase space degrees of freedom, whose dynamics can be investigated in considerable detail.

Ansatz I and the Effective Action
A reasonably simple, yet non-trivial configuration is constructed by introducing two separate collective time-dependent functions multiplying the fuzzy four-and two-sphere matrices.
where r(t) and y(t) are real functions of time. In this ansatz, we consider a single spinj = N −1 2 IRR of SU (2) as the fuzzy S 2 configuration, while taking direct sums of fuzzy two spheres with different IRRs of SU (2), i.e. forming Z i as a block-diagonal matrix composed of a direct sum of different IRRs of SU (2) with spin less than j remains an open possibility. We will briefly discuss Z i 's made up of spin 1 2 fuzzy spheres in section 4. Z i exist at every matrix level, while this is not so for Y a . Fuzzy four spheres exist at the matrix levels 4, 10, 20 · · · as given by the dimension N = 1 6 (n + 1)(n + 2)(n + 3) of the IRR (0, n) of SO (5). Accordingly the fuzzy two spheres are taken at the matrix levels matching these dimensions of the fuzzy four sphere. In what follows, we initially keep the mass parameters µ 2 1 and µ 2 2 as unspecified in some of the key equations, but will investigate the detailed dynamics for µ 2 1 = −16 and µ 2 2 = −2, which emerge as the limiting values for the static solutions of (2.10) and (2.11) and subsequently will also briefly investigate the consequences of taking another set of values for the masses, namely µ 2 1 = −8 and µ 2 2 = 1 on the chaotic dynamics of the reduced models. Substituting the (2.12) configuration in the action (2.8), we perform the trace over the fuzzy four-and two-sphere matrices at each of the matrix levels N = 1 6 (n + 1)(n + 2)(n + 3) for n = 1, 2, · · · , 6. Using Matlab to evaluate the traces we obtain the Lagrangian of the reduced models in the form where the coefficients c β = c β (n) (β = 1, 2, 3) depend on n and their values (given up to one digit after the decimal point at most) for n = 1, 2, · · · , 7 are listed in the table 1 given below. We suppress the label n of the coefficients c β (n) in (2.13) in order not to clutter the notation. Coefficients given for n = 7 in table 1 are evaluated by obtaining a polynomial function of n approximating 2 c β (n) for n = 1, 2, · · · , 6 and interpolating this result to n = 7.  Table 1 The corresponding Hamiltonian is easily obtained from (2.13) and it is given by where V n (r, y) denotes the potential function.
To explore the dynamics of the models governed by H n , we first evaluate the Hamilton's equations of motion. These take the forṁ Taking the mass parameter values as µ 2 In order to investigate the dynamics of these models governed by (2.16) in detail, it is quite useful to start the analysis by determining the fixed points corresponding to the equations of motion in (2.16) and addressing their stability at the linear order. Fixed points of a Hamiltonian system are defined as the stationary points of the phase space [30][31][32][33], and can thefore be determined by the set of equations specified as These polynomial functions are given as Using (2.17) in (2.16) leads to four algebraic equations, two of which are trivially solved by (p r , p y ) ≡ (0, 0), which means that all the fixed points are confined to the (p r , p y ) ≡ (0, 0) plane in the phase space. The remaining two equations are and have the general set of solutions given as where h 1 and h 2 are given in terms of c β as (2.20) Clearly, only real solutions of (2.18) are physically acceptable. From table 1 it is straightforward to compute that both h 1 and h 2 are real except at n = 1. For n > 1 the set of fixed points are given as where the values of h 1 and h 2 are presented in the table 2 below  Table 2 while for the n = 1 model we only have as the fixed points. Let us note that (2.19) corresponds to the critical points of the potential V n , since (2.18) are the equations determining the extrema of the latter. From the eigenvalues of the matrix ∂ 2 V n ∂g i ∂g j , (with the help of the notation (g 1 , g 2 ) ≡ (r, y)), we see that the points (±1, 0) and (0 , ±1) are local minima, (0, 0) is a local maximum, while (±h 1 (n), ±h 2 (n)), (±h 1 (n), ∓h 2 (n)) are all saddle points of V n . Evaluating V n at the local minima, we find V n (±1, 0) = −8c 1 and V n (0, ±1) = −c 2 . Using table 1, we easily conclude that (±1, 0) is the absolute minimum of V n for n = 1, 2, 3, while (0, ±1) is the absolute minima for the models with n = 4, 5, 6. From now on we add to the Hamiltonian's, H n , the constant term 8c 1 for n = 1, 2, 3, and c 2 for n = 4, 5, 6, respectively, to shift the minimum value of the potentials V n to zero and use the notation (2.23)  Table 3 2

.2. Linear Stability Analysis in the Phase Space
It must be clear that the properties of extrema of V n does not provide sufficient information to decide on the stability of the fixed points. We now perform a first order stability analysis around the fixed of H n given in (2.21) and (2.22). Together with the Lyapunov spectrum and the Poincare sections that will be determined in the next section, this analysis, will allow us to comment on the outset and variation of chaos, that is, the increase and decrease in the amount of chaotic orbits in the phase spaces of H n , w.r.t. energy. For the phase space coordinates it is useful to introduce the notation (g 1 , g 2 , g 3 , g 4 ) ≡ (r, y, p r , p y ) ,

25)
From g α andġ α , we may form the Jacobian matrix

26)
whose eigenvalue structure allows us to decide on the stability character of the fixed points [30][31][32][33]. Written in explicit form, we have

27)
where J 31 and J 42 are (2.28) Eigenvalues of J(r, y) at the fixed points (2.21) are easily evaluated and listed in the table 4 below.

Fixed Points
Eigenvalues of J(r, y)

Table 4
General criterion of the linear stability analysis [30][31][32][33] states that a fixed point is stable if all the real eigenvalues of the Jacobian are negative, and unstable if the Jacobian has at least one real positive eigenvalue. It may that all the eigenvalues of the Jacobian matrix are imaginary. This is called the borderline case and an analysis beyond first order is necessary to decide if the system is stable or unstable at such a point. Accordingly, we see that (0, 0, 0, 0) and (h 1 (n), h 2 (n), 0, 0) are all unstable fixed points as the corresponding Jacobians have at least one real positive eigenvalue, as readily seen from the table 4. From the same table and the values of c α (n) in table 1, we see that (±1, 0, 0, 0) and (0, ±1, 0, 0) are borderline cases. We are not going to explore the structure of these borderline fixed points any further, as we expect that their impact on the chaotic dynamics should be rather small compared to those of the unstable fixed points which we just identified at the linear level. Our numerical results on the Lyapunov spectrum indeed corroborates with this expectation as will be discussed shortly.

Lyapunov Spectrum
In order to probe the presence and analyze the structure of chaotic dynamics of the models described by the Hamiltonians H n , we will first examine their Lyapunov spectrum. As it is well known, Lyapunov exponents measure the exponential growth in perturbations and therefore give a reliable and decisive means to establish the presence of chaos in a dynamical system [31][32][33]. Suppose that we denote the perturbations in phase space coordinates g(t) ≡ (g 1 (t), g 2 (t) , · · · , g 2n (t)) by δg(t). The system is chaotic if δg(t) deviates exponentially from its initial value at t = t 0 : |δg(t)| = e λt |δg(t)|, and λ > 0 are called the Lyapunov exponents.
There are 2n of them for a phase space of dimension 2n. This description is essentially equivalent to 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. Thus, 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, 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 the Liouville's theorem. These facts are well-known and their details may be found in many of the excellent books on chaos [31][32][33].
In order to obtain the Lyapunov spectrum for our models we run a Matlab code, which numerically solves the Hamilton's equations of motion in (2.16) for all H n (1 ≤ n ≤ 7) at several different values of the energy. We run the code 40 times with randomly selected initial conditions satisfying a given energy and calculate the mean of the time series from all runs for each of the Lyapunov exponents at each value of n. In order to give certain effectiveness to the random initial condition selection process we developed a simple approach which we briefly explain next. Let us denote a generic set of initial conditions at t = 0 by (r(0), y(0), p r (0), p y (0)). For H n<4 we take y(0) = 0 and for H n≥4 , r(0) = 0 as part of the initial condition and subsequently generate three random numbers ω i (i = 1, 2, 3) and E for a given energy E of the system, so that E = Ω 2 i = Ω 2 1 + Ω 2 2 + Ω 2 3 . Subsequently, we take positive roots in the expressions 3 p r (0) = 4c 1 Ω 2 1 , p y (0) = 4c 2 Ω 2 2 , (3.1) and the real roots of

2)
where in the last step of the process our code randomly selects from the available real roots of the equations (3.2). In the computations we use a time step of 0.25 and run the code for a sufficient amount of computer time to clearly observe the values that the Lyapunov exponents converge to. Below, we present sample plots for the time series of Lyapunov exponents at each value of n. From the figures 1 chaotic dynamics of the models are clearly observed, as in each case (except in figure 1b) a positive Lyapunov exponent is present. We also observe that the properties of Lyapunov spectrum for Hamiltonian systems summarized at the end of the first paragraph of this section are readily satisfied. Let us immediately note that for the model at n = 1 have distinct features from the rest. This is already observed from the first two plots (figures 1a and 1b); for E = 30 there is a positive Lyapunov exponent, while at E = 500 all the Lyapunov exponents appear to be converging to zero indicating that very little chaos remains at this energy. The distinct features of the n = 1 model will also be seen in the ensuing discussions.  In order to see the dependence of the mean largest Lyapunov exponent(MLLE), (denoted as λ n in figures and tables) to energy, we obtain the MLLE at several different values of energy in a range, which appears to be best suited to observe the onset and progression of chaotic dynamics in these models. As may be expected, the energies determined for the unstable fixed points in the previous section are of central importance here. From the figures 2 we see that appreciable amount of chaotic dynamics starts to develop once the energy of the systems exceeds the energy E F of the models at the fixed points (±h 1 (n), ±(∓)h 2 (n), 0, 0). The onset of the chaotic dynamics as observed from progression of the MLLE values with increasing energy is highlighted by the blow-up figures provided in the insets of the plots in 2. Error bars at each data point is found by evaluating mean square error using MLLE and the LLE values of each of the 40 runs.  Several observations can be made from these numerical results. Firstly, we see that in the n = 1 model the MLLE acquires a peak value of about ≈ 0.55 at an energy around ≈ 32 and rapidly decreases toward zero with increasing energy. From the profile of MLLE w.r.t. energy in figure 2a as well as the time series plot (1b) of the model at E = 500, we conclude that this model is not chaotic for energies E ≥ 500. These conclusions are also fully supported by the Poincaré sections given in the figures 3a -3d.
In order to elaborate on the data obtained for the MLLE values, it is useful to explore the dependence of the LLEs at a given level n with respect to the energy. We find that the function λ n (E) = α n + β n 1 √ E , gives not perfect but essentially very good fits to our data as can be seen from the figures 2.  ing these fits, we find that the energies at which the MLLE vanish are given approximately as 25, 58, 122, 352, 799, 1760 for n = 2, 3, 4, 5, 6, 7, respectively. We find that these are somewhat less than the E F 's given in table 3. As previously noted, the latter are marking the onset of chaotic dynamics in our models and the comparatively lower values of energy found for vanishing MLLE is consistent with this fact. Results of this extrapolation appear to be also consistent with the Poincaré sections obtained at nearby energies as can be seen from the figures 3. We may use (3.3) to compare the relative rate at which the chaotic dynamics tends to develop once the systems reach their respective fixed point energies in table 3. Defining as the quantity measuring this relative rate, we find that R n takes on the values 0.017, 0.015, 0.010, 0.005, 0.002, 0.0012 for n = 2, 3, 4, 5, 6, 7, respectively. Thus, as n increases, R n values indicate a slow decrease in the rate at which MLLE increases with energy. This suggests that the models at low values of n become chaotic somewhat more rapidly. It is also interesting to note that MLLE values show essentially the same functional relationship with the energy, as that was found in the Yang-Mills 5-matrix models with a mass term, which was studied in [12]. Let us also remark on a feature of the models H n≥4 which is observed from the plots in figure 2d -2g at energies E 2500. We see that for a range of energy values in these models there is an observable decrease in the value of MLLE and the mean square errors appear to be considerably larger than those computed for the rest of the data points. A closer analysis of the Lyapunov time series at these data points reveal that LLE values of less than a quarter of the 40 initial conditions approach to zero, leading to the observed decrease in MLLE values and the increase in the mean square errors. From a physical point of view, approach of some the LLE values to zero implies that the systems' development in time, starting from these initial conditions are of either periodic or quasi-periodic type and not chaotic. Nevertheless, the overall MLLE values are still quite large and the sample Poincaré sections taken at one of these energies are densely chaotic, showing no sign of KAM tori signaling the presence of quasi-periodic orbits. Thefore, we are inclined to think that such periodic or quasi-periodic orbits occur only at comparatively very small regions of the phase space and evaluated the MLLE values at these energies by excluding the initial conditions leading to vanishing LLE. The results of MLLE obtained this way are given in the plots in the figures 2d -2g in yellow color for comparison and used to obtain the fits given in the table 5.

Poincaré Sections
In order to supplement the results of the previous section, we have obtained the Poincaré sections at several different values of the energy. For each model we plot the Poincaré sections at energies below, around and above the energy of unstable fixed points (±h 1 (n), ±(∓)h 2 (n), 0, 0) to visualize how the phase space trajectories develop and capture the onset of chaotic dynamics. As in the calculation of Lyapunov spectrum we use 40 randomly selected initial conditions for each model at a given energy and the same procedure as in the analysis of the Lyapunov spectrum is followed. For H n<4 , y(0) = 0 is a part of the initial conditions and for H n≥4 , r(0) = 0 is so, therefore it is convenient to look at the Poincaré sections on the r − p r -plane and the y − p y -plane, respectively in these cases. Our plots for the models n = 1, 2, 4 are presented below in the figure 3, while the results for the n = 3, 5, 6, 7 are carried over to the appendix and given in the figures 7.

Other Mass Values
We now consider assigning different values to the mass parameters and their impact on the dynamics of our models. Let us first note that µ 2 1 = −16 and µ 2 2 = −2 is a suitable and immediate guiding choice for the mass values due to the reasons discussed around equations (2.10) and (2.11), but surely not canonical in the sense that they are not enforced on us by the full set of equations of motion (2.9). Considering, Chern-Simons-type terms, i.e. a cubic term in X i 's as in the BMN model and a fifth order term in X a 's [?] alter the equations of motion given in (2.10) and (2.11) and lead to solutions with µ 2 1 > −16 and µ 2 2 > −2. It is not our aim in this section to provide a detailed analysis of such possibilities, but simply confine ourselves to examining another choice for the mass values and to serve this purpose we take µ 2 1 = −8 and µ 2 2 = 1. With the Lagrangian and Hamiltonian given in the form (2.13) and (2.14) and the corresponding Hamilton's equations given as in (2.15), we find that the fixed points of the phase space are  Linear stability analysis shows that (0, 0, 0, 0) is an unstable fixed point while (± 1 √ 2 , 0, 0, 0) are of the borderline type that we encountered previously and play no significant role in the numerical analysis that follows next. For convenience of comparison to numerical data for n = 1 , · · · 6, the fixed point energies are given in the same order as 12, 21, 32, 45, 60, 77.
Following the same steps of the numerical analysis for the Lyapunov spectrum as outlined previously, we find that, in this case too, the models exhibit chaotic dynamics for n = 2, 3, 4, 5, 6, while the model at the level n = 1 is essentially not chaotic for energies E 100, but retain some chaos only in a narrow band of energy from around E ≈ 12 (i.e. the fixed point energy) to E 100 (See the figures 4). The transition to chaos for n = 2, 3, 4, 5, 6 appears to happen around the fixed point energies as can be clearly seen from the blow-up insets of these figures. This fact is also captured by inspecting the Poincaré sections taken at energies somewhat below and above those of the fixed points.
In contrast to the previous case, MLLE values fluctuate and the mean square errors are comparatively larger at a narrow band of energies after E F . Due to this reason, for each model, we consider fits to the data starting at end of this transient band, where the change in MLLE w.r.t. energy starts to settle in a steady pattern of development and do not attempt to extrapolate them all the way to zero MLLE value as it would clearly be misleading to do so. We find that a logarithmic fit of the form λ n (E) =α n +β n log(E) ,  Table 6 appears to be well-suited to model the dependence of MLLE to energy after transient band is passed. Let us also note that within the transient band of energies, chaotic and quasi periodic motion coexist, this is clearly seen from those of the Poincaré sections in 9 that are taken at energies somewhat above the E F . Using (3.7) we havẽ measuring the relative rate of the development of chaos at E F . We find that R n takes the values 0.016, 0.011, 0.008, 0.009, 0.007 for n = 2, 3, 4, 5, 6, respectively, which essentially indicates that there is not any significant difference at the rate in which chaos develops in these models once the systems have energies around those of the fixed point.

Ansätz II
We would like to briefly discuss the consequences of changing the fuzzy two sphere part of ansätz I. Namely, we consider the configurations X a = r(t) Y a , X i = y(t) Z i , X 9 = 0 , (4.1) at the matrix levels N = 1 6 (n + 1)(n + 2)(n + 3) for n = 2, 3, 5. In (4.1), Y a are the same as in ansätz I, while we take Z i = ⊕ K n k=1 Σ i (k), with Σ i (k) spanning the spin 1 2 UIR of SU (2), in the k th block of the direct sum and K n is the number of 2 × 2 blocks, which are 5, 10, 23 for n = 2, 3, 5, respectively. Thus, we have a direct sum of 2 × 2 fuzzy two spheres. 5 We call this the ansätz II.
The Linear stability analysis reveals that (0, 0, 0, 0) is the only unstable fixed point in these models.
We have numerically studied the Lyapunov spectrum at several different values of the energy in these systems and determined that for an interval starting around the fixed point energies E F = 48, 84, 180, respectively and going up to ≈ 500 for the first two cases and ≈ 1000 for n = 5, there is a positive Lyapunov exponent indicating the presence of chaotic motion. From the times series plots in figure 5 and the Poincaré sections in figure 6 taken at energies within these ranges, we see that chaos is not dense, coexists together with quasi periodic motion and remains local in the phase space. At higher energies chaos ceases to exist and the phase space becomes dominated by periodic and/or quasi periodic orbits.

Conclusions and Outlook
We have studied the emergence of chaos in a massive deformation of a YM matrix theory, which has the same matrix content as that of the bosonic part of the BFSS model. Using an ansatz configuration involving fuzzy two-and four-spheres as backgrounds and assuming collective time dependence of the matrices, we were able to obtain a family of effective models descending from tracing over the fuzzy spheres at matrix levels N = 1 6 (n+1)(n+2)(n+3), for n = 1 , · · · , 7. We have performed a detailed numerical analysis and demonstrated the development of chaotic dynamics in these reduced models by obtaining their Lyapunov spectrum and Poincaré sections. From our results, were able to see that the onset of chaotic motion is at the energies which are at or around the lowest of that of the unstable fixed points and modeled, by fitting curves to the data, how the largest Lyapunov exponents change as a function of the energy for two different set of the mass values. The similarities and differences in these two cases are also discussed. In the penultimate section, we have considered another ansatz, which leads to chaos only with in a narrow interval of energies in each model.
As a final remark, we want to look at the results from ansatze I and II from a broader perspective. In the former case, fuzzy two sphere part is taken as the maximal one for the given matrix size, while in the latter, it is made up of blocks of lowest level spin 1/2 fuzzy spheres. Thus, we may conclude that fragmentation of the maximal fuzzy sphere into smallest fuzzy spheres results in the plummeting of the chaotic motion. Put it in another way, chaotic motion is favored in configurations with a big blob of a fuzzy two sphere than those in which small bubbles of fuzzy spheres proliferate. Although in both cases, Z i 's describing block diagonal and the maximal fuzzy sphere are both sparse, i.e. mostly filled with zeros and become even more so with increasing n, in the former case, interactions among D0branes in some directions remain confined to distinct pairs only, while for the maximal fuzzy sphere, they interact on a chain of nearest neighbors and this makes a difference. Viewing these facts in the reverse order suggests that as small bubbles of fuzzy spheres coalesce to form fuzzy spheres of higher spin, it also triggers the entire system to evolve from phase space filled with quasi periodic motion to that dominated by chaos.
A possible mechanism for splitting or fragmentation of fuzzy spheres into fuzzy spheres of smaller spin was suggested some time ago in [35] and it involves the introduction of an appropriate SU (2) equivariant co-product and an associated Hopf algebra structure. We suspect that it may be possible to make use of this co-product in the present context to get a handle on the topology change of the fuzzy background and its impact on the dynamics of D0-branes and therefore on chaos. We hope to report on any developments in these directions elsewhere.

A.1. Fuzzy S 2
Let us denote the SU (2) representations in the spin-j = N −1 2 irreducible representation by L i (i = 1, 2, 3). They satisfy the commutation relations The fuzzy two-sphere, S 2 F at the matrix level N is defined in terms of the rescaled generators of SU (2) and their commutation relations are simply given as

They satisfy
As N goes to infinity, the standard commutative S 2 is recovered.

A.2. Fuzzy S 4
Fuzzy four-sphere, S 4 F can be constructed as follows. We introduce the γ-matrices in five dimensions that are associated to the isometry group SO(5) of S 4 . They are 4x4 matrices satisfying the anticommutation relations γ a act on the four-dimensional spinor space C 4 which is the carrier space of the (0, 1) IRR of SO(5). Let us define the Hilbert space H n , as the carrier space of the (0, n) IRR of SO (5), which may be formed as the n-fold symmetric tensor product of C 4 as H n = C 4 ⊗ · · · ⊗ C 4 Sym .
(A. 8) and has the dimension N = dim(0, n) = 1 6 (n + 1)(n + 2)(n + 3) . (A.9) We may as well form the n-fold symmetric tensor product of the γ-matrices X a = 1 2 (γ a ⊗ 1 4 ⊗ ... ⊗ 1 4 + ... + 1 4 ⊗ ... ⊗ 1 4 ⊗ γ a ) (A.10) which carry the (0, n) IRR of SO(5) and naturally act on H n . X a are N × N Hermitian matrices, which satisfy the defining relations for the fuzzy four-sphere, which are given as X a X a = 1 4 n(n + 4)1 N , abcde X a X b X c X d = (n + 2)X e , (A.11) where X a X a is an SO(5) invariant operator.Commutators of X a yields [X a , X b ] = iM ab , (A. 12) where M ab are the generators of SO(5) in the (0, n) IRR and satisfy the same commutation relations as G ab . We also have that X a transform as vectors under the SO(5) rotations generated by M ab , i.e.
[M ab , X c ] = i(δ ac X b − δ bc X a ) . (A.13) From (A.12) we see that, as opposed to the construction of S 2 F , the algebra of X a 's defining the S 4 F do not close. A more detailed analysis yields that S 4 F may be understood as a squashed CP 3 F with S 2 F fibers In other words, we may state that S 4 F comes attached with a fuzzy two sphere at its every point. Further details of these features may be found in [16].

B. Poincaré Sections
Below we give the Poincaré sections taken for the models at the levels n = 3, 5, 6, 7 at various energies.