Chaos from equivariant fields on fuzzy S4

We examine the 5d Yang-Mills matrix model in 0 + 1-dimensions with U(4N) gauge symmetry and a mass deformation term. We determine the explicit SU(4) ≈ SO(6) equivariant parametrizations of the gauge field and the fluctuations about the classical four concentric fuzzy four sphere configuration and obtain the low energy reduced actions(LEAs) by tracing over the SF4s for the first five lowest matrix levels. The LEAs so obtained have potentials bounded from below indicating that the equivariant fluctuations about the SF4 do not lead to any instabilities. These reduced systems exhibit chaos, which we reveal by computing their Lyapunov exponents. Using our numerical results, we explore various aspects of chaotic dynamics emerging from the LEAs. In particular, we model how the largest Lyapunov exponents change as a function of the energy. We also show that, in the Euclidean signature, the LEAs support the usual kink type soliton solutions, i.e. instantons in 1+ 0-dimensions, which may be seen as the imprints of the topological fluxes penetrating the concentric SF4s due to the equivariance conditions, and preventing them to shrink to zero radius. Relaxing the Gauss law constraint in the LEAs in the manner recently discussed by Maldacena and Milekhin leads to Goldstone bosons.


JHEP12(2018)015
Matrix models associated to M-theory and string theories [1][2][3] have been under investigation from various perspectives ever since their discovery over twenty years ago. The broad span of interest on the subject is reflected in the literature (for a recent review, see [4] and references therein.) Among these, the BFSS model [1] is a supersymmetric U(N ) gauge theory consisting of nine N × N matrices in its bosonic part, whose entries depend on time only and it also goes by the common name of matrix quantum mechanics in the literature. It is associated to type II-A string theory [4][5][6] and appears as the DLCQ (discrete lightcone quantization) of M-theory on flat backgrounds [4][5][6]. The massive deformation of the BFSS theory, preserving the supersymmetry, is known as the BMN model and describes the DLCQ of M -theory on pp-wave backgrounds [2,3]. These matrix models describe systems of N coincident D0-branes, respectively in flat and spherical backgrounds. The latter is due to the fact that fuzzy 2-spheres appear as vacuum configurations in the BMN model. At large N and strong coupling/low temperature limit, D0-branes form a black brane, i.e. a string theoretic black hole [4,5]; the structure of this gravity dual is discussed in varying amount of detail in several references [4,5]. In [7], a model on how the fuzzy spheres in the BMN model collapse to form a black hole is discussed. The perspectives gained from the matrix models have recently started to motivate numerous investigations oriented to acquire new information on the properties of black holes, such as their thermalization, evaporation processes as well as their microscopic constituents via the study of BFSS and BMN models at large N , using both analytical insights and ever increasingly numerical techniques [8][9][10][11][12][13][14][15][16][17]. In some of these works, for the BMN and BFSS models in the large N and high temperature limit numerical evidence for the fast thermalization is obtained [11][12][13][14]. From a more general perspective, the latter is an example of the fast scrambling conjecture, which has been proposed by Sekino and Susskind [18] and which may be stated as the fact that in black holes, chaotic dynamics set in faster than in any other physical system and the rate at which this occurs is logarithmic in the number degrees of freedom of the black hole, that is, it is proportional to the logarithm of its Bekenstein-Hawking entropy. Even though the aforementioned limit is distinct from the one in which the gravity dual is obtained, it is the natural limit in which the classical mechanics provides a good approximation to the quantum mechanical matrix model. 1 This limit is free from fermions as the latter contributes to the dynamics of the bosonic matrices only at low temperatures. It has been also noted that [14], since, numerical studies performed so far do not show a phase transition occurring between the low and the high temperature limits of these matrix models [8-10, 19, 20], it is reasonable to expect that features like fast scrambling of blackholes in the gravity dual could survive at the high temperature limit too. In fact, chaotic dynamics in the BFSS model is studied in [14] in this classical limit by calculating the Lyapunov spectrum, where it was also demonstrated that a classical analogue of fast scrambling is valid for this system as the scrambling time JHEP12(2018)015 is found to be proportional to log N 2 . In [21] simulations of the BFSS model is performed at intermediate temperatures to numerically study the black hole horizon.
Even the matrix models at small values of N appear to be highly non-trivial manybody system whose complete solution evades us to date. Recently, there also has been some interest in examining the chaotic dynamics emerging from such models [15,22] (See, also [23], for earlier attempts). These studies also provided some qualitative implications on the black hole phases in such models; for instance, in [15] it has been argued that the edge of the chaotic region found in a SU(2) YM matrix model with only two matrices, i.e. the smallest matrix model with non-trivial dynamics, corresponds to the end of the black hole phase. In [24] a novel approach has been developed to estimate the ground state energy of this smallest matrix model. Authors of [22], considered simple ansatzes for the BMN model at N = 2, 3 satisfying the Gauss law constraint to probe the chaotic dynamics.
Fuzzy two sphere and its direct sums are not the only compact spherical geometries appearing in M-theory. In fact, it has been known for a quite long time that fuzzy four spheres make their appearance in matrix models as longitudinal five branes [25]. For the purposes of this paper however, presence of fuzzy four sphere solutions in Yang Mills 5matrix model with massive deformation term plays the central role. Pure YM 5-matrix model has been known in the literature for quite a while [26] and may be obtained as the reduction of the YM theory in 5 + 1-dimensions to 0 + 1-dimensions keeping only the time dependence of the matrix elements. Together with the mass term it can also be conceived as a deformation of a subsector of the bosonic part of the BFSS model, as we will explain in more detail in the next section. Contrary to the fuzzy two sphere solutions in the BMN model, fuzzy four spheres in the mass deformed YM-matrix model are classical solutions for negative mass squared (µ 2 = −8), which may be an indication of tachyonic instabilities. Nevertheless, it was recently shown by Steinacker [27] that in pure YM 5-matrix model, one-loop quantum corrections stabilizes the radius of the fuzzy four sphere and prevents its collapse by shrinking to zero radius. In this paper, we mainly focus on a mass deformed U(4N ) YM 5-matrix model and consider the exact parametrization of SU(4) equivariant fluctuation modes about the four concentric fuzzy four sphere configurations. Using the equivariant parameterization of the gauge field and the fluctuations, we perform the traces over the fuzzy four spheres at first five lowest lying levels and obtain the corresponding low energy reduced actions (LEAs). We demonstrate that the potentials of all of these reduced effective actions are bounded from below, from which we infer that the negativity of µ 2 does not actually cause any instabilities under equivariant fluctuations. As we will briefly discuss in section 3, this feature of our treatment may also be viewed as a consequence of the fact that the equivariant parametrization of the fluctuations introduces topological fluxes through the fuzzy four sphere, preventing it to shrink to zero radius.
Equivariant parametrizations breaks the U(4) symmetry of the concentric S 4 F configuration down to U(1) × U(1) × U(1) and this is further reduced to only U(1) × U(1) in LEAs as one of the gauge fields completely decouple after tracing over S 4 F . The gauge fields in the reduced actions are not dynamical, and their equations of motion lead to the constraints, which may, in fact, be seen as the residue of the Gauss law constraint on the matrix model enforcing the physical states to be the singlets of the gauge symmetry group.

JHEP12(2018)015
In the LEAs the latter condition simply translates to the requirement that the two complex fields appearing in the LEAs be real, that is, uncharged under the abelian gauge fields. This breaks the U(1) × U(1) symmetry further down to Z 2 × Z 2 .
We utilize the LEAs to explore the chaotic structure emerging from the matrix model with the fuzzy four sphere background. For the reduced action obtained after tracing at the matrix levels 4N = 16, 40, 80, 120 and 224 corresponding to the first five S 4 F levels, n = 1, · · · , 5, we numerically solve the Hamilton's equations of motion and compute the Lyapunov spectrum at several different energies, revealing the chaotic dynamics. We explore various features of the chaotic dynamics using our data. We show that the Largest Lyapunov Exponents (LLE) have a dependence on energy, which fits very well with the functional relation λ n (E) = α n + β n 1 √ E . The data on LLEs also enables us to probe the onset of chaos in the LEAs' dynamics. In fact, we are able to estimate the energies at which appreciable amount of chaotic dynamics is observed at each matrix level (n = 1, · · · , 5) and also compute the rate at which LLEs change to be proportional to E − 3 2 . Except at n = 1 the phase space of the LEAs are all ten dimensional, meaning that there are ten Lyapunov exponents associated to the LEAs at the matrix levels n ≥ 2. At n = 1, however, out of five of the generalized coordinates and corresponding velocities in the LEA, three of them combine to appear only in a particular form thereby reducing the dimension of the phase space to six. At low energies n = 1 model exhibit different features compared to those for n ≥ 2 and this is discussed in through detail section 4.2. Plots of the time development of all the Lyapunov exponents at several different values of energy at all of the matrix levels (n = 1, · · · , 5) is given in the appendix B.4 and exhibit, in particular, that all the Lyapunov exponents at a given energy sum up to zero, as is expected to happen in all Hamiltonian systems. Some technical features of the computation of the Lyapunov spectrum is outlined in section 4.2 for completeness.
The paper is organized as follows. In section 2, for completeness, we give a brief review of S 4 F and how it appears in Yang-Mills matrix models. In section 3, we first determine the exact parametrization of the gauge field and the fluctuations, which are restricted to transform as a scalar and vectors, respectively, under the combined adjoint action of SO(5) isometry of S 4 F and the SU(4) gauge symmetry generators in SO (5). In this section we also obtain the LEAs by tracing over the S 4 F at several different matrix levels, and elaborate on their structure. In section 4 we focus on the dynamical structure of the LEAs. After discussing the implications of the Gauss law constraint, we present our results exhibiting the chaotic dynamics emerging from the LEAs. In section 5, we examine the properties of the LEAs in Euclidean signature, and make evident through a number of examples that, they possess Z 2 kink solutions, i.e. instantons in 1 + 0-dimensions. These may be seen as the imprints of the non-trivial topological fluxes piercing the S 4 F , which were mentioned above. Motivated by the recent work of Maldacena and Milehkin [28] on relaxing the Gauss law constraint in BFSS and BMN models (see [29] for supporting numerical work), in section 6, we revisit the gauge symmetry of the LEAs and present a concise treatment on the consequences of not imposing the Gauss law constraint, which leads us to conclude the presence of massless excitations (Goldstone bosons) associated to these LEAs. We close the paper by summarizing our findings. Appendices collect reference formulas, intermediate

JHEP12(2018)015
steps of some of the analytic calculations, explicit form of the LEAs for n ≥ 2, the whole sets of the corresponding minima of the associated potentials, as well as all the time series plots of the Lyapunov spectrum at all matrix levels (n = 1, · · · , 5) and at several different values of the energy.
2 Fuzzy S 4 in Yang Mills matrix models

Basics
We launch the developments in this section by considering the Yang-Mills 5-matrix model in Minkowski signature and with U(4N ) 2 gauge symmetry, whose action may be given as [4,26,27] where X a (a : 1, . . . 5) are 4N × 4N Hermitian matrices transforming under the adjoint representation of U(4N ) as and T r stands for the normalized trace T r1 4N = 1. For future reference we write out the potential part of L YM separately as S YM is also invariant under the global SO(5) rotations of X a . It can be obtained from the dimensional reduction of the U(4N ) gauge theory in 5 + 1-dimensions to 0 + 1-dimensions, where the SO(5, 1) Lorentz symmetry of the latter yields to the global SO(5) of the reduced theory. There are two distinct deformations of S YM preserving its U(4N ) gauge and the SO(5) global symmetries. One of these is obtained by adding a fifth order Chern-Simons term to S YM (i.e. a Myers like term) which is given as [26] while the other is a massive deformation term of the form

JHEP12(2018)015
It is useful to write out the potential terms for S 1 = S YM + S Mass and S 2 = S YM + S CS explicitly: (2.7) In this paper, our main interest is on the massive deformations and hence we will focus on the dynamics emerging from S 1 , but will, although very briefly, also consider the consequences of some of the developments presented in the paper for S 3 = S YM + S Mass + S CS .
We may as well note that S 1 , S 2 and S 3 may be thought as deformations of a subsector of the bosonic part of the BFSS [1] matrix quantum mechanics. As is already well-known, BFSS model can be conceived to emerge from the dimensional reduction of the YM theory in 9 + 1 dimensions to 0 + 1-dimensions [4] with the SO(9, 1) of the former yielding to the global SO(9) of BFSS on the nine matrices X I (I : 1 . . . 9), while the latter may be further broken to SO(5) × SO(4), via the addition of S Mass and/or S CS terms. To be more precise, these deformations terms, involving only X a (a : 1, . . . 5), spontaneously break SO(9) down to SO(5) × SO(4) and naturally split the X I to an SO(5) vector X a and a SO(4) vector X α (α : 1 . . . 4). Then, S 1 , S 2 and S 3 emerges by focusing on the sector of X a 's only. V 2 is extremized by the matrices fulfilling This equation has two immediate solutions [26], one of which is the diagonal matrices while the other is given by a fuzzy four sphere S 4 F , where λ is forced to take on a fixed value λ = 2 n+2 , which depends on the matrix level of S 4 F . This latter fact requires attention, since it implies that the direct sum of fuzzy S 4 solve (2.8) if and only if they are fuzzy spheres of the same matrix level. For the configuration (2.9) the potential take the value zero as is readily observed from S CS , whereas a simple calculation shows that [26] it, in fact, has a lower (negative) value for the fuzzy S 4 solutions. Thus at the classical level, fuzzy S 4 appears to be a more stable solution than the diagonal commuting matrices. Further numerical studies have, however, revealed that the fuzzy S 4 is not a minimum of V 2 , but instead a saddle point [30]. 3 As for the massive deformation the potential part of S 1 is extremized by the matrices fulfilling Fuzzy four spheres S 4 F and their direct sums (even from different matrix levels) are solutions of this equation for µ 2 = −8. In a recent article Steinacker [27] showed that quantum corrections in the pure YM 5-matrix model stabilizes the radius of the fuzzy four sphere.

JHEP12(2018)015
We will see that superficial instability implied by negativity of µ 2 does not actually lead to a problem when we consider the equivariant fluctuations in S 1 about the S 4 F backgrounds. The reason for this is essentially that the potential of the emergent equivariantly reduced action is bounded from below at any finite matrix level. We may as well interpret this outcome as being due to the fact that the equivariant parametrization of the fluctuations introduces topological fluxes through the S 4 F stabilizing its radius. We will further elucidate on these points in the sections 3 and 4. Non-trivial topological fluxes leaves its imprints as kink type solutions of the reduced action in Euclidean signature, as we will exhibit in section 5.
For the action S 3 = S YM + S CS + S Mass the potential is extremized by the matrices fulfilling [26] [ whose fuzzy four sphere solutions need to satisfy the relation λ = 8+µ 2 4(n+2) , which includes the previous two cases of S 1 and S 2 as µ 2 = −8 and µ 2 = 0, respectively. The equation of motion for S 3 isẌ in the gauge A = 0 and subject to the Gauss law constraint a [X a ,Ẋ a ] = 0. Equations of motion for S 1 and S 2 are readily inferred from (2.12) by the remark following (2.11).
Being static, S 4 F configurations satisfying any one of the equations (2.8), (2.10), (2.11) also satisfy the corresponding equation of motions as well as the Gauss law constraint.

Models in the Euclidean signature
Wick rotating to the Euclidean signature we make the changes, t → −iτ , ∂ t → i∂ τ , A → iA, D t → iD τ and L → −L. Euclidean action (in 1 + 0-dimensions) is then (2.13) In section 5, we will consider the kink-type solutions, that is to say, the instantons of the low energy effective actions in 1 + 0-dimensions, which we obtain from the equivariant reduction of S E 1 .

Brief review of fuzzy S 4
In this subsection we collect some of the main features of the fuzzy S 4 construction [4,25,26,[31][32][33]. To start with we note that S 4 is embedded in R 5 as (2.14) Construction of fuzzy S 4 proceeds as follows. Let us denote by Γ a , (a : 1, · · · , 5) the Hermitian 4 × 4 gamma matrices associated to SO(5) fulfilling the defining anticommutation relations

JHEP12(2018)015
For concreteness we take them to be given in the form We introduce the n-fold tensor product acting on the n-fold completely symmetrized tensor product space which is the carrier space of the (0, n) IRR 4 of SO (5). Obviously the latter is equivalent to the completely symmetric tensor product Sym n (0, 1) of the fundamental 4-dimensional spinor representation of SO(5) acting on C 4 . Dimension of this representation and hence that of H n is given as X a are then N × N Hermitian matrices satisfying the relations which are the defining relations for the fuzzy four sphere, S 4 F . In fact (2.20a) may be seen as the SO(5) invariant condition giving the radius of S 4 F as r n := n(n + 4). This construction appears to be quite analogous to that of the fuzzy two sphere [34] with SO(3) of the latter replaced with SO(5), nevertheless only to the extent until one recognizes that the commutation relations of X a do not close but instead they are given by where G ab are the ten generators of SO(5) in its (0, n) IRR satisfying the commutation relations, G ab are anti-hermitian by the definition (2.21). Under the SO(5) transformations generated by G ab , X a transform as vectors (i.e. in the (1, 0) IRR of SO(5)) of SO (5) since

JHEP12(2018)015
We have that It is also useful note that, G ab and X a together generate the SO(6) in its (n, 0, 0) IRR. This structure may be compactly expressed by writing the 15 generators of SO(6) as with the commutation relations taking the usual form We now observe that the SO(5) invariant condition (2.20a) can be expressed as the difference of the quadratic Casimir operators: as (n, 0, 0) of SO(6) branches solely to (0, n) of SO(5) and hence they are of the same dimension N . The relation (2.20b) can also be expressed equivalently in the form (2.28) Another noteworthy feature of S 4 F is that there is a S 2 F attaching to each point of S 4 F . In other words, there is a S 2 F bundle over S 4 F with the total space being CP 3 F . This fact is reflected in the fields defined on S 4 F carrying an intrinsic spin, whose rank can be up to n [26,32].
We may also record that the commutative limit is achieved by taking n → ∞. The intricate structure of S 4 F with S 2 F fibers leads to where x a , with x a x a = 1, are the coordinates of S 4 ⊂ R 5 and ω ab is antisymmetric in its indices, satisfy x a ω ab = 0, as seen by taking the commutative limit of the first equality in (2.28) and generate, in fact, the S 2 in the fibration S 2 → S 4 → CP 3 . Detailed discussion on this may be found in [35].
In the commutative limit, adjoint action of X a and G ab become the differential operators [26] ad where the derivatives with respect to ω ab 's are given via We may note for future use that

Symmetry constraints and parametrization of the fields
We now focus on the YM model with mass term, whose action is given by S 1 The potential V 1 has an extremum given by four concentric S 4 F , satisfying the equation (2.10). We observe that U(4N ) gauge symmetry of the action S 1 is broken down to U(N ) × U(4), and the commutant of (3.2) is just U(4). Let us denote with F a the fluctuations about (2.10). We may write where the r.h.s. is introduced as a self-evident short hand notation, which will be used in the ensuing developments. We are interested in finding the equivariant fluctuations about the configuration (3.2). 5 To be more precise, we would like to concentrate on those F a , which are invariant under the SO(5) rotations of S 4 F up to SU(4) ⊂ U(4) gauge transformations. For this purpose we proceed as follows. We introduce the equivariant symmetry generators as are the generators of SO(5) in the 4-dimensional fundamental spinor representation (0, 1). They constitute the subset of generators Σ AB : in the fundamental (1, 0, 0) spinor representation of SO (6). Evidently, W ab satisfies the SO(5) commutation relations and its SO(5) representation content is given by the decomposition of the tensor product (0, n) ⊗ (0, 1) into a direct sum of IRRs as

JHEP12(2018)015
We digress a moment to note that, this structure can be lifted to SO(6) group by writing W AB = (W ab , W a6 ) with the representation content (n, 0, 0) ⊗ (1, 0, 0) ≡ (n + 1, 0, 0) ⊕ (n − 1, 1, 0) , (3.6) whose branching under the SO(5) simply yields (3.5). Same branching under SO(5) also holds for the complex conjugate representation (0, 0, n)⊗(0, 0, 1) ≡ (0, 0, n+1)⊕(0, 1, n−1). Let us examine the adjoint action of W ab . It consists of two terms, first of which generates the infinitesimal SO(5) rotations of the S 4 F , while the second term is responsible for generating the SO(6) gauge transformations in SO (5). Adjoint action of W ab carries the tensor product of the reducible representation (3.5) with itself and using general formulas on tensor product of SO(5) IRRs, it has the following decomposition in terms of SO(5) IRRs with the coefficients in bold indicating the multiplicities of the respective IRRs. For the case of n = 1, the decomposition takes the form To study the equivariant fluctuations about the configuration (3.2), we impose the symmetry constraints, that is the equivariance conditions, first of which means that the gauge field A is simply required to transform as a scalar of SO(5) under the adjoint action of W ab , which is quite natural as it does not carry any SO (5) index, while the second requires that the fluctuations F a introduced in (3.3) transform as a vector of SO(5) to comply with the equivariance condition. We infer from the decomposition (3.7) that the space of rotational invariants that may be constructed form the intertwiners of (0, n) and (0, 1) IRRs of SO (5) is threedimensional, while the space of vectors that may be constructed in term of the intertwiners and X a is of dimension seven. The aforementioned intertwiners may be introduced via the projection operators to the IRRs appearing in the r.h.s. of (3.5). We may express these three projections as where the factors of two in front of Casimirs are due to the unrestricted sum over a's and b's. P I are projections to the IRRs of SO(5) in the order given in the r.h.s. of (3.5) and C 2 (λ I ) = − 5 a<b M ab M ab stand for the quadratic Casimirs of SO(5) in the IRRs labeled by λ I ≡ ((0, n + 1), (1, n − 1), (0, n − 1)). Using (3.10) and the fact that idempotents can JHEP12(2018)015 be given as Q I = 1 4N − 2P I we may list the intertwiners as the idempotents By construction we do have Q 2 Let us also note that Q I are not all independent from each other as we have I Q I = −1 4N . A straightforward, but a long calculation, whose details are provided in appendix A, gives and will be useful in what follows. Adjoint representation of SO (6) branches under SO(5) as 15 → 5 ⊕ 10, or in the Dynkin notation: ( Thus, further insight on how SU(4) ≈ SO(6) generators sits in these intertwiners is gained by observing that Q I contain, ten of these generators as Σ ab , and the remaining five as Γ a , as seen from (3.12). Using (2.29), commutative limit of (3.12) takes the form (3.14) Consequently, we find for q I := lim n→∞ Q I : We may argue that, equivariant parametrization of the fluctuations introduces topological fluxes through the concentric S 4 F s, preventing the latter to shrink to zero radius. Without going into any technicalities regarding the S 2 fibre coordinates ω ab over S 4 , using (3.15b), this reasoning is supported by the fact that in the commutative limit the topological flux piercing the S 4 F may be linked to the second Chern number on S 4 : where p 2 = 1 2 (1 − q 2 ) stand for the rank four projectors generating the projective module over the algebra of functions over S 4 .

JHEP12(2018)015
We may solve the constraints given in (3.9a) and (3.9b) as follows. To satisfy (3.9a), we may choose to parameterize the gauge field A as where we have introduced α µ ≡ α µ (t) are real functions of time only, and eliminated Q 2 in favor of 1 4N using I Q I = −1 4N . From this form of the gauge field, it is readily observed that the U(4) gauge symmetry is broken down to U(1) × U(1) × U(1). However, later on we will see that term proportional to identity matrix in (3.17) decouples after the dimensional reduction and the gauge symmetry of the reduced actions will eventually be U(1) × U(1). For the fluctuations satisfying (3.9b), a convenient parameterization that befits the ensuing developments turns out be where we have introduced φ µ = φ µ (t) and χ ν = χ ν (t) as real functions of time only and used the notationX The 1 n factors appearing in the last three terms of F a via,X a andΓ a are naturally expected for the convergence of F a . Similar analysis performed in [36][37][38][39] on equivariant parametrization of fluctuations over S 2 F and S 2 F × S 2 F , shares the same features. Indeed, as n → ∞, we find (3.20) Demanding the fluctuations f a to be tangential to S 4 means that x a f a = 0 has to be fulfilled. Since x a ∇ a = 0, as noted after (2.31), the latter condition is satisfied if and only if φ 3 (t), χ 3 (t) and φ 4 (t) all vanish in this limit.

Reduction over S 4 F and the low energy effective action
We are now in a position to exploit the explicit parameterizations given (3.17) and (3.18) to determine the low energy effective action that emerges from the action S 1 by tracing over the concentric S 4 F 's. After a straightforward but a long calculation, with some intermediate steps relegated to the appendix B, we may write down the covariant derivatives in the form where we have used for the fields φ i (t) and χ i (t) for (i : 1, 2) only.
Although (3.23) does not appear to be manifestly positive definite, a simple analysis using Mathematica confirms that it is so, as it should be by construction. Thus, it is possible and useful to make a linear field redefinition in the sector spanned by φ 3 , φ 4 and χ 3 , i.e. convert to a basis, in which the kinetic term (3.23) is diagonalized. In the next section, we will naturally work with such linearly redefined fields, which diagonalize the kinetic term for specific values of n from n = 1 to n = 5.

(3.24)
A few remarks are in order. Firstly, we see that there are terms which are linear in the fields φ 3 , φ 4 and χ 3 . For finite values of n, which is essentially going to be our main focus in the next section, these terms cause no harm. In the n → ∞ the coefficients of these terms diverge linearly with n. Thus, alluding to our previous remark following (3.20), for the finiteness of the large n limit, it will suffice to assume that φ 3 , φ 4 and χ 3 vanish faster than 1 n . The mass terms then converges to −2µ 2 (|φ| 2 + |χ| 2 ), while the kinetic term is given by |D t φ| 2 + |D t χ| 2 in this limit. Let us also note that, we have µ 2 = −8, since we are inspecting this term about the S 4 F configuration satisfying (2.10). The exact form of the constant term C(n) is immaterial; in the next section we will adjust the overall constant term in the reduced Lagrangians so as to set the minimum value of the potential to zero.
Analytic calculation of the trace of the interaction term 1 4 T r[X a , X b ] 2 in the action S 1 for the equivariant parametrizations (3.17) and (3.18) turns out to be quite a formidable task as it involves large number of rather complicated traces. As an alternative approach, we have evaluated the traces for this term for the values of n = 1, 2, 3, 4, 5 using Mathematica. These already correspond to reasonably large span of matrix sizes 4N = 40, 80, 140, 224, respectively and gives us ample information for exploring the dynamics of the low energy reduced action. This is what we take up in the next section.

Gauge symmetry & the Gauss law constraint
The equivariantly reduced action obtained from S 1 is invariant under the U(1)×U(1) gauge transformations with the remaining fields φ 3 , φ 4 and χ 3 being real and thus uncharged under U(1) × U(1). We observe this manifestly from (3.23) and the mass term (3.24). The interaction term has the same gauge symmetry too by construction and can be manifestly seen from (6.1) at the level n = 1.
As the time derivatives of U(1) gauge fields α 1 (t) and α 3 (t) appear nowhere in the action, these fields have no dynamics of their own. Their equations of motion give us the Gauss law constraints: 1 2i (4.2)

JHEP12(2018)015
We make the gauge choice α 1 (t) = 0 and α 3 (t) = 0, which essentially amounts to the reality conditions φ * = φ and χ * = χ. We can be more precise by first noting that the Gauss law constraints does not break the U(1) × U(1) gauge symmetry completely, but a residual Z 2 ×Z 2 remains. Indeed, writing φ ≡ (φ 1 , φ 2 ) = |φ|(cos θ, sin θ) and χ ≡ (χ 1 , χ 2 ) = |χ|(cos σ, sin σ), we may express the constraints in the form Therefore, the remaining Z 2 × Z 2 symmetry is encoded in the gauge functions as Λ 1 (t) = Λ 0 1 + πk and Λ 3 (t) = Λ 0 3 + πk, where Λ 0 1 and Λ 0 3 are constants and k ∈ Z 2 . This indicates that, for both of the gauge functions, Λ 1 and Λ 3 , we have more generally Due to (4.3), we have θ(t) = θ 0 + πk and σ(t) = σ 0 + πk, and (4.4) holds for both θ(t) and σ(t), as well. Having noted these points, in what follows we set φ 2 (t) and χ(t) to zero (i.e., we have both θ 0 and σ 0 set to zero). Then, the Z 2 × Z 2 symmetry is implemented by In section 5 we consider the structure of the LEAs in the Euclidean time τ . Due to the Z 2 × Z 2 symmetry, we will be able to explore kink type solutions of the LEAs by imposing appropriate boundary conditions on φ 1 (τ ) and χ 1 (τ ) as τ → ±∞. Availability to impose topologically non-trivial boundary conditions on the latter can be attributed to the property (4.4) of the restricted gauge functions, which holds the same in the Euclidean signature.
There is also the possibility of not imposing the Gauss law constraint as recently discussed for BFSS and BMN matrix models in [28]. This leads to presence of Goldstone bosons for the LEAs that we have obtained, as we will briefly discuss and demonstrate in section 6.

Dynamics of the reduced action and chaos
The explicit form of the equivariantly reduced Lagrangian for n = 2 is given below while the explicit form of L (n) for n = 3, 4, 5 are given in appendix B. Let us note that in L (n) for n = 2, 3, 4, 5, we have i) performed the linear transformation among the fields φ 3 → φ 3 , φ 4 → φ 4 , χ 3 → χ 3 which diagonalizes the kinetic term (3.23), and dropped the 's in the final form, ii) have set µ 2 = −8 as we have already remarked to do so in the paragraph after (3.24), iii) have imposed the Gauss law constraints as we just discussed, iv) adjusted the constant in the final form of each L (n) , so that the associated potentials, V (n) , take the value zero at their minima, v) introduced an over-dot,(. . .), to denote the time derivatives and vi) set the coupling constant g to one, as it has no effect on the classical physics save for determining a global normalization in the energy unit.
For n = 1, the reduced action takes a simpler form as compared to the cases for n ≥ 2, which is given as where we have introduced Φ = φ 3 +φ 4 −χ 3 , which is the only combination of the constituent dynamical variables φ 3 , φ 4 and χ 3 upon which L (n=1) depends. This is to be expected in view of (3.8). For L (n=1) too, all items following (4.5) are performed as well, except the item i), which is already taken care of with the introduction of Φ(t). An important feature of the reduced Lagrangians is that their potentials are all bounded from below. Therefore, we infer that at any level n the equivariant fluctuations about the concentric S 4 F solutions, with µ 2 = −8 do not cause any instability. The absolute minima of the potentials are given in (B.7), (5.6), (B.8) and (B.9).
We find that the reduced Lagrangians, L (n) , have chaotic dynamics. One of the basic tools to probe the presence of chaos in a dynamical system is to compute the Lyapunov exponents, which measures the exponential growth in perturbations [44]. If, say, x(t) is a phase space coordinate, in a chaotic system the perturbation in x(t), denoted by δx(t), deviates exponentially from its initial value at t = 0; |δx(t)| = |δx(0)|e λ L t , λ L being the Lyapunov exponent corresponding to the phase space variable x(t).
The phase space corresponding to the LEA are 10-dimensional, except for the n = 1 case, and spanned by where p i are the corresponding conjugate momenta and the Hamiltonians, H (n) , are obtained from L (n) in the usual manner using H = p iqi − L. To be more precise, we ran a Matlab code, which calculated the mean of the time series for 20 runs with randomly selected initial conditions for all of the Lyapunov exponents at each n and for the energies given above. The initial conditions are randomly selected by the code from a sector of the 10-dimensional phase space for (n = 2, 3, 4, 5). The latter is specified by giving the initial values of the eight of the phase space coordinates, while the code randomly selects an initial value for one phase space coordinate and calculates the last one to satisfy the given value of the energy. We have checked that, increasing the number of the randomly selected coordinates somewhat increases the computation time, but does not have any significant impact on our results. For n = 1, dimension of the phase space reduces to 6 as easily observed from L (1) in (4.6). A similar analysis to the one described above is also performed at this level, whose marked differences and similarities with the rest will also be pointed out in what ensues. Table 1 summarizes our findings for the largest Lyapunov exponents λ max at (n = 2, 3, 4, 5) at the listed values of the energy. Table 2 and 3 give the LLE data at a larger set of energies for n = 1 to probe especially the low energy region, which appears to have different features. It is especially interesting from our data to explore the dependence of the LLEs at a given level n with respect to the energy. We find that our data for LLEs fit very well with JHEP12(2018)015     Table 4. α n & β n values for the fit in (4.9). the functional relation λ n (E) = α n + β n 1 √ E (4.9) where λ n (E) denotes the LLE as a function of energy at fixed n. The plots for the data and corresponding fits are given in figures 1, 2, 3, and 4. We have also calculated the standard error for the largest Lyapunov exponents from the standard deviation of the final mean value of λ n (E) using the Largest Lyapunov exponents from each of the 20 runs at each level n and at the energies listed above. The errors are quite small, typically remain around 0.0050 with the span being from ±0.0018 to ±0.011 and thus appear very small in the figures. Our findings appear to be quite novel in the sense that, to the best of our knowledge, they constitute the only result in the literature within the context of matrix Yang-Mills theories at zero temperature, in which the dependence of the LLEs on the energy is predicted from numerical data at several of the lowest lying matrix levels. The data and the JHEP12(2018)015  At the level n = 1 we find a markedly different behaviour of LLEs with increasing energy for E ≤ 50 as can be seen from the data and corresponding plots given in figures 5 and 6. Sudden decrease in the value of the Lyapunov exponents around E = 40 can be attributed to the fact that the potential V (1) takes the value 40 at its local maximum. Although the value of LLE decreases around this nonstable equilibrium point it resumes to attain the growing profile for E ≥ 45 and the same functional form which befit those of n ≥ 2 also fits well with the numerical results as can be seen from figure 5. We have found the suitable fit to be λ 1 (E) = 1.78 − 12.9 1 √ E for E ≥ 45. These results also enable us to probe the onset of chaos in our models. To do so, let us first remark that the numerical calculations are performed for a finite computation time only, therefore even at low energies numerical values of the LLEs are small numbers compared to the characteristic scale of LLE values within a given model, but may not be seen to vanish in finite time. Secondly, it is well-known that in Hamiltonian systems periodic dynamics and chaotic dynamics can coexist and that there is in general no sharp passage from one to another [45]. Keeping these two facts in mind, it thus appears reasonable to set a critical lower bound on the LLE value at and above which the models have appreciable amount of chaotic dynamics, (i.e. there are comparable number of chaotic and periodic trajectories.) Inspecting the data from table 1 and the time series plots of the Lyapunov spectrum given in the appendix B.4, we see that a reasonable choice would be to take this bound to be about 0.1. Then, we find that for n = 2, 3, 4, 5, critical energies for the onset of chaos turns out to be E ≈ 16, 22, 23, 24, respectively, with the corresponding LLE values being λ 2 (E = 16) = 0.1167, λ 3 (E = 22) = 0.1039, λ 4 (E = 23) = 0.0996, λ 5 (E = 24) = 0.1178. Lyapunov exponents become smaller below these energies as can be inspected from table 1. In fact, using equation (4.9), we predict that they get smaller at an increasing rate of dλn(E) dE = − 1 2 β n E − 3 2 , (β n < 0), as E decreases. From the fits we find that for n = 2, 3, 4, 5, LLEs vanish at the energies ≈ 12.6, 16, 17.4, 17.3, respectively, which JHEP12(2018)015  Table 5. γ E and δ E fit values.
is consistent with the small values obtained for LLE at E = 15. Below E = 15, part of the initial conditions become too small for the numerical integrator built into Matlab to handle and we can not obtain any healthy data on the LLE in this region. However, there appears no reason to expect that any significant amount of chaos remains below E = 15.
It is also worthwhile to explore the change in LLE values as n takes on the values n = 2, 3, 4, 5 at fixed value of the energy. Figure 7 depicts this at several different values of the energy. It turns out that logaritmic functions of the form provide a good fit to the data. Here λ E (n) stand for LLE as a function of n at fixed energy. The dashed lines in figure 7 are provided for visual guidance as n takes on only the integer values. Values of γ E , δ E are also provided below for convenience. Only important feature we infer from these fits is that at a given energy E, decrease in λ E (n) appears to be logarithmically slow, suggesting that, even for n > 5, we may expect to have chaotic dynamics at moderate values of the energy.

JHEP12(2018)015
Before closing this section we want to note that the LEAs treated in this paper have only Z 2 × Z 2 symmetry after gauge fixing as already explained in section 4.1. This residual gauge symmetry is spontaneously broken by the degenerate vacuum configurations of zero energy, which are given section 5.2 and in appendix B.5 (the potential V (n) and hence the corresponding vacua are the same both in real time and Euclidean time), as well as by the instanton solutions given in the next section. We see that the configurations that respect the Z 2 × Z 2 symmetry are only those for which both φ 1 and χ 1 are vanishing, as these are the only two fields changing sign under the Z 2 × Z 2 action. We thus easily infer that those configurations respecting this Z 2 × Z 2 symmetry and possessing the smallest energy are the static solutions of the equations of motion minimizing V (n) (φ 1 = 0, χ 1 = 0, φ 3 , φ 3 , χ 3 ). Using Mathematica we have calculated that these configurations have the energies E = 10, 15.1, 18.0, 19.7, 20.8, respectively for n = 1, 2, 3, 4, 5. At each of n = 2, 3, 4, 5 there are eight distinct configurations at the given energy, while at n = 1 there are only two.

Kink solutions
We now consider the matrix model in the Euclidean signature S E 1 given in (2.13) and the corresponding LEAs. The latter have multiple degenerate vacua supporting kink solutions, i.e. instantons in 1 + 0-dimensions, whose features are sketched out in what follows.

Kinks at level n = 1
The Lagrangian is given as where stands for derivatives with respect to the Euclidean time τ . We have also scaled out an unimportant factor of 5 2 compared to (4.6). We can easily see from (5.1) that there are three different pairs of vacua, which are given by the configurations Since either φ 1 or χ 1 vanish in these vacua, we infer that, the kink solutions could be of the type with topological indices (±1, 0) or (0, ±1) ∈ Z 2 ⊕ Z 2 . These are the familiar kink solutions [46,47]. Indeed, we find that the equations of motion are of the form

JHEP12(2018)015
which have the following solutions

Kinks at levels n ≥ 2
For n = 2, 3, 4, 5, the number of degenerate vacua increases. This may be expected due to the larger number degrees of freedom in the LEAs. Nevertheless, a similar structure in vacuum configurations to that of n = 1 is observed, and allow for the kink solutions. At n = 3, for instance, we have eight pairs of degenerate vacua, which are given as The equations of motion for L (n=3) are coupled non-linear differential equations, which are not easy to solve. We may look at the linearized system of equations around one of the minima. For notational simplicity, let us write (φ 1 , χ 1 , φ 3 , φ 4 , χ 3 ) ≡ (S 1 , S 2 , S 3 , S 4 , S 5 ) := S and also introduce S = S 0 + s, where S 0 is one of the vacuum configurations in (5.6) and s stands for the fluctuations. The linearized system of equations is given by where no sum over repeated indices is implied on th l.h.s. and C i can be easily read off from (B.4). For, say, S 0

JHEP12(2018)015
The leading order profiles of the solutions of these equations which are regular as τ → ∞ are given below, while their complete forms are given in appendix B. We have where c i (i : 1 · · · 5) are arbitrary constants. (5.9) provides the asymptotic profile of the kink solutions near S 0 + .

Kinks in the presence of the Chern Simons term
Here we will confine our discussion only to the level n = 1. In the Euclidean signature, the reduced action obtained from S 3 takes the form (5.10) Last term in (5.10) is the contribution coming from the reduction of S CS as can be clearly observed from the fifth order terms it contains. Due to presence of this term the potential of (5.10) does not have an absolute minimum, in fact, naturally, it is not bounded from below. Nevertheless, for µ 2 < 0, it is a matter of a simple calculation to see that the potential still have degenerate local minima. This still allows for kink solutions, which are, however, only metastable and will decay under sufficiently large perturbations. As a concrete example, we have, for instance for µ 2 = −1, the local minima occurring at and a kink solution to the equations of (5.10) is given by

Gauge symmetry revisited
In [28] Maldacena and Milekhin considered the BFSS model without imposing the Gauss law constraint, i.e. without the SU(N )-singlet condition on the physical states. This is based on the fact that the BFSS model with A = 0 is still well-defined even in the absence of the Gauss law constraint, at the expense that the SU(N ) is no longer a local but a global gauge symmetry group in this situation.

JHEP12(2018)015
We observe that such a possibility is also valid and applicable to the low energy reduced actions that we have obtained in this paper. The latter have U(1) × U(1) gauge symmetry and the Gauss law constraint was given in (4.2), imposing the U(1) × U(1) singlet condition meaning that the complex fields φ(t) and χ(t) are uncharged, i.e. real, under the gauge fields α 1 (t) and α 3 (t), respectively. If we do not impose the constraint we can still set the gauge fields α 1 (t) and α 3 (t) to zero in the LEAs at any level n. In this case, the LEAs have a global U(1) × U(1) symmetry, which is spontaneously broken by several different vacuum configurations, and hence imply the existence of Goldstone bosons in these LEAs. For instance, at the level n = 1, we have the action with three different vacuum configurations, as easily recognized from (5.2) spontaneously breaking the U(1) × U(1) symmetry. Thus, we immediately infer that the fluctuations around each of these vacuum configurations give rise to one Goldstone boson in the usual manner that it arises in an abelian gauge theory with degenerate vacua. As a concrete example, we have, with φ = 2 + σ 1 + iσ 2 , Φ = − 1 2 + ρ and denoting the fluctuation around χ = 0 still with χ, potential part of L (n=1) takes the form V (n=1) (σ 1 , σ 2 , χ, ρ) = 1 4 (σ 2 1 + σ 2 2 + 2σ 1 + |χ| 2 ) 2 + 3 4 (σ 2 1 + σ 2 2 + 2σ 1 + 4ρ 2 ) 2 + (2 + σ 1 ) 2 + σ 2 (|χ| 2 + 6ρ) , showing that σ 2 is massless, i.e. it is the Goldstone boson associated with this particular vacuum configuration. Similar analysis show the existence of Goldstone bosons for the other two vacuum configurations in (6.2). Finally, we note that for n = 2, 3, 4, 5 we infer from (B.7), (5.6), (B.8) and (B.9) that, at each value of n, there are eight distinct vacuum configurations determined by the values of φ, χ and the real fields φ 3 , χ 3 and φ 4 , each of which comes with a Goldstone boson.

Conclusions
In this paper we have studied the 5d mass-deformed Yang-Mills matrix model with U(4N ) gauge symmetry. We have found the exact form of the SU(4) ≈ SO(6) equivariant parametrizations of the gauge field and the fluctuations about the classical four concentric fuzzy four sphere configuration and used them to calculate the LEAs by performing traces over the S 4 F s for the first five lowest matrix levels. The LEA's obtained in this manner have potentials bounded from below, which indicates that the equivariant fluctuations about the S 4 F configurations with a tachyonic mass term (µ 2 = −8) do not lead to any instabilities.

JHEP12(2018)015
We have showed through detailed numerical computations that these reduced systems have chaotic dynamics and exhibited its various features. In particular, based on the numerical calculations of the Lyapunov spectrum we deterined how the LLE behaves as a function of energy, and also were able to comment on the aspects of the onset of chaos in these models.
In the Euclidean signature, we have demonstrated that the LEAs support the usual kink type solutions, i.e. instantons in 1 + 0-dimensions. The latter may be viewed as the residual topologically non-trivial configurations, linked to the topological fluxes penetrating the concentric S 4 F s due to the equivariance conditions, and preventing them to shrink to zero radius.

Acknowledgments
We thank O. Oktay for discussions and sharing his computer code with us for the calculation of the Lyapunov exponents. We also thank K. Baskan for his help on our computer code at various instances. S.K. thanks H. Steinacker for discussions during the annual workshop of COST Action MP1405 held in Sofia, Bulgaria, February 2018. The work of S.K. and G.C.Toga is supported by the Middle East Technical University under Project GAP-105-2018-2809.

JHEP12(2018)015
A.2 Computation of (G · Σ) 2 Here we present the details of the calculation leading to (3.12). We have Multiplying both sides of (A.6) with abcdm and using ijklm γ i γ j γ k γ l = 24γ m as can be inferred from (2.20b) one obtains To handle the left hand side of (A.7) we make use of the determinant Left hand side of (A.7) then reads where the second line follows after making use of {Γ a , Γ b } = 2δ ab for rearrangements. Therefore, (A.7) can be cast into the form Employing G ab = 1 2 [X a , X b ], we can expand the first term in r.h.s. of (A.2) as after using (2.20b). Substituting G ab G ab = −4n(n + 4)1 4N for the second term and simplifying the third term as and finally, combining all the terms together we get (G · Σ) 2 = 12 Γ a Γ b G ab + 8(n + 2)X a Γ a + 8n(n + 4)1 4N . (A.12)

JHEP12(2018)015
B Details on the dimensional reduction of S 1

B.1 Useful identities
Some useful identities among Q 1 , Q 3 and X a , which simplify the analytic calculations are listed below:  .17) and (3.18) in D t X a = ∂ t X a − i[A, X a ], we find

(B.2)
With the help of the identities listed in (B.1), this simplifies to                                    The absolute minima of the potentials V (n) associated to L (n) are given below.