Chaotic Dynamics of the Mass Deformed ABJM Model

We explore the chaotic dynamics of the mass-deformed Aharony-Bergman-Jafferis-Maldacena model. To do so, we first perform a dimensional reduction of this model from $2+1$ to $0+1$ dimensions, considering that the fields are spatially uniform. Working in the 't Hooft limit and tracing over ansatz configurations involving fuzzy 2-spheres, which are described in terms of the Gomis-Rodriguez-Gomez-Van Raamsdonk-Verlinde matrices with collective time dependence, we obtain a family of reduced effective Lagrangians and demonstrate that they have chaotic dynamics by computing the associated Lyapunov exponents. In particular, we focus on how the largest Lyapunov exponent, $\lambda_L$, changes as a function of $E/N^2$. Depending on the structure of the effective potentials, we find either $\lambda_L \propto (E/N^2)^{1/3}$ or $\lambda_L \propto (E/N^2 - \gamma_N)^{1/3}$, where $\gamma_N(k, \mu)$ are constants determined in terms of the Chern-Simons coupling $k$, the mass $\mu$, and the matrix level $N$. Noting that the classical dynamics approximates the quantum theory only in the high-temperature regime, we investigate the temperature dependence of the largest Lyapunov exponents and give upper bounds on the temperature above which $\lambda_L$ values comply with the Maldacena-Shenker-Stanford bound, $ \lambda_L \leq 2 \pi T $, and below which it will eventually be not obeyed.


