Chaos in the BMN matrix model

We study classical chaotic motions in the Berenstein-Maldacena-Nastase (BMN) matrix model. For this purpose, it is convenient to focus upon a reduced system composed of two-coupled anharmonic oscillators by supposing an ansatz. We examine three ansätze: 1) two pulsating fuzzy spheres, 2) a single Coulomb-type potential, and 3) integrable fuzzy spheres. For the first two cases, we show the existence of chaos by computing Poincaré sections and a Lyapunov spectrum. The third case leads to an integrable system. As a result, the BMN matrix model is not integrable in the sense of Liouville, though there may be some integrable subsectors.


Introduction
One encounters non-linear differential equations everywhere in theoretical physics. String theory is no exception, as a matter of course. For example, when a string is moving on a curved space, the classical equation of motion becomes non-linear generally. If the curved background has a special coset-structure like AdS 5 ×S 5 , the string theory is classically integrable in the sense that the Lax pair exists [1]. Then classical solutions can be studied analytically, for example, by using the classical inverse scattering method. However, in most cases, one cannot solve analytically non-linear differential equations. The systems are non-integrable and exhibit chaos.
Recently, chaotic string solutions have been well studied in various ways [2][3][4][5][6][7][8][9][10][11][12]. Motivation lies on potential applications to the AdS/CFT correspondence [13][14][15] beyond integrability [16], but it has not been clarified yet what corresponds to chaotic strings in the gauge-theory side. Note that chaotic motions cannot be described with analytic functions. This is a characteristic of chaos. So it may be intriguing to ask whether the correspondence holds even in the non-analytic region or not. One can test it at a deeper level.
Along this direction, one may focus upon the M-theory dynamics. A matrix model description of light-front M-theory is proposed by Banks-Fischler-Shenker-Susskind [17]. It is called the BFSS matrix model. Chaos in this model was studied in [18] by following the way in classical Yang-Mills (YM) theories [19,20]. However, there are some unsatisfactory points. In the BFSS matrix model, there are flat directions which may cause JHEP06(2015)191 unbounded motions. Then the motions may conflict with the boundedness of chaos during long-time numerical computation of classical trajectories. Another point is a technical problem. In [18], the largest Lyapunov exponents were computed by reducing the system to simple models. However, it does not seem that the results converge on definite values. In addition, the time scale of (would-be) convergence is also much larger than the inverse of the Lyapunov exponent. Thus, it would be worth revisiting chaotic motions in the BFSS matrix model. A more appropriate system is the Berenstein-Maldacena-Nastase (BMN) matrix model [21]. It contains mass terms 1 and flat directions are lifted up. Hence there is no subtlety with the boundedness of chaos. In addition, the Myers term may give rise to a separatrix structure in the phase space. It means a bifurcation of the behavior of classical solutions and may be a source of chaos. In total, the BMN matrix model is a good laboratory to study chaotic motions.
In this paper, we will study classical chaos in the BMN matrix model. For this purpose, it is useful to reduce the full system to a simple model composed of two-coupled anharmonic oscillators by supposing an ansatz. We examine three cases: 1) two pulsating fuzzy spheres, 2) a single Coulomb-type potential, and 3) integrable fuzzy spheres. For the first two cases, we show the existence of chaos by computing Poincaré sections and a Lyapunov spectrum. The third case leads to an integrable system. As a result, the BMN matrix model is not integrable in the sense of Liouville, though there may be some integrable subsectors.
This paper is organized as follows. In section 2 we give a brief review of chaos. In particular, the characteristics of chaos in the Hamiltonian system. Section 3 considers the Hénon-Heiles system as an example of non-integrable Hamiltonian systems. In section 4 we consider chaotic motions and an integrable subsector in the BMN matrix model. Section 5 is devoted to conclusion and discussion. In appendix A, the computation scheme of Lyapunov spectrum is summarized.

What is chaos?
Probably, anybody has said the term, "chaos". For example, when you see the top of the desk in a mess, you would say "It's chaos". So it seems likely that "chaos" implies that some things are completely messed up. In fact, this point captures an aspect of chaos, randomness.
There would be no universally accepted definition of chaos. The following is a list of characteristics of chaos commonly seen in almost any books: 1. High sensitivity to initial conditions, 2. Randomness of motions (continuous spectrum), 3. Non-periodic and bounded motions (non-analyticity).
The item 1) is significant and it is sometimes called the Butterfly effect. One may know the following famous sentence, "Does the flap of a butterfly's wings in Brazil set off JHEP06(2015)191 a tornado in Texas?" Here the flapping wing corresponds to a small change of the initial condition and it may cause a chain of events leading to large scale phenomena.
Suppose a continuous dynamical system described by a set of non-linear differential equations. The system is deterministic because the dynamics is determined by fixing an initial condition. But one cannot predict the motions generally after long time has passed because of the non-linearity. That is, it is impossible to forecast weather in the far future. This indicates that some information is lost dynamically and it is measured by the Kolmogorov-Sinai (KS) entropy. This entropy production is concerned with the randomness [the item 2)].
The item 3) can be intuitively understood. If the motion is periodic, then it does not satisfy the other properties. It exhibits the line spectrum and no sharp sensitivity to initial conditions. To see that the boundedness is necessary, let us consider a simple systemẋ = x . The orbits run away exponentially to infinity, hence it seems to exhibit high sensitivity to initial conditions. However, it just means an attracting fixed-point at infinity and does not indicate chaotic motions. The recurrence of motions should be required implicitly so as to exclude fixed points. Non-analyticity follows from the continuous spectrum. The smooth behavior in the Fourier space indicates a random, non-analytic motion in the real space (i.e., the motion is continuous but not differentiable everywhere). Hence chaotic motions do not behave in an analytic way in comparison to (quasi-)periodic ones. Indeed, if the motion is described by an analytic function, then one can predict the motion for an arbitrarily long time and it does not satisfy the item 1).
There are three methods to see chaos, 1) Poincaré sections, 2) Lyapunov spectrum and 3) power spectrum. Since some readers are not familiar to these issues, it would be helpful to introduce each of them below.

Poincaré sections
To study complex trajectories in the phase space, it is helpful to use a Poincaré section, which is a surface of section intersecting the flow in the phase space. It is defined so that all of the possible orbits cross the section repeatedly. The dimension of the flow (i.e., the dynamical system) is the number of the first-order ordinary differential equations, 2 and the dimension of Poincaré sections is one lower than it. On a Poincaré section, an orbit is mapped to dots, which are the intersection points of the orbit and the Poincaré section. Since the dynamics is governed by a set of differential equations, the original trajectory can be reconstructed from the dots on the Poincaré section. In this sense, the flow and the Poincaré section are equivalent, while there is a degree of freedom to choose the surface in practice. It is the usual to pick up only the dots at which the trajectory intersects from only the one side to the other side. This constraint allows any dots on a Poincaré section to uniquely determine the next point of the trajectory in autonomous Hamiltonian systems. Suppose a Hamiltonian system with a four-dimensional phase space JHEP06(2015)191 (q 1 , q 2 ; p 1 , p 2 ) . A Poincaré section is a two-dimensional surface, for example, with q 2 = 0 and p 2 > 0 .
By changing initial values, a set of trajectories is generated and it contains several kinds of motions. The Poincaré section is helpful to figure out their behaviors. Periodic orbits are projected onto the Poincaré section as finite numbers of dots on closed contours. A periodic trajectory comes back to a point on the Poincaré section after its period. Quasi-periodic orbits are also projected as dots on closed contours, but the resulting portrait becomes dense. They are understood as non-resonant tori. Completely chaotic orbits result in randomly scattered dots on the Poincaré section. These behaviors depend on initial values and generally chaotic motions coexist with (quasi-)periodic motions on the section.
Suppose that a small perturbation makes an integrable system be non-integrable. Then the Kolmogorov-Arnold-Moser (KAM) theory [23][24][25] claims that trajectories incline to be chaotic gradually as the perturbation becomes larger. When the perturbation is small enough, most of the orbits are (quasi-)periodic. These are called KAM tori, which are the remnants of Liouville tori in the unperturbed integrable system. In other words, almost all of the Liouville tori survive the small perturbation, though the tori are distorted and translated. As the perturbation becomes larger, chaotic trajectories appear and coexist with KAM tori. Then the region of chaotic trajectories enlarges while the KAM tori tend to get broken. When the perturbation is larger, KAM islands and islets are in the sea of chaos.

Lyapunov spectrum
The next is to introduce the Lyapunov spectrum. It is very useful to measure the dynamical information loss.
Let us first introduce a Lyapunov exponent. We basically follow the description in [26]. Suppose that a solution x(t) is chaotic and consider a nearby point x(t) + δ(t) . Then the deviation δ(t) grows exponentially like (2.1) Here λ is a positive constant called a Lyapunov exponent. This behavior indicates that the deviation of the initial value δ(0) is amplified exponentially even if it is very tiny. This is nothing but a realization of sensitive dependence on initial conditions. In general, the phase space is higher-dimensional. Then there are N different Lyapunov exponents for an N -dimensional phase space. An intuitive explanation is the following. Consider time evolution of an infinitesimal sphere of perturbed initial configurations. During its evolution, the sphere tends to be distorted into an infinitesimal ellipsoid. Let δ k (t) (k = 1, . . . , N ) denote the length of the k-th principal axis of the ellipsoid. Then the deviations behave as λ k : Lyapunov exponents.
The set of the exponents λ k generates a spectrum called the Lyapunov spectrum. For large t , the diameter of the ellipsoid is controlled by the most positive λ k , and it is called the largest Lyapunov exponent.

JHEP06(2015)191
The Lyapunov spectrum is also related to the dynamical entropy production. For example, the sum of all of the positive Lyapunov exponents is identified with the KS entropy when the system is ergodic. It would be very interesting to look for an application of the KS entropy in the context of string theory and M-theory.

Power spectrum
The power spectrum is defined as where y(t)y(s) is a two time correlation function. It is quite helpful when one needs to distinguish quasi-periodic motions and chaotic motions. When the ratio of periods is irrational, the former leads to dense orbits in the phase space (the denseness). On the other hand, the latter is related to a stochastic motion (the randomness). In the Hamiltonian system, there is no attractor because the phase space volume has to be preserved due to Liouville's theorem. Hence it is often difficult to distinguish them at a first glance.
However, the power spectrum shows clear distinction. For quasi-periodic motions, the spectrum contains sharp peaks associated with periods (the line spectrum) . On the other hand, the continuous spectrum appears for chaotic motions.
In the following, we will focus upon the first two methods, Poincaré sections and a Lyapunov spectrum.

Warming up -the Hénon-Heiles system
In this section, let us consider the Hénon-Heiles system [27]. This is a famous example of non-integrable Hamiltonian systems containing chaotic motions. This model was originally proposed by Hénon and Heiles (1964) in the study of celestial mechanics.

Hénon-Heiles system. The Hamiltonian is given by
where a is an arbitrary real constant parameter. The system is composed of two particles with the same mass (normalized to 1) and the interaction is third-order. The positions of the particles are described by q 1 and q 2 , and the conjugate momenta are p 1 and p 2 . It is well-known that this system is non-integrable. Before going to the detailed analysis, it may be interesting to see that the form of the interaction term in (3.1) is very important. If the interaction term is slightly modified like then the system becomes classically integrable. The resulting system is often called "anti Hénon-Heiles system".

JHEP06(2015)191
The system with (3.1) is equivalent to a periodic lattice with three particles, The unperturbed Hamiltonian H 0 describes a system of three connected harmonic oscillators. Hence the equations of motion for H 0 are exactly solved by using the normal coordinates. By using the normal coordinate for H 0 , one can rewrite the total Hamiltonian H . Then, after dropping the total momentum, the Hénon-Heiles system is reproduced.
Here one may sum up higher-order interaction terms with appropriate numerical coefficients and obtain an exponential-type potential. The resulting system is nothing but a Toda lattice system [28]. Hence, inversely speaking, the Hénon-Heiles system may be regarded as a truncation of the Toda lattice system.
Poincaré sections and Lyapunov spectrum. Let us see chaotic motions in the Hénon-Heiles system (3.1) . It is a good exercise to understand chaotic motions in the Hamiltonian system. Here we will consider Poincaré sections and a Lyapunov spectrum. These quantities can be evaluated numerically and the results are shown in figure 1. Figures 1 (a), (b) and (c) are Poincaré sections for E = 0.0833, 0.1250 and 0.16667 , respectively. The sections are taken with q 1 = 0 andq 1 > 0 . We set that a = 1 . In figure 1  (a), most of the solutions are (quasi-)periodic (KAM tori). In figure 1 (b), chaotic motions appear but still (quasi-)periodic motions coexist. In figure 1 (c), most of the solutions are chaotic. Figure 1 (d) is a Lyapunov spectrum with E = 0.16667 . It is computed by taking an initial condition randomly from the region around q 1 = 0.0, p 1 = 0.35, q 2 = −0.115, p 2 = −0.443 with the radius 5.0 × 10 −3 . The largest Lyapunov exponent is λ 1 ∼ 0.15 . This positive value indicates high sensitivity to initial conditions and provides a positive support for the existence of chaos in this system. As a feature of the Hamiltonian system, the spectrum is symmetric and the total sum of all of the exponents becomes zero. The remaining two exponents become zero due to the energy conservation.
The features presented in figure 1 are characteristics of chaos in the Hamiltonian system. In the next section, we will show similar results in the case of the BMN matrix model.

Chaotic motions in the BMN matrix model
We study here classical chaotic motions in the BMN matrix model. It is useful to reduce the full theory to a simple model so as to extract a mechanism which induces chaotic motions. We examine three reduction ansätze and argue whether chaotic motions exist or not.

Setup
The BMN matrix model [21] was originally proposed as a realization of M-theory on the maximally supersymmetric pp-wave background [29]. That is, it can be regarded as a deformation of the BFSS matrix model [17]. This model can also be derived as the matrix  regularization of a supermembrane theory on the pp-wave [30,31] by following a seminal work [32]. The BMN matrix model is a one-dimensional U(N ) gauge theory composed of matrixvalued variables X r (t) (r = 1, . . . , 9) , a gauge field A and 16 fermions. For the present purpose, we will focus upon the bosonic sector.
The bosonic part of the classical action is given by where i, j = 1, 2, 3 and a, b = 4, . . . , 9 . The covariant derivative D t is defined as The symbol "·" denotes the derivative with respect to t . The real constant parameter µ

JHEP06(2015)191
measures the deformation. When µ = 0 , the system is reduced to the BFSS matrix model. When µ = 0 , mass terms and the Myers term [33] are turned on. 3 In the BMN matrix model, the classical trajectories are bounded due to the mass terms. Hence there is no flat direction in comparison to the BFSS case. This characteristic would be so important as to study chaotic behavior definitely for sufficiently long time.
Due to the presence of the Myers term, there exist fuzzy sphere vacua as well as the trivial vacuum X r = 0 . Thus the BMN matrix model has the rich vacuum structure. Furthermore, the Myers term leads to a separatrix structure in the phase space and may provide a new source of chaotic motions.
We will work with the Weyl gauge A = 0 hereafter. Then the classical equations of motion are given by In addition, one has to take account of the Gauss law constraint, This comes from the equation of motion of A .
In the next subsection, we will study classical solutions of the equations of motion (4.2) equipped with the Gauss law constraint (4.3) .

Reduction to simple models
To argue classical chaotic motions in the BMN matrix model, it is sufficient to study a simple reduced system by imposing an ansatz. We will examine three examples of the reduction ansatz below.

JHEP06(2015)191
With the ansatz (4.4) , the classical equations of motion (4.2) are reduced to 0 =r + µ 3 2 r − µr 2 + 2r 3 + 2rx 2 , 0 =ẍ + µ 6 2 x + 2x 3 + 2xr 2 . (4.5) The resulting system (4.5) is a non-linear dynamical system in which two anharmonic oscillators are coupled. Note that the system (4.5) can also be obtained from the Lagrangian, . This is the onset of reordered motions. This phenomenon is quite natural because motions tend to be harmonic oscillation as the energy becomes much higher.
The Lyapunov spectrum for E = 10 and µ = 2.0 is presented in figure 2 (f). It is computed by taking an initial condition randomly from the region around r = 0.0, p r = 3.8, x = −0.50, p 2 = 2.33 with the radius 5.0 × 10 −3 . The largest Lyapunov exponent is approximately 0.302 , which is significantly non-zero. This result gives a support for the existence of chaos.
(2) A single Coulomb potential (N = 3). As the next example, let us consider the following 3 × 3 ansatz:

JHEP06(2015)191
Here x(t) and y(t) are real functions to be determined. The angle variable ϕ(t) is cyclic and the associated constant of motion is L ≡ y 2φ . This ansatz (4.7) was originally proposed by Arnlind and Hoppe [34]. Then the resulting differential equations are given by The system (4.8) can also be obtained from the Lagrangian, This system can be seen as a system composed of two-coupled anharmonic oscillators with a repulsive Coulomb-type potential. The Lyapunov spectrum in figure 3 (d) also supports the existence of chaos. It is computed by taking an initial condition randomly from the region around x = 0.95, p x = 0.0, y = 0.5, p 2 = 1.44 with the radius 5.0 × 10 −3 .
(3) Integrable fuzzy spheres (N = 4). Let us concentrate on the SO(3) subsector and consider the following ansatz: Here σ i are the standard Pauli matrices and hence the matrix size of X i is 4 × 4 . Then r(t) and x(t) are real functions to be determined. The ansatz (4.10) satisfies the Gauss law constraint (4.3). With the ansatz (4.10) , the equations of motion (4.2) are boiled down to x + −2µr + 6r 2 x .  One can derive this system (4.11) from the following Lagrangian: (4.12) Note that the system (4.11) is completely integrable. This means that there exists the second integral of motion as well as the Hamiltonian. In fact, the integrability of this system was shown in [35] . In order to show the integrability, it is helpful to shift r as r = y + µ/6 . Then the resulting Lagrangian is where the constant term has been dropped off. Now it is easy to compare the expression (4.13) with the results in [35] , and the second integral is given by (4.14)

JHEP06(2015)191
It may be interesting to rewrite the Hamiltonian H by subtracting I . One can easily show the following expression: Here R ≡ x+y and P is the conjugate momentum for R . Thus the Hamiltonian H = H R +I is nothing but a system of two independent oscillators. Finally, it should be remarked that the following orthogonal transformation makes the integrability apparent: It recasts X i into the block-diagonal form that indicates that the system is composed of two independent pulsating fuzzy spheres with the radial coordinates r + x and r − x , respectively. Although this ansatz seems rather trivial, it would be worth to mention that the BMN matrix model contains an integrable subsystem.

Conclusion and discussion
In this paper, we have studied classical chaotic motions in the BMN matrix model. We have examined three ansätze: 1) two pulsating fuzzy spheres, 2) a single Coulomb-type potential, and 3) two coupled fuzzy spheres. For the first two cases, we have shown the existence of chaos by computing Poincaré sections and a Lyapunov spectrum. The third case has led to an integrable system. As a result, the BMN matrix model is not integrable in the sense of Liouville, though there may be some integrable subsectors. It would be worth to comment on preceding works on the integrability of the BMN matrix model [36][37][38][39][40]. The integrability argument is based on the fact that the BMN matrix model is obtained as a consistent Kaluza-Klein (KK) reduction of the N =4 super Yang-Mills (SYM) theory on R × S 3 [36]. The dilatation operator of the N =4 SYM on R 1,3 is identified with the Hamiltonian of the N =4 SYM on R × S 3 via the operatorstate mapping. It is now realized that the dilatation operator is closely connected with an integrable spin chain Hamiltonian [41][42][43]. Hence the integrability of the dilatation operator might be inherited to the BMN matrix model through the KK-reduction. Note here that we have considered classical chaos in the BMN matrix model, which is concerned with the classical Hamiltonian of the N =4 SYM on R × S 3 (or the classical dilatation operator D cl of the N =4 SYM on R 1,3 ) . On the other hand, the spin-chain description is realized at the quantum level and the quantum dilatation operator D qu plays a central role. 4 4 More precisely, Dqu is divided as follows: The first term Dtree is concerned with the free Hamiltonian of the BMN matrix model, in which there is no chaotic motion as a matter of course.

JHEP06(2015)191
Thus our analysis has indirect relevance with the preceding one. It may be interesting to study whether the classical chaos can survive the quantization. The BMN matrix model was originally proposed in the context of M-theory and it can be regarded as a deformation of the M(atrix) theory [17]. It would be very intriguing to try to understand the chaotic motions from the point of view of the D0-brane dynamics. With some generalizations, those may capture an aspect of black hole physics such as the fast scrambling process [44]. For interesting works along this direction, see [45,46]. In the very recent, a bound on a Lyapunov exponent has been argued [47].
It seems promising that non-linear dynamical perspectives of string theory and Mtheory would be more and more important. The way and result presented here would be a key ingredient in tackling the fundamental problems like information loss due to a black hole formation.

Acknowledgments
It is a pleasure to acknowledge helpful discussions with Sinya Aoki, Masanori Hanada, Takeshi Matsumoto and Shin-ichi Sasa. The work of YA was supported by the Japan Society for the Promotion of Science (JSPS). This work is supported in part by the JSPS Japan-Hungary Research Cooperative Program. A part of the numerical calculations was carried out on SR16000 at YITP in Kyoto University.

A How to compute the Lyapunov spectrum
In section 2, we have explained the notion of the Lyapunov spectrum in an intuitive way. More practically, in order to numerically evaluate the Lyapunov spectrum, we have to follow a mathematical formulation. The basic method is based on the expansion rate of the volume [48,49]. With this method, all of the N Lyapunov exponents can be computed.
Consider a system of the first differential equations described bẏ Then, by taking a variation, the following relation is obtained, Hence a deviation vector δx (t) can be expressed as Here U t is a time evolution operator and δx (0) is the initial value of the deviation vector. This is a solution of (A.2). The linearity of U t leads to the following relation:  Then the Lyapunov exponent for this deviation vector is defined as Here in the last expression T is divided into n small steps like T = n∆t . The convergence of this limit is proven by Oseledec [50]. Let us next introduce a fancy trick to compute the Lyapunov spectrum composed of N exponents. The first is to prepare an orthonormal basis set {v 1 0 , v 2 0 , . . . , v N 0 } to span the tangent vector space along the trajectory. Then we make each basis v i 0 evolve by U ∆t . The resulting set of the vectors {u 1 1 , u 2 1 , . . . , u N 1 } cease to be one of the orthonormal vectors, where u i 1 = U ∆t v i 0 . This is depicted in figure 4. Then it is necessary to rebuild a new set of the orthogonal vectors { v 1 1 , v 2 1 , . . . , v N 1 } from {u 1 1 , u 2 1 , . . . , u N 1 } by following the Gram-Schmidt process and compute the expansion rate a k for each vector v k 1 : Lastly, the orthogonal vector set { v 1 1 , v 2 1 , . . . , v N 1 } is normalized as Thus a single cycle has been completed. Then the same cycle should be repeated n times. This means that the set {v 1 i , v 2 i , . . . , v N i } is used as an initial orthonormal vector set for the (i + 1)-th cycle. Now the largest Lyapunov exponent λ 1 is defined as

JHEP06(2015)191
By construction, v 1 i heads toward the direction which has the highest sensitivity to initial conditions at sufficiently late time. Thus the expansion rate a 1 i has the largest value in this region. Then it is obvious from (A.5) that the largest Lyapunov exponent is given by the sum of logarithm of the expansion rate.
In a similar way, the k-th Lyapunov exponent λ k is defined as Thus the set of the exponents {λ 1 , . . . , λ N } is obtained. This set is called the Lyapunov spectrum. Note that the computation time can be estimated by the inverse of the smallest positive Lyapunov exponent. This is because all of the exponents converge within this scale.
A characteristic of the Lyapunov spectrum in the Hamiltonian system is that the total sum of the exponents becomes zero, i.e., Liouville's theorem says that the phase-space volume is conserved in the Hamiltonian system (i.e., V (t) = V (0)). Thus the relation (A.10) must be satisfied. As a side note, for the dissipative system, the following inequality is satisfied: This is just because the volume V (t) monotonically decreases.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.