Introduction
Studies on exploring the structure of chaotic dynamics emerging from the matrix quantum mechanics have been continuing with growing interest for quite sometime [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Early investigations on the chaotic dynamics of Yang-Mills (YM) gauge theories dates back to the 1980's [15][16][17] and in the context of the Banks-Fischler-Shenker-Susskind (BFSS) model [18] to the work Arefeva et al. [19]. Recent studies are especially motivated by a result due to Maldacena-Shenker-Stanford (MSS) [6], which briefly states that, under rather general conditions met by a physical system, the largest Lyapunov exponent (which is a measure of chaos in both classical and quantum mechanical systems) for quantum chaotic dynamics is controlled by a temperature-dependent bound and given by λ L ≤ 2πT . It is conjectured that systems which are holographically dual to the black holes are maximally chaotic, meaning that they saturate this bound. This is already demonstrated for a particular fermionic matrix model, namely, the Sachdev-Ye-Kitaev [7] model and expected to be so for other matrix models which have a holographic dual such as the BFSS [18] model. The latter and the Berenstein-Maldacena-Nastase (BMN) model [20] are supersymmetric SU (N ) gauge theories, describing the dynamics of the N -coincident D0-branes in the flat and spherical backgrounds, respectively, and also appear in the Discrete Light Cone Quantization (DLCQ) of M theory in the flat and the pp-wave backgrounds [18,[20][21][22][23][24][25]. It is well known that the gravity dual of the BFSS model is obtained in the 't Hooft limit and describes a phase in which D0-branes form a so-called black brane, i.e., a string theoretical black hole [24][25][26].
Classical dynamics of YM matrix models provide a good approximation of the hightemperature limit of the quantum theory. Although this regime is distinguished from that in which the gravity dual is obtained (i.e., the low-temperature limit), early numerical studies conducted in Refs. [27,28] gave no indication of an occurrence of a phase transition between the low-and high-temperature limits, which makes it quite plausible that some features like fast scrambling [1] of black holes in the gravity dual and temperature dependence of chaotic dynamics may be retained to a certain extent at the high-temperature limit too. For instance, the numerical results obtained in Ref. [2] by exploiting the classical dynamics of the BMN model results in a fast thermalization. In Ref. [4], classical chaotic dynamics of the BFSS models is studied, and there it is found that the largest Lyapunov exponent is given as λ L = 0.2924(3)(λ t Hoof t T ) 1/4 . Therefore, the MSS bound is not obeyed only at temperatures below the critical temperature T c ≈ 0.015 , while it remains parametrically smaller than 2πT for T > T c . In Ref. [14], we have studied chaos in massive deformations of the SU (N ) Yang-Mills gauge theories in 0 + 1-dimensions, with the same matrix content as that of the bosonic part of the BFSS model, by making use of ansatz configurations involving both fuzzy 2-and 4-spheres. Our numerical results have shown very good agreement with the λ L ∝ (E/N 2 ) 1/4 -type functional behavior of the largest Lyapunov exponent with energy, which, together with the application of the virial and the equipartition theorems, allowed us to put upper bounds on the critical temperature, T c . Depending on the values of the mass parameters, our estimates for T c were around twice or about an order of magnitude larger than that obtained for the BFSS model in Ref. [4]. In the present paper, we extend and apply the methods we have developed in Ref. [14] to another interesting gauge theory, namely, the massive deformation of the Aharony-Bergman-Jafferis-Maldacena (ABJM) model [29] [30,31]. Before focusing our attention in this direction, let us also note that not only the BFSS and the BMN matrix models, but even their subsectors at small values of N are quite nontrivial many-body systems, which escape a complete solution to this day. Nevertheless, the chaotic dynamics of the smallest YM matrix model composed of two 2 × 2 Hermitian matrices with SU (2) gauge and SO(2) global symmetries has recently explored in Ref. [5] (see also Refs. [32,33] in this context) with the chaotic phase, corresponding to a toy model for a black hole, being controlled by the angular momentum associated to the rigid SO(2) symmetry. In Ref. [34], two of us explored the minimal Yang-Mills-Chern-Simons matrix model and analyzed the effect of the Chern-Simons (CS) coupling on the chaotic dynamics.
As is well known, the ABJM model is a 2 + 1-dimensional N = 6 supersymmetric SU (N ) × SU (N ) CS gauge theory at the CS level (−k, k) [29] and describes the dynamics of N coincident M 2-branes [29,35]. This model consists of four complex scalar fields C I (I : 1, 2, 3, 4), as would be expected due to the eight transverse directions to the M 2-branes, and four Majorana fermions ψ I to match the bosonic and fermionic degrees of freedom as a minimal requirement for the presence of supersymmetry. These fields are coupled bifundamentally to the SU (N ) CS gauge fields A µ andÂ µ , i.e. they carry the (N,N ) representation of the SU (N ) × SU (N ) group. The model has the R-symmetry group U (1) × SU (4) under which both the complex scalars and the fermions transform in the four-dimensional fundamental representation of SU (4) and carry +1 charge under the U (1) factor. The ABJM model is dual to type-IIA string theory on AdS 4 × S 7 /Z k (this becomes AdS 4 × CP 3 in the k → ∞ limit) via the AdS/CFT correspondence [29,35], and it possesses a massive deformation due to Hosomichi et al. [30] and Gomis et al. (GRVV) [31], preserving all the supersymmetry, but breaking the R symmetry. It is this model on which we focus our attention in the present paper. The vacuum configurations in this model are given by the GRVV matrices, which describe fuzzy 2-spheres as a somewhat intricate analysis demonstrates [31,35]. This feature is similar and comparable to the BMN model, which also has fuzzy spheres as the vacuum solutions. Our aim here is to explore the chaotic dynamics emerging from this model at the classical level as an approximation to the quantum theory in the high-temperature regime using both analytic and numeric techniques and determine upper bounds on the temperature of the system at consecutively higher matrix levels, above which the MSS bound is satisfied and below which it will eventually not be obeyed. The latter is naturally expected to occur, since, as we already noted, the classical treatment of the model could approximate the dynamics of the full quantum theory only at sufficiently high temperatures. Toward this aim, we first perform a dimensional reduction of this model from 2 + 1-to 0 + 1 dimensions by considering that the fields are spatially uniform. We work in the 't Hooft limit and focus on two distinct ansatz configurations both involving fuzzy 2-spheres, which are described in terms of the GRVV matrices. These configurations have collective time dependence, which are introduced by real functions of time multiplying the latter. Tracing over these configurations yields a family of reduced effective Lagrangians and we demonstrate that they have chaotic dynamics by computing their Lyapunov exponents. In particular, we direct our attention to examine how the largest Lyapunov exponent, λ L , changes as a function of E/N 2 . It turns out that, depending on the structure of the effective potentials, we find that either is a constant determined in terms of the CS coupling k, the mass µ, and the matrix level N . This power-law response of λ L to energy is also anticipated and supported by the exact scaling symmetry possessed by the model in the massless limit as will be discussed in the next section. Making use of our numerical results, and evoking the virial and equipartition theorems, we explore the implications for the aforementioned MSS conjecture in the context of this model. The main outcomes are the upper bounds we obtain on the temperatures above which largest Lyapunov exponents comply with the MSS bound and below which it will eventually be not obeyed. At the same time, we demonstrate that with increasing matrix level, i.e., with better numerical approximation of the 't Hooft limit, estimates of T c display a decreasing trend; put differently, the temperature range in which the MSS bound is valid expands gradually.
The paper is organized as follows. In Sec. 2, we outline and review the various features of the massive deformation of the ABJM model and obtain its reduction from 2 + 1 to 0 + 1 dimensions by assuming spatially uniform fields. In section 3, we introduce our first ansatz configuration, obtain the reduced effective actions, and present the results of the numerical analysis leading to the modeling of the energy dependence of the largest Lyapunov exponent. This is followed by the discussion explaining how we extract the temperature dependence and relate our findings to the MSS conjecture. In Sec. 4, results of an analysis following mainly the same steps of Sec. 3 are presented for another ansatz configuration. Several details of the calculations are relegated to the Appendixes A and B. We conclude in Sec. 5 by briefly summarizing our results and indicating some directions for future studies.
2 Reduction of mass-deformed ABJM mode to 0 + 1 dimensions We start with writing out the action for the bosonic part of the mass-deformed ABJM model. This is given as [30,31] where A µ andÂ µ (µ : 0, 1, 2) are two distinct gauge fields transforming under the SU (N ) k and SU (N ) −k gauge transformations, respectively. The subscripts ±k ∈ Z label the level of the Chern-Simons terms associated to these gauge fields. The potential term is given as

4)
and µ stands for the mass of the fields (Q α , R α ). e iS ABJM is invariant under the gauge group The latter is the level quantization of the Chern-Simons couplings in the action. In (2.3), we use the notation

5)
and likewise for R α 's. The ABJM model has the global SU (4) R × U (1) R R-symmetry group, which is broken down to SU (2) × SU (2) × U (1) A × U (1) B × Z 2 by the mass deformation terms given in M α and N α . If there is no mass deformation, it is suitable to formulate the theory in terms of the complex scalar fields C I , which transform under the four-dimensional fundamental representation of the SU (4) R factor and carry U (1) R charge +1. In the mass-deformed model, Q α transform under the first and R α transform under the second of the SU (2) factors of the R-symmetry group, and under U (1) A , they have the charges 1 , −1, respectively, while under U (1) B , they both have charge 1, and the Z 2 factor serves to exchange Q α and R α .
To dimensionally reduce S ABJM to 0 + 1 dimensions, we declare that all fields are independent of the spatial coordinates and depend on time only. Consequently, all partial derivatives with respect to the spatial coordinates vanish. We may introduce the notation A µ ≡ (A 0 , X i ),Â µ ≡ (Â 0 ,X i ) with (i = 1, 2). Spatial and time components of the covariant derivative are then 6) and the action takes the following form: Expressing the Chern-Simons parts of this action in terms of the covariant derivatives D 0 X i : we may as well write In (2.7) and (2.8), it is understood that all fields depend on time only. We have readily written the reduced action in the 't Hooft limit. The latter is defined as follows. While reducing from 2 + 1 to 0 + 1 dimensions, we have integrated over the two-dimensional space whose volume may be denoted, say, by V 2 . Therefore, we may introduce 1 λ t Hoof t := N V 2 and require that it remains finite in the limit V 2 → ∞ and N → ∞. In S ABJM −R , we have scaled λ t Hoof t to unity. If needed, it is possible to restore λ t Hoof t back in S ABJM −R by performing the scalings Let us note also that S ABJM −R is manifestly gauge invariant under the SU (N ) k × SU (N ) −k gauge symmetry. 2 The ground states of this reduced model are the same as that of the original model and given by configurations minimizing the potential V in (2.2). Since the latter is positive definite, its minimum is zero and is given by the configuration There are two immediate solutions to (2.9), which are given as where G α are the GRVV matrices [31,35] defining a fuzzy 2-sphere [36] at the matrix level N and c = kµ 4π . Let us note in passing that c = 0 gives a trivial solution in which both the fields Q α and R α vanish and is of no interest to us in what follows. Explicitly, G α are given as [31] 1 Note that in the original model, i.e., in 2+1 dimensions, the 't Hooft coupling is identified as λ t Hoof t = N k held fixed with N, k → ∞ [35]. In the reduced model, too, we may defineλ t Hoof t := N kV 2 = λ t Hoof t k held fixed with N, V 2 → ∞, while, in contrast, k can remain finite. In this case, scalingλ t Hoof t to 1 k is the same as scaling λ t Hoof t to unity. 2 In particular, let us note that pure CS action is indeed manifestly gauge invariant in 0 + 1 dimensions as opposed to the non-Abelian CS action in 2 + 1 dimensions, which is not. The latter gives rise to the level quantization of the CS coupling, i.e., k ∈ Z. In fact, after the reduction of the CS terms to 0 + 1 dimensions but prior to introducing the 't Hooft parameter λ t Hoof t , the effective CS coupling is simply κ := 1 4π kV 2 and is no longer an integral multiple of 1 4π due to the arbitrary volume V 2 of the two-dimensional compact space we have integrated over. This is consistent with the fact that CS term in 0 + 1 dimensions is gauge invariant and therefore its coupling is not level quantized. The latter also follows from the fact that π 1 (SU (N )) = 0 and the general considerations on the gauge symmetry properties of e iS CS which may be found, for instance, in Refs. [37,38].
with m, n = 1 , · · · , N , and they fulfill the relation We may notice at this stage that it is possible to work in the gauge with A 0 = 0 and A 0 = 0. Evaluating the variations of S ABJM −R with respect to A 0 andÂ 0 , we find the Gauss-law constraint is given by the two equations It is also useful to note that the Hamiltonian takes the form where are the conjugate momenta associated to Q α and R α , respectively. It is straightforward to see that the Hamiltonian for the CS part of the action vanishes identically as expected [38]. Let us consider the scaling transformation where ρ is an arbitrary positive constant. Under this transformation, we have Therefore, the energy scales as E → ρ −3 E. Since the Lyapunov exponent has the dimensions of inverse time, we see that it scales as in the massless limit. In the ensuing sections, we will see that this scaling of the Lyapunov exponents with energy is essentially preserved after taking the mass deformations into account.
We are now in a position to propose ansatz configurations, through which we will be able to explore the emerging chaotic dynamics. We will consider two different ansatz configurations involving the GRVV matrices and satisfying the Gauss-law constraints. Both of these ansatz configurations involve collective time dependence and are introduced via real functions of time.

Ansatz I and the effective action
The first matrix configuration we focus on is specified as where (A i ) m , (B i ) m are constants and i = 1, 2, m = 1, 2, ..., N and α = 1, 2. Thus, X i and X i are taken as diagonal matrices. No sum over the repeated index α is implied in the last line of (3.1). Here φ α (t), α(t), β(t) are real functions of time, and the Gauss-law constraint given in the Eq. (2.13) is easily seen to be satisfied by this choice of the matrices.
Evaluating the equations of motion for α(t) and β(t), we find that the emerging coupled equations have only one possible real solution and that is the trivial solution given simply as α(t) = β(t) = 0. This result is proved in Appendix (A). Henceforth, setting X i andX i to zero, inserting last line of (3.1) in the action (2.8), and performing the traces over the GRVV matrices at the level of N × N matrices, we obtain the reduced Lagrangian The corresponding Hamiltonian is where V N (φ 1 , φ 2 ) introduced in the second line denotes the potential of this reduced system and defined by the relevant expression in the first line. For k > 0, we see that is clearly positive definite, while for k < 0, this is not manifest, but it is indeed so since and t → ρ t, as can be readily expected in view of the discussion given at the end of the previous section.
To explore the dynamics of the model, we calculate the Hamiltonian equations of motion.
These take the forṁ In what follows, we will explore the dynamics emerging from the equations at µ = 1 at several different matrix levels N and the CS coupling k.
To gain some immediate information on the system, it is useful to explore its fixed points and also investigate the stability around these points at the linear level. Details of this analysis are provided in Appendix B. In brief, for k > 0, the only fixed point of this Hamiltonian system is given as (φ 1 , φ 2 , p φ 1 , p φ 2 ) ≡ (0, 0, 0, 0) with a vanishing fixed-point energy. Analysis in Appendix B shows that this fixed point is of borderline type, meaning that the linear level analysis is inconclusive to identify it as either stable or unstable character, and we do not attempt to perform a higher-order analysis. For k < 0, we find that there are several fixed points, some of which are still of borderline type. Nevertheless, the set of fixed points given as ±(∓) (for kµ < 0), as calculated in Appendix B. This may be taken as the first indication to expect the dynamics of H N to be chaotic, since the latter is usually associated to the presence of unstable fixed points in the phase space [39][40][41][42]. In fact, the systems do not tend to exhibit any appreciable chaos at energies below that of the unstable fixed points, and depending on the structure of the potential, even at energies exceeding the latter, phase space may have comparable number of quasiperiodic and chaotic trajectories; i.e., at low energies, chaos and quasiperiodic motion can coexist. In particular, a randomly picked initial condition may correspond to either a quasiperiodic or a chaotic trajectory. Therefore, it is important to pay attention to this fact in the computation of the Lyapunov exponents, and we will do so in the ensuing sections.

Chaotic dynamics and the 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 [39][40][41][42]. 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, λ are called the 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, the 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 the 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 [39][40][41][42]. The phase space for the Hamitonians H N considered in this paper are all four dimensional. From the general considerations summarized above, it is clear that the emerging chaotic dynamics of these models are governed by the largest (and only) positive Lyapunov exponent at given values of the parameters k, µ, and N .
To give a certain effectiveness to the random initial condition selection process, we adapt and use the simple approach we have developed in Ref. [14]. We briefly explain this next. Let us denote a generic set of initial conditions at First of all, we generate four random numbers and denote three of them as ω i (i = 1, 2, 3) and where E is the energy of the system. We denote the last random With the help of this relation and the energy functional, i.e., Hamiltonian given in (3.3), initial conditions are picked randomly in the form which one of the two we select is immaterial. In our calculations, we take φ 2 (0) = ω 4 ; then, φ 1 (0) is given by the solution of For some randomly picked values of ω 4 , Eq. (3.6) may not have a real solution. However, the code runs and randomly picks another ω 4 until a real solution for φ 1 (0) is obtained. Similar to the analysis performed for the Yang-Mills matrix models with massive deformations presented in Ref. [9], we set up and run a MATLAB code, which numerically solves the Hamilton equations of motion given in (3.4) at different matrix levels. We run this code 40 times for k ≥ 1 and 100 times for k ≤ −1 with randomly selected initial conditions at a given energy value E and matrix level N and calculate the average for each and every Lyapunov exponent from all runs a the final time. In the simulation, we take a time step of 0.25 and run the code from time 0 to 3000. Our code checks if the largest Lyapunov exponent has a value below a certain threshold at t = 3000 and does not include it in the averaging over the initial conditions. In our computations, we picked this threshold as 0.05 after a number of numerical trials. 3 Let us note that Lyapunov exponents below this threshold at large time (t = 3000 in our simulations) correspond essentially to the quasiperiodic trajectories in the phase space, which do not exhibit chaos but may have comparatively small or large periods and therefore usually have very small but nonvanishing Lyapunov exponents at large time, 4 and in the manner just described, we exclude them in order to obtain more precise values for the largest Lyapunov exponents of the chaotic trajectories in the phase space. In particular, we focus on H N for N = 5, 10, 15, 20, 25 at several different values of the energy.

Dependence of the largest Lyapunov exponent on energy
Since we are working in the 't Hooft limit, it is useful to consider the dependence of the largest Lyapunov exponent, λ L , on E/N 2 rather than on E, to capture the main features of the chaotic dynamics emerging from the family of Hamiltonians H N and subsequently relate it to the temperature of these systems via the use of virial and equipartition theorems.
In this case, to capture the λ L ∝ E 1/3 dependence of the Lyapunov exponent anticipated by the scaling argument given in Sec. 2, we find that it is sufficient to choose E/N 2 in the interval (0, 100). Lyapunov exponent data and the best-fitting curves of the form λ L = α N ( E N 2 ) 1/3 are given in Fig. (1) for k = 1 and N = 5, 10, 15, 20, 25. For k = 2, at sufficiently low matrix levels, the E/N 2 interval can still be taken as (0, 100), while for N > 20, it turns out to be better to stretch it to a wider range, and in Fig. 2c, we take it to be (0, 500). Lyapunov exponent data and the best-fitting curves of the form λ L = α N ( E N 2 ) 1/3 are given in the Fig. (2).
If we further increase the k value, we also need to inspect the dependence of λ L to E/N 2 in a sufficiently large range of the latter. For instance, in Fig. (3), we depict the Lyapunov data for 0 ≤ E/N 2 ≤ 500, at the matrix level N = 10 and for k = 5, 10. At higher matrix levels N and/or larger values of k, it is necessary to further increase the range of E/N 2 in order to clearly observe the E 1/3 dependence of λ L via the best-fitting curves in the form λ L = α N ( E  Tables 1, 2, and 3 given in the next subsection. In this case, we seek best-fitting curves of the form λ L = α N ( E N 2 −γ N ) 1/3 to the Lyapunov data. Motivations for considering this function of E/N 2 are twofold. For one, as noted earlier, the energies of the unstable fixed points are nonvanishing in this case and given by 108π . Since no significant chaos is present for E ≤ E F , we expect λ L 's to vanish at energies below E F (indeed, all our numerical computations show that λ L are vanishingly small for E ≤ E F ). This suggests then that γ   comparable, and for the evaluation of the coefficients α N of the fitting curves, we use the latter as they tend to produce slightly better fits.
Both for k = −1 and k = −2, we find that the (0, 500) interval for E/N 2 is sufficiently well suited to capture the energy dependence of λ L . This is corroborated by the best-fitting curves of the form λ L = α N ( E N 2 − γ N ) 1/3 given in Figs. (4) and (5). Coefficients α N for the fitting curves in Figs. (4) (5) and the respective values of γ N are provided Tables 4 and 5 in the next subsection.

Temperature dependence of the Lyapunov exponent
In Ref. [4], temperature dependence of the Largest Lyapunov exponent of the BFSS matrix model in the 't Hooft limit was determined using dimensional analysis to be of the form λ L ∝ (λ t Hoof t T ) 1/4 since λ t Hoof t and the temperature are the only dimensionful parameters of the model. Let us note that this result is consistent with the fact that for the BFSS model the potential is purely quartic and the system has a scaling symmetry implying that λ L ∝ E 1/4 and hence λ L ∝ T 1/4 temperature dependence by evoking the equipartition theorem. In Ref. [14], we focused on a mass-deformed Yang-Mills matrix theory, with the same matrix content as the bosonic part of the BFSS model, and a similar analysis is considered where the effects of the mass deformations were taken into account in a simple way. For the present matrix model, a similar approach can also be followed. We may expect that λ L ∝ (λ t Hoof t T ) α , where α is a constant that needs to be determined. As we noted in Sec. 2, we may take λ t Hoof t = N V 2 , where V 2 is the volume of the two-dimensional space we have integrated over in going from 2+1 to 0+1 dimensions. From this definition of λ t Hoof t , we see that it has the dimension of [Length] −2 , while we already know that each of λ L and T have the dimension [Length] −1 . Putting these fact together, we have the equation for α,

7)
and therefore we find α = 1/3. Thus, we may expect that λ L ∝ (λ t Hoof t T ) 1/3 . In view of the equipartition theorem, this is consistent with the λ L ∝ E 1/3 based on the scaling symmetry as discussed in the previous section. Shortly, we will see the circumstances under which the remaining dimensionful parameter, namely, the mass, may effect the relation between the energy and temperature upon the application of the virial and the equipartition theorems. To prepare for the latter, let us first obtain the total number of independent degrees of freedom (d.o.f.) of the ABJM matrix model described by the action (2.1). Here, we have both X i and X i as N × N Hermitian matrices with i = 1, 2. Each has N 2 real degrees of freedom, and therefore these give 4N 2 d.o.f. in total. We also have the fields Q α and R α with α = 1, 2 as N × N complex matrices, and each has real degrees of freedom. Subtracting these from the latter, we find the independent d.o.f. count to be 12N 2 − 2N 2 − 8 = 10N 2 − 8.
For the ansatz I, in (3.1) there is a further reduction of the independent d.o.f., which comes about as follows. Since X i andX i are null matrices due to vanishing of α(t) and β(t) as the only admissible on-shell solution as argued in Sec. 3   A, we need to subtract out a factor of 4N 2 . Additionally, because R α = 0 in this ansatz, we need to subtract another factor of 4N 2 d.o.f.. Finally, we also need to note that two equations of the Gauss-law constraint reduce to the same equation upon integrating by parts and taking the Hermitian conjugate of one or the other. Thus, the Gauss-law constraint imposes only N 2 real relations in this case. These facts bring the total number of d.o.f. count to

and demonstrated in Appendix
Let us now apply the virial theorem to the Hamiltonian in (3.3). Since the potential V N (φ 1 , φ 2 ) is not a homogeneous polynomial of its arguments, there is no exact proportionality relation linking the average kinetic and potential energies. Instead, we have

8)
where the relevant expression in the first line provides the definition ofṼ N (φ 1 , φ 2 ) introduced in the second line. The latter is positive definite for k > 0 (assuming that µ > 0, too, indeed we set µ = 1), but this is not so for negative k. In fact, the minimum ofṼ N (φ 1 , φ 2 ) is given as if both k and µ have the same sign , if k and µ have the opposite sign .

(3.9)
Applying the equipartition theorem to the kinetic energy yields where the approximation is valid at large N .
Case i: k ≥ 1: In this case,Ṽ N (φ 1 , φ 2 ) is positive definite, and therefore we have from (3.8) the inequality K ≥ V N . This and (3.10) together imply that We can express this inequality in the form

11)
where we have also dropped the brackets on energy for ease in notation.
Since we expect that λ L ∝ E 1/3 due to the scaling properties of the model and also that λ L ∝ T 1/3 as implied by the pure dimensional analysis, albeit both holding exact only in the massless limit, we are, nevertheless, led to examine the best-fitting curves of the form to profile the variation of the largest Lyapunov exponent as a function of E/N 2 . The fitting curves are plotted in Fig. (1) for k = 1, N = 5, 10, 15, 20, 25, in Fig. (2) for k = 2 at the matrix levels N = 5, 15, 25 and in Fig. (3) for k = 5, 10 at the matrix level N = 10. We observe that these fits represent the variation of the largest Lyapunov exponent with respect to E/N 2 quite well.
In fact, to evaluate the goodness of the fits in explaining the variation of λ L with respect to E/N 2 , we may inspect the square of the multiple correlation coefficient, R squared (we use R sq for short in what follows), and the residual sum of squares (SSE), which are usual statistical measures used for this purpose. The former takes a value between 0 and 1, with R squared close to 1 indicating better fits; i.e., a greater portion of the variance in the data is accounted for by the fitting curve, while the latter could take any positive value, with values close to zero indicating better fits. In the present context, R sq measures the correlation between the λ L values of the data and those predicted from the fitting curve, while SSE represents the total deviation of the predicted values from the fitting curve to the data. We find that for k ≥ 1 the fits given in Figs. (1), (2), and (3) have R sq ≥ 0.97 and with an average R sq ≈ 0.979, i.e., the fitting curves accounting for the variation of the λ L with   We are now in a position to compare and relate our result to the MSS bound λ L ≤ 2πT on the Largest Lyapunov exponent for quantum chaos [6]. This bound is conjectured to be satisfied in systems which are holographically dual to gravity, and it is shown in Ref. [7] that it is saturated for the Sachdev-Ye-Kitaev fermionic matrix model. In Ref. [4], chaotic dynamics of the BFSS matrix models are studied at the classical level, which provides an approximation to the quantum theory only in the high-temperature limit, and it was shown that the largest Lyapunov exponent disobeys the MSS bound only at sufficiently low temperatures. The authors of Ref. [4] estimated the latter to be ≈ 0.015. In Ref. [14], we study a deformation of the bosonic sector of the BFSS via two mass terms and investigating the chaos in this model via reduced effective Lagrangians, we were able to put upper bounds on the critical temperature above which MSS inequality is satisfied, and below which it will eventually be not obeyed. In Refs. [11,12], a so-called Gaussian state approximation(GSA) is introduced to investigate the quantum chaotic dynamics of the BFSS and related Yang-Mills matrix models. Results obtained in Refs. [11,12] demonstrate that the largest and all the other Lyapunov exponents tend to zero at a nonzero value of the temperature and therefore comply completely with the MSS conjecture at all temperatures. Nevertheless, it remains an open problem to show if and how the BFSS model saturates the MSS bound.
As we noted in the Introduction, the ABJM model has a gravity dual [35] via the AdS/CFT correspondence. Therefore, we may expect the MSS conjecture to hold for quantum chaotic dynamics of the ABJM model too. In this article, we are investigating the dynamics of the mass-deformed ABJM model only at the classical level, as an approximation of the quantum theory in the high-temperature limit, so we should expect that the MSS bound eventually be not obeyed at sufficiently low temperatures. 5 In other words, we expect the classical chaotic dynamics to comply with the MSS bound to a very large extent, while we also expect it to be insufficient to capture all the quantum features at low temperatures. Indeed, using (3.11) and (3.12), we find that there is a critical temperature, which we may denote as T c and is given by solving the equation which yields (3.14) From this result, we understand that for T ≥ T c the present model complies with the MSS bound on λ L , while for T ≤ T c , there is a temperature at and below which MSS bound is not respected. Thus, we may say that T c is an upper bound for the critical temperature at or below which MSS bound will eventually not be obeyed by our model. The estimated T c values at the matrix levels N = 5, 10, 15, 20, 25 are given in Tables 1 and 2 for k = 1 and k = 2, respectively, and in Table 3 at N = 10 level for k = 5, 10. We observe from the values of T c in these tables that with increasing matrix size, their values tend to decrease, which is in agreement with the fact that the 't Hooft limit is better emulated with increasing matrix size. From Table 3, we also infer that T c values tend to decrease with increasing values of k; i.e., the models with larger CS coupling tend to comply with the MSS conjecture within a wider range of the temperature.
Case ii: k ≤ −1 : In this case,Ṽ (φ 1 , φ 2 ) N is not positive definite as we have already noted; its minimum is negative and given by the expression in the second line of (3.9). Adding and subtracting M in(Ṽ 2 ) to the (3.8), we may write

15)
which implies that We may therefore write (3.17) 5 Let us note in passing that the presence of mass terms may keep the system away from saturating the MSS bound even if the full quantum dynamics could be studied. However, mass deformations lead to non-trivial vacuum solutions in the form of fuzzy sphere matrix configurations and provide us a good departure point to probe the chaotic dynamics as we do in the present paper.
Using K ≈ 3 2 N 2 T at large N , this leads to the inequality In view of this relation, we conjecture to use best-fitting curves of the form These curves are given in Fig. (4) for N = 5, 10, 15, 20, 25 at k = −1 and in Fig. (5) for N = 10, 15, 25 at k = −2. Similar to the previous case, they represent the variation of the largest Lyapunov exponent with respect to E/N 2 quite well, with R sq ≥ 0.96, with an average R sq ≈ 0.972 for the fitting curves in Figs. (4) and (5) and average SSE values ≈ 1.005.  By the same line of reasoning discussed in the previous case, using (3.18) and (3.19), we find that the critical temperature is given again as (3.14), and the numerical estimates using the α N values of the fitting curves at several different matrix levels are listed in Tables 4 and  5 for k = −1 and k = −2, respectively.
Viewing the results of the cases i and ii together, we conclude that T c values decrease with increasing N and/or |k|; i.e., the MSS bound is respected in a wider range of the temperature at matrix levels which better capture the 't Hooft limit and/or at larger values of the CS coupling.
We would like to introduce another ansatz configuration with nonzero R α and Q α matrices and examine the ensuing dynamics. We consider the ansatz while we still take X i andX i as arbitrary diagonal matrices as given in (3.1). This configuration satisfies the Gauss-law constraints given in (2.13), as can easily be checked, and the equations of motion for α(t) and β(t) yield the only real solution as the trivial solution α(t) = β(t) = 0 as shown in Appendix A. Thus, in this case, too, we set X i = 0 =X i in what follows. Substituting the matrix configuration (4.1) into the action (2.8), and performing the trace over the GRVV matrices, we obtain the effective Lagrangian as The corresponding Hamiltonian is where V N (q, r) is the effective potential defined by the relevant terms in the first two lines of (4.3). A few remarks regarding the structure of V N (q, r) are now in order. Let us first note that this potential is not positive definite for either k > 0 or k < 0, while its minimum is at zero. Next, we easily see that V N (q, r) is symmetric under the exchange of q and r. For k ↔ −k, the two terms which are proportional to 1 k change sign, but this can be compensated by exchanging q and r. Thus, we conclude that the dynamics due to this potential is independent of the sign of k.
Hamilton's equations of motion are easily obtained and are given below: To gain more insight about this Hamiltonian system, we explore its fixed points and their stability at the linear order. The details of this analysis are relegated to Appendix B. We find that for real values of µ either the set (0, ± √ kµ 2 √ 3π , 0, 0) or the set (± √ −kµ 2 √ 3π , 0, 0, 0) gives unstable fixed points for kµ > 0 and kµ < 0,respectively, while the remaining fixed points are of borderline type. The corresponding energies in either case are and the system is likely to exhibit dynamical evolution which is chaotic at and above these energies. This suggests that we may consider an offset γ (1) for the fitting curves of the form λ L = α N ( E N 2 − γ N ) 1/3 . In the next subsection, we compare this with the values of γ N implied upon the use of the virial and equipartion theorems.

Dependence λ L on energy and temperature
To obtain the profile of the mean largest Lyapunov exponent λ L with respect to the variation of E/N 2 , we numerically solve the Hamilton equations (4.4) and evaluate the mean of λ L by averaging out the largest Lyapunov exponents over 100 runs of the code with randomly selected initial conditions. For this ansatz, numerical aspects of the initial condition selection turn out to be somewhat more conveniently handled by setting q(0) = 0. Using three random numbers ω i (i = 1, 2, 3) and writing Ω i = ω i √ ω 2 i √ E as in the case of ansatz I, we generate the initial conditions in the form where the last equation in (4.6) takes the explicit form and its real roots are used to pick r(t) at t = 0, i.e., the r(0) value. Applying the virial theorem, we find that Evaluating the minimum ofṼ N (q, r), we find that it is given as (4.10) Following the same line of development and steps as in Sec. 3.3., we have From this consideration as well as the energies of the unstable fixed points, we are led to consider best-fitting curves of the form to the λ L versus E/N 2 data. Values of γ (1) N and γ (2) N are comparable, and in what follows we use the latter as they tend to work slightly better with the fitting curves.
Let us recall once again that, depending on the structure of the potential, chaos and quasiperiodic motion can coexist and may fill comparable hypervolumes of the phase space at a given energy. For the model emerging from ansatz II, we also let our code check if the largest Lyapunov exponent has a value below a certain threshold at the final time (here, we continue to use a time step of 0.25 and run the code from time 0 to 3000) and do not include it in the averaging over the initial conditions. From numerics, we found that roughly ≈ 1/5 to ≈ 1/10 of the initial conditions lead to quasiperiodic orbits at low energies, but their number too also tends to zero with increasing energy. Applying this process allows us to evaluate the average λ L value at a given energy with high precision, which is otherwise only obtained with relatively large root-mean-square errors. In our computations, we picked this threshold as 0.05 after a number of numerical trials. 6 For all cases of interest, it appears sufficient to use E/N 2 in the range (0, 100). For |k| = 1, the data points and the fitting curves are depicted in Fig. 6 for the matrix levels N = 5, 10, 15, 20, 25 and the β N coefficients of the fits are given in Table 6.
To profile the variation of λ L at larger values of CS coupling, we first inspect the case k = ±2. Data points and fitting curves are given in Fig. (7), and the corresponding β N values are listed in Table 7.
At k = ±5 , ±10 the data for λ L and the corresponding best-fitting curves are provided in Fig. 8 with β N coefficients listed in Table 8.
From the plots provided in Figs. 6, 7, and 8, we observe that the fitting curves represent the Lyapunov data almost perfectly. We find that all the fitting curves given in these figures R sq ≥ 0.99, while the SSE values generally vary around ≈ 0.002 to ≈ 0.007, except for a few cases (k = ±2, N = 25, and k = ±5 , ±10, N = 10) for which they are around ≈ 0.02. Thus the fitting curves capture the variation of the λ L with respect to E/N 2 around ≈ 99%.   To obtain critical upper bound temperatures, T c , using the coefficients β N of the fitting curves, we need to count the independent degrees of freedom of the matrices given in (4.1).  In contrast to our first ansatz (3.1), in this case, R α matrices are no longer zero. Thus, we need to note that R α (α = 1, 2) contribute 4N 2 degrees of freedom in total, while the two equations of the Gauss-law constraint (2.13) imply the same condition upon integration by parts of one or the other equation. Thus, the Gauss-law imposes only N 2 real constraints in this case, too. We therefore have n d.o.f. = 7N 2 − 8, which in the large-N limit is given as n d.o.f. ≈ 7N 2 . Therefore, Eq. (4.11) immediately leads to the inequality We find that the critical temperature is obtained by solving 14) and this yields Our estimates for the critical temperatures are given in Tables 6, 7, and 8. Let us note that classical chaotic dynamics of the family of effective Hamiltonians, H N , comply with the MSS bound for T > T c , while they will eventually not obey it at or below T c values. Similar to the result obtained for ansatz I, we notice that with increasing matrix size and/or CS coupling values k critical temperatures decrease. In particular, it is interesting to note that T c ≈ 0.0167 for N = 10 and k = ±10, which is comparably close to ≈ 0.015 found in Ref. [4] for the BFSS model, although the two models are quite different in terms of the power-law dependence of λ L s on energy (∝ E 1/3 for the ABJM model and ∝ E 1/4 for the BFSS model).

Conclusions and Outlook
In this paper, we performed a detailed study of the chaotic dynamics of the mass-deformed ABJM model. Working in the 't Hooft limit, and assuming that all the fields are spatially uniform and introducing ansatz configurations involving fuzzy spheres in the form of GRVV matrices with collective time dependence, we have obtained effective models and computed their Lyapunov exponents using numerical algorithms. Our results clearly indicate that these models possess chaotic dynamics. In particular, we directed our attention to the profile of the largest Lyapunov exponent and found that, depending on the form of the effective potential, either ) is a constant determined in terms of the Chern-Simons coupling k, the mass µ, and the matrix level N . They represent the result of the numerical findings considerably well as it is observed from Figs. (1)- (8) and also further corroborated by the λ L ∝ E 1/3 power-law dependence due to the scaling symmetry of the model in the massless limit. Upon the use of the virial and the equipartition theorems, we were able to examine the temperature dependence of the λ L 's and derived critical upper bounds, T c , on the temperature above which the MSS inequality, λ L ≤ 2πT , is respected and below which it will eventually not be obeyed. Our numerical finding for these T c values are presented in the tables given in Secs. 3 and 4, from which it is also observed that the T c values display a decreasing trend with increasing matrix size, i.e., with the better numerical emulation of the 't Hooft limit, as well as with the increasing values of the CS coupling k. We strongly feel that the next step is to devise new methods to go beyond the classical analysis presented in the present paper and explore the quantum dynamics of these models. The latter appears to be quite a formidable task. Nevertheless, inspired by the methods used in quantum chemistry in approaching many-body problems, recently a new real-time method, 7 which can be named the Gaussian state approximation, was developed and thoroughly applied to the BFSS model [11,12]. In its simplest form, GSA aims at incorporating the quantum corrections by considering a larger but a truncated set of observable whose Heisenberg equations of motion are obtained via the use of a Gaussian density matrix. Application of this method to the BFSS model demonstrated that all the Lyapunov exponents tend to zero at a nonvanishing temperature, implying that the quantum description of the BFSS model within the GSA approximation is fully compliant with the MSS inequality. However, given that it is still only an approximation of the full quantum dynamics, it falls short of providing an explicit saturation of the MSS bound by the largest Lyapunov exponent, in contrast to the result for the Sachdev-Ye-Kitaev model obtained in Ref. [7] and expected for all models with holographic duals according to the MSS conjecture. We think that it will be extremely useful to attempt to apply the GSA to the ABJM model as well to the family of effective Hamiltonians introduced in the present manuscript, not only to test the usefulness of GSA beyond the BFSS model but also to probe the quantum chaotic dynamics of the ABJM model. We hope to report on the possible developments along this direction elsewhere. and working in the A 0 = 0,Â 0 = 0 gauge, we immediately see that CS part of the action vanishes identically: Next, we evaluate Tr |D i Q 1 | 2 + Tr |D i Q 2 | 2 . We have where the indices a, b : 1 , · · · , N . Using (2.11), we obtain, for α = 1, where the last line follows since X i andX i are diagonal. The corresponding Hermitian conjugate is These give Similarly, for α = 2, with the Hermitian conjugate given as These give (A.10) We therefore have, Since S 1 S 2 − S 2 3 = 0, as one can see readily see by inspection (this is also verified using Mathematica at several different choices of N ), therefore, the only solution to the equations of motion is the trivial solution, α(t) = 0 , β(t) = 0 , (A. 15) as we intended to show.
B Fixed points of the reduced Lagrangians and their stability 2.1. Ansatz I Fixed points of a Hamiltonian system are defined as the stationary points of the phase space [14]. For the reduced dynamical system obtained using ansatz I, these points are given by the solutions of the equations (φ 1 ,φ 2 ,ṗ φ 1 ,ṗ φ 2 ) = (0, 0, 0, 0) . (B.1) Combining (B.1) and the equations of motion given in (3.4) leads to four algebraic equations, two of which are trivially solved by (p φ 1 , p φ 2 ) ≡ (0, 0). The remaining two give us the coupled algebraic equations, which are expressed as Fixed points may be determined by solving these equations. Linear stability of the system may be inspected around a given fixed point, in order to determine whether it is a stable or unstable fixed point. A similar analysis was performed in Ref. [14] and we follow it in what follows. For simplicity, let us introduce the notation (q 1 , q 2 , q 3 , q 4 ) ≡ (φ 1 , φ 2 , p φ 1 , p φ 2 ) .