Chaotic spin chains in AdS/CFT

We consider the spectrum of anomalous dimensions in planar $\mathcal{N}=4$ supersymmetric Yang-Mills theory and its $\mathcal{N}=1$ super-conformal Leigh-Strassler deformations. The two-loop truncation of the integrable $\mathcal{N}=4$ dilatation operator in the SU$(2)$ sector, which is a next-to-nearest-neighbour deformation of the XXX spin chain, is not strictly integrable at finite coupling and we show that it indeed has Wigner-Dyson level statistics. However, we find that it is only weakly chaotic in the sense that the cross-over to chaotic dynamics is slower than for generic chaotic systems. For the Leigh-Strassler deformed theory with generic parameters, we show that the one-loop dilatation operator in the SU$(3)$ sector is chaotic, with a spectrum that is well described by GUE Random Matrix Theory. For the imaginary-$\beta$ deformation, the statistics are GOE and the transition from the integrable limit is that of a generic system. This provides a weak-coupling analogue of the chaotic dynamics seen for classical strings in the dual background. We further study the spin chains in the semi-classical limit described by generalised Landau-Lifshitz models, which are also known to describe large-angular-momentum string solutions in the dual theory. We show that for the higher-derivative theory following from the two-loop $\mathcal{N}=4$ SU$(2)$ spin chain, the maximal Lyapunov exponent is close to zero, consistent with the absence of chaotic dynamics. For the imaginary-$\beta$ SU$(3)$ theory, the resulting Landau-Lifshitz model has classically chaotic dynamics at finite values of the deformation parameter.


Introduction
The realisation that in planar N = 4 supersymmetric Yang-Mills (SYM) theory, the dilatation operator acting on single-trace operators can be mapped to a long-range integrable spin chain [1][2][3], has led to an impressive range of non-perturbative results, see [4][5][6][7] for reviews. These results make fundamental use of the holographic duality [8] which relates the anomalous dimensions of singletrace operators to the energies of non-interacting closed strings in an AdS 5 ×S 5 background. While many of these results have been extended to other theories, such as to deformations of N = 4 SYM or to theories which share the feature of being highly supersymmetric like ABJM theory [9], such theories are very special and in general one will not have at hand the tools provided by integrability.
For classical strings moving in more general geometries there is in fact significant evidence of chaotic dynamics. By computing the power spectrum, Lyapunov characteristic exponents and Poincaré sections, in [10] it was shown that the motion of classical circular strings in the AdS 5 -Schwarzschild background is chaotic. This was subsequently extended to a wide range of backgrounds such as the AdS soliton background [11], AdS 5 ×T 1,1 [12], other AdS 5 × Einstein 5 spaces [13] and p-brane backgrounds [14][15][16], see also . Moreover, the string motion in the gravity dual of β-deformed N = 4 SYM theory was shown to be chaotic for imaginary values of β [39]. In this work we provide a weak-coupling analogue of this classical string chaos by considering the spectral statistics of the perturbative β-deformed dilatation operator in the dual gauge theory.
On the gauge-theory side, chaotic behaviour in classical Yang-Mills theory has been long known, see [40][41][42] for reviews of early work. It was shown that upon restricting to homogeneous field configurations or by dimensionally reducing on S 3 , the resulting Hamiltonian has non-integrable dynamics. This can be related to the computation of anomalous dimensions as, by using a conformal transformation, the dilatation generator of a conformal theory defined on R 1,3 becomes the Hamiltonian for the theory on R × S 3 . In the case of N = 4 SYM, reducing the theory on S 3 produces the BMN matrix theory [43] which is known to have classically chaotic dynamics [44]. More generally, Yang-Mills is a theory of many interacting degrees of freedom and so we expect it to demonstrate the property of quantum many-body chaos. Spectral statistics can provide a key signature of quantum chaos in many-body systems, see [45,46] for reviews. In particular, the BGS conjecture [47] states that quantum-chaotic systems have spectra whose fluctuations are described by Random Matrix Theory (RMT). An analysis of the spectrum of the planar dilatation operator of N = 4 SU(N ) SYM theory expanded in the 't Hooft coupling λ was performed in [48] and [49], and it was shown that, for rank-one sectors of the one-loop planar theory, the spacing between adjacent levels is described by the Poisson distribution as is characteristic of integrable systems [50]. However, when non-planar corrections are included, the spectrum becomes chaotic and the level spacings closely follow the Wigner-Dyson distribution for the Gaussian Orthogonal Ensemble (GOE). Additional analysis of spectral observables and multi-trace eigenstates provided further evidence that, at finite values of N , the theory is described by GOE RMT. For operators whose dimensions scale as N 2 , and so are sensitive to non-planar effects, the theory has also been shown to exhibit fast scrambling [51].
In this work, we examine a different set of questions and consider quantum-chaotic behaviour that also occurs at weak coupling but in the planar limit. First we consider the two-loop planar dilatation operator of N = 4 SYM theory restricted to a closed SU (2) sector. The theory is explicitly integrable at one-loop order, where it corresponds to the XXX 1/2 -Heisenberg spin chain, and is believed to be integrable also at finite coupling. At loop order k, the spin chain has interactions involving k + 1 adjacent spins and so the all-order chain has infinite interaction range. When truncated to just the two-loop terms involving next-to-nearest neighbour interactions, the spin chain is only perturbatively integrable and higher charges are only conserved to leading order in the coupling. We investigate the level statistics of this two-loop spin chain at finite coupling and find that it is indeed described by GOE RMT. The two-loop deformation is nonetheless non-generic and only gives rise to a weakly chaotic system in the sense that the transition between the integrable and chaotic regimes is slower and only occurs for larger values of the deformation than in the generic, strongly chaotic case. This is consistent with the work [52] which showed that deformations of the XXZ spin chain constructed via the boost method of [53,54] are weakly chaotic and the fact that the two-loop dilatation operator is similarly constructed. An interesting corollary of this analysis is that the SU(2) preserving, next-to-nearest deformation of the XXX spin chain is only weakly chaotic. We find similar results for twisted versions of the theory which describe a rank-one sector of the integrable β-deformed theory.
To understand quantum chaos in more generic planar theories, we consider the one-loop dilatation operator in the SU(3) sector of the marginally deformed N = 4 SYM theory preserving N = 1 supersymmetry [55]. This theory is parametrised by two additional complex parameters, h and q. There are specific values for which the theory is integrable [56], see also [57][58][59][60], and where the level statistics are Poisson. In general this is not the case and we show that for generic values the spectrum of the theory is described by Gaussian Unitary Ensemble (GUE) RMT. In the case h = 0 and q ∈ R, which corresponds to imaginary values of β as studied in [39], it is again GOE. By examining the transition from integrable points to the chaotic regime, we see that these deformations are strongly chaotic, thus suggesting a weak-coupling analogue of classical string chaos. Interestingly, quantum chaos with GOE statistics has previously been seen for strings in space-times dual to confining theories [61].
To make the connection between spin-chain quantum chaos and classical string chaos more manifest, we analyse the large-length limit of the spin chain. In this limit the low-energy states of the Heisenberg spin chain are known to have a semi-classical description in terms of a twodimensional Landau-Lifshitz (LL) type sigma model [62]. The same sigma model action can be found from a large angular momentum, or fast-string, limit of the AdS 5 ×S 5 world-sheet theory [63] which thus provides a bridge between the two sides of the duality. This identification has been extended to larger sectors [64][65][66][67][68], including a PSU(2, 2|4) LL model describing the complete oneloop N = 4 SYM dilatation generator [69]. The semi-classical limit of the higher-loop, longer-range Hamiltonians gives, after integrating out the short wavelength modes, a generalised LL action with higher-derivative interactions. The two-loop effective LL action with four-derivative terms was found from the spin chain [70] and matches the corresponding limit of the string theory. However, in general the LL action following from the spin chain and string theory disagree. The allowed terms in a generalised LL action can be combined with arbitrary coefficients which then parametrise the space of possible theories. The agreement at one-and two-loops between the gauge-theory LL model valid at weak coupling and the string LL model valid at strong coupling is presumably a result of a non-renormalisation theorem which fails for the λ 3 , six-derivative terms [71][72][73][74][75]. Nonetheless, even where the two models disagree, having a common sigma-model description is extremely useful for understanding the duality. We show that while the two-loop LL model with higher derivative interactions has more erratic behaviour, the dynamics do not appear to be strongly chaotic. In particular, we compute the Lyapunov characteristic exponents for a specific family of solutions, and show that they are not clearly larger than zero. This is consistent with the weak chaotic behaviour of the spin chain and the classical integrability of the string theory.
LL models for deformed backgrounds have also been studied and in [76] the LL model corresponding to the β-deformed theory, i.e. the marginally deformed theory with h = 0 and q = exp(iβ) with β ∈ R, was considered. In the gauge theory this corresponds to considering the long-wavelength limit of the twisted spin chain, while on the string side it describes fast strings moving in the dual geometry constructed in [77] by using the TsT-transformation. In this case the classical string theory is integrable [78] and there are all-order twisted Bethe equations for the spectrum [60], see also [79][80][81][82][83] for reviews. The more general dual background for q = exp(κ + iβ) was also constructed in [77] by taking additional S-duality transformations. In this case the geometry is not an exact string solution and there are likely α corrections. Moreover, the string theory is no longer integrable and as mentioned has been shown explicitly to be chaotic [39]. One can still take the fast string limit of the world-sheet theory and study the long-wavelength limit of the one-loop dilatation operator. This was done in [84] for a holomorphic SU(3) sector where, after a careful accounting of the relevant BPS states, the same deformed LL model was found in both cases. In this work we study the dynamics of this deformed LL model and show, by computing the Lyapunov exponents for specific solutions, that in case of real q = exp(κ) it has classically chaotic dynamics. This demonstrates that there is classical chaos in the semi-classical limit of the spin chain consistent with our quantum spin-chain results and also in agreement with the classical chaos found in the string theory. We further study the transition from integrable to chaotic dynamics as κ increases from zero. We find that, for all the configurations we consider, there are finite values of κ for which the largest Lyapunov exponent remains zero. Such behaviour is familiar in classically chaotic systems but is surprising in the quantum theory where we expect a sharp transition to chaotic dynamics in the infinite volume limit and is likely due to the fact that the limit considers only the low-energy states.

Weakly chaotic SU(2) spin chains
In N = 4 SYM theory, the planar dilatation operator acting on single-trace gauge-invariant operators composed of two types of complex scalars, X and Z, can be written as an SU(2) invariant spin-chain Hamiltonian. The map between gauge-theory operators and spin-chain states is Tr(ZXZZX . . . ) = |↑↓↑↑↓ . . .
(2.1) with each Z replaced by |↑ , each X by |↓ , and where the cyclicity of the trace must be imposed as an additional condition on the spin-chain state. This corresponds to restricting to states that are invariant under cyclic shifts of the lattice sites or, equivalently, that have zero lattice momentum. The perturbative dilatation operator acting on such operators at the origin can be computed as an expansion in the rescaled 't Hooft coupling g 2 , (2.2) In the planar limit, the spectrum of the dilatation operator can be reproduced by a spin-chain Hamiltonian with periodic boundary conditions which is written as a sum of terms each acting on contiguous sets of spin-chain sites where L corresponds to the length of the spin chain. The leading term, D (0) i = 1 i , is simply the identity operator and counts the number of fields/spins and gives the classical scaling dimension. The one-loop and two-loop terms [1,2] are is the permutation operator acting on the a-th and b-th sites. Writing the permutation operator in terms of the usual spin operators, S = 1 2 σ, we can write the one-and two-loop terms as a special case of the Heisenberg spin-chain Hamiltonian with next-to-nearest-neighbour (NNN) corrections with the choice It is believed, and there is extensive evidence, that planar N = 4 SYM theory is integrable at all orders in λ (see [4][5][6] for reviews). Specifically, the perturbative L → ∞ spin-chain spectrum can be reproduced by asymptotic Bethe equations which have a finite radius of convergence, |g| < 1/4 [85]. On the other hand, while the unperturbed XXX spin chain is integrable, it is known that the XXX spin chain with NNN interactions e.g. (2.6) with 0 = 0, λ 1 = 1 and λ 2 = δ, is chaotic [86,87].
Since the all-order planar dilatation operator is integrable, as is the one-loop term on its own, it is interesting to see the appearance of chaos when we truncate to include only the two-loop terms. In particular we are interested in the behaviour of the theory as we vary the ratio of the strength of the next-to-nearest-neighbour term to the nearest-neighbour term which controls the strength of integrability breaking. We can write this ratio in terms of the planar gauge coupling as where we see that δ increases monotonically from g 2 = 0 until g 2 = 1/4 where it becomes negative and then asymptotically approaches δ = −1/4. We thus expect that the system becomes chaotic as we increase the coupling g 2 from 0 to 1/4 and then become somewhat less chaotic. In order to explore this behaviour, we numerically compute the statistical properties of the energy levels in the following.

Level-spacing statistics
For quantum-chaotic systems it is generally believed 1 that, while the overall dependence of the spectral density on the energy is specific to a given system, the fluctuations of energy levels about this average trend are universal and described by Random Matrix Theory. This is often referred to as the BGS conjecture [47]. For RMT, the distribution of spacings, s, between adjacent levels is well-described by the Wigner-Dyson distribution where the parameter α depends on the Gaussian ensemble: α = 1 for the orthogonal, α = 2 for the unitary and α = 4 for the symplectic ensemble. The choice of the appropriate ensemble for a given Hamiltonian is determined by the presence or absence of discrete symmetries, and in particular the structure of the correlations between energy levels is influenced by time-reversal symmetry. In general the time-reversal operator is an anti-unitary operator which can be written in the form T = KC, where K is a particular unitary operator and C takes the complex conjugate of any operator or state it acts upon. It is known, see e.g. [88], that for systems which are timereversal symmetric and additionally rotationally invariant, there is a choice of states for which the Hamiltonian is real and symmetric and thus the appropriate ensemble is the orthogonal one. On the other hand, for integrable systems energy levels are essentially uncorrelated and it has been shown that in a variety of integrable systems, e.g. [50,87,86], the level-spacing distribution is Poisson Before computing the level-spacing distribution for various parameter configurations of the Hamiltonian (2.6)-(2.7), we must first de-symmetrise the spectrum by focussing on a subspace of states all of which have the same values of the global charges. States with different charges cannot mix and so the Hamiltonian is block diagonal, leaving eigenvalues in different sectors uncorrelated, and thus a level-statistics analysis is only meaningful within individual blocks. For the Hamiltonian (2.6)-(2.7) de-symmetrisation requires considering chains with fixed length L, and fixed total spin S z = L i=1 S z i or equivalently fixed impurity number M = L/2 − S z . As the Hamiltonian is translation invariant, we must also fix the total momentum P and to connect with the cyclicity of the trace defining gauge-invariant operators we choose P = 0. Furthermore, there is a parity charge corresponding to the symmetry transformation which reverses the ordering of spins in a spin-chain state, e.g. P : |↓↑↑↓↓↑ → |↑↓↓↑↑↓ , (2.12) and here we consider the larger sector of positive-parity states. Finally, as there is a global SU(2) symmetry, we only consider lowest-weight states satisfying S − |Ψ = 0. After computing the energy eigenvalues, {e i }, for a particular de-symmetrised subspace, we must unfold the spectrum to remove the overall energy dependence and fix the average spacing to be one. This we do by splitting the cumulative level number where m is the total number of states in the sector, into an average and fluctuation part n(e) = n av (e) + n fl (e) . (2.14) The average part is determined by first ordering the energy levels, then discarding some (usually around 5-10%) of the low-and high-energy levels, and then fitting the spectrum to a relatively high-order (usually degree 15) polynomial which defines n av . The unfolded spectrum is then given by and captures the physics of the spectral fluctuations. Finally, we compute the level spacings s i = ε i+1 − ε i , divide the spacings range into bins (usually [15][16][17][18][19][20] and approximate the distribution, P (s)ds, by counting the number of spacings in each bin. The details of the unfolding and binning procedure can in certain cases have significant effects, so we choose values where the results do not vary strongly for small changes in the parameters. We show the results of these numerical calculations for the sector of positive-parity primary states with L = 20 and M = 7 in Fig. 1, where the Poisson distribution for g 2 = 0 and the GOE Wigner-Dyson distribution for g 2 = 0.2 is manifest. More quantitatively, we can define a generalised Brody distribution, 17) which interpolates between the Wigner-Dyson distribution (2.10) for arbitrary α at ω = 1 and the Poisson distribution at ω = 0. For g 2 = 0 and assuming α = 1, we find the best fit is for a value of ω close to zero, in particular ω = 0.03, i.e. the nearest-neighbour levels are uncorrelated which is consistent with the integrability of the g 2 = 0 spin-chain model. For g 2 = 0.2 and assuming α = 1, we find the best fit is for a value of ω close to unity, ω = 0.88. If we assume ω = 1 instead, the best fit is for α = 0.82, and fitting both parameters at the same time we find (α = 1.1, ω = 0.81) which is close to the GOE values and consistent with a quantum-chaotic nature of this spin-chain system. We can formally continue the coupling beyond g 2 = 1/4 such that the effective coupling δ becomes negative and we see that in this case the distribution becomes intermediate between Poisson and Wigner-Dyson with a best fit of ω = 0.33 when assuming α = 1.
That we find a GOE spectrum is consistent with the symmetries of the Hamiltonian in (2.6). General spin-chain Hamiltonians with SU(2) symmetry and arbitrary interaction range [57,58] where K is the unitary operator which, up to a minus sign, simply flips each spin. This symmetry guarantees that the Hamiltonian is real and symmetric, which argues that the spectrum of the XXX spin chain with a long-range deformation of the type (2.18) should be described by the orthogonal ensemble. However, because the symmetry involves a spin flip, the action of the time-reversal operator will relate states in different sectors (except for the case where the M = L/2) and thus does not quite provide the explanation for GOE statistics within a sector of fixed M . This is important, as in order to see the Wigner-Dyson statistics, we worked within such sectors of fixed M . Fortunately, there is an additional symmetry: If we take K to be the spin-chain parity operator P which takes σ A i σ B i+s → σ B i σ A i+s , the combined operator PC acts within each of the fixed M sectors and leaves the spinchain Hamiltonian invariant. Now by working with eigenvectors of P, which is also necessary to fully desymmetrise the spectrum, we can choose a basis of states such that P = +1 and the Hamiltonian is real and symmetric when restricted to these sectors. This construction can be generalised to a time-reversal symmetry of the one-loop non-planar dilatation operator of N = 4 SYM theory by interpreting P as the operator which reverses the order of fields within individual traces of gauge-invariant operators and so explains the GOE statistics seen in [48,49].

Transition from integrable to chaotic dynamics
The smooth transition from Poisson to Wigner-Dyson statistics is quite generic for finite-size systems, e.g. [91][92][93][94][95], and is naturally explained in the Hamiltonian system (2.6)-(2.7) by the integrable structure no longer being exact at finite g 2 . At the same time, deformations which preserve integrability to leading order in the deformation parameter and only break it at higher orders are special and have been shown to be "weak" [52] in the sense that the onset of chaos is slower, i.e. the cross-over from Poisson to Wigner-Dyson occurs at larger couplings than for generic deformations. In [52], the authors considered the XXZ spin chain with an integrability-breaking current deformation which gives rise to long-range interactions. These current deformations were generated by the spin-chain boost method [53,54] which was originally introduced in the context of the integrable planar N = 4 long-range spin chains. This method guarantees that the deformation can be completed into a longrange integrable theory non-perturbatively in the deformation parameter. This similarity suggests that the two-loop planar dilatation correction is weak in the same sense, and we will explore this in the following.
We consider the transition from integrable to chaotic level statistics for the two-loop Hamiltonian (2.6) with parameters (2.7) as we vary g 2 from 0 to 1/4 and the same transition for a Hamiltonian where a twist parameter β ∈ R is introduced into the NNN-deformation as Note that here we are not twisting the leading NN term, and thus we cannot assume that this deformation corresponds to an integrable deformation for generic values of β even at O(g 2 ) and so it is not weak in the above sense.
For a range of values of β and g 2 ∈ [0, 0.25] we compute the level-spacing distributions and, assuming α = 1, find the values of ω that give the best fits of P B to the data. As can be seen from the bottom right plot in Fig. 1, the undeformed two-loop spin chain transitions from the integrable, ω = 0, case to the chaotic, ω = 1, distribution consistently slower than the generic twisted deformation which is compatible with it being a weak deformation. It is worth noting that although we are using the parameters (2.7), this Hamiltonian is for intents and purposes the XXX spin chain with a NNN interaction. It is somewhat surprising to find that this well-known chaotic system is weak in any sense 2 . However, there were some previous suggestions that this is indeed the case, e.g. the deformation was noted by [86] to be weaker than the interchain perturbation. More recently, [96,97] found that for the XXX chain with NNN deformation there are quasi-conserved charges . This is to be contrasted with the next-to-NNN (range-four) deformation for which there are no quasi-conserved charges at all.
One subtlety with our comparison of the SU(2)-invariant chain and the twisted chain is that the twisted deformation, in addition to breaking integrability, breaks the global SU(2) and parity symmetries. As a result it mixes larger subsectors (for L = 20 and M = 7 the dimension is m = 3876) than the undeformed theory and so the comparison is not between equivalent sectors of eigenstates. When g 2 is taken to 0 in the twisted case, the symmetries are restored and the subspaces with different charges decouple. The combined distribution of the corresponding sectors is neither Poisson nor Wigner-Dyson as it has many more spacings near 0. This explains the negative values of ω near g 2 = 0 for the twisted deformation. It is thus interesting to turn to the β-deformed gauge theory which is described by a twisted spin chain even at leading order and where the NNN deformation does not change the global symmetries.

Marginally deformed SU(2) spin chain
The β-deformed N = 4 SYM theory is a Leigh-Strassler deformation [55] preserving N = 1 supersymmetry, but breaking the SO(6) R-symmetry of N = 4 SYM to U(1) 3 . In terms of N = 1 superfields, Φ i with i = 0, 1, 2, the superpotential can be written as where in general the parameters g and β are allowed to be complex, though in this section we restrict them to be real. Additionally, requiring conformal invariance provides a constraint relating g , β, the gauge-coupling constant g and the rank N of the gauge group, as well as constraining the gauge group to SU(N ) rather than U(N ). The Lagrangian may be obtained from the undeformed Lagrangian by replacing all products of fields by a Moyal-like star-product corresponding to a noncommutativity in the internal SU(4) R-symmetry directions [77]. As a result of Filk's theorem [98], planar Feynman diagrams in the deformed theory are generally the same as those in N = 4 SYM theory up to an overall phase depending on the charges of the external fields 3 . Consequently, the perturbative dilatation operator is just a twisted version of the N = 4 spin chain and the all-loop spectrum of asymptotically long operators is given by twisted Bethe equations [60,100].
If we again focus on an SU(2) β sector comprising operators made of two complex scalars, X and Z, the resulting spin-chain Hamiltonian can again be written in terms of spin operators acting on each site, S i . The two-loop Hamiltonian was constructed in [76] by means of a unitary transformation which acts on the spin operators as and so we can find the SU(2) β Hamiltonian from the N = 4 spin chain by making the replacement The unitary transformation corresponds to a twist of the boundary conditions, so the spectrum remains integrable for finite β and λ and is given by the same Bethe equations but with an additional boundary factor. For later convenience we define a Hamiltonian with different twists in the NN and NNN contributions where the two-loop planar dilatation operator of the β-deformed theory corresponds to the parameters (2.7) and β 2 = β. Further, with this choice of parameters, there is a long-range integrable completion of this Hamiltonian and so the two-loop deformation, while not integrable at finite g, should be weakly chaotic in the sense of [52]. This can be seen in the numerical analysis of the level spacings in Fig. 2. The analysis is essentially the same as in the undeformed case with the simplification that, as the twist breaks the global SU(2) and parity symmetries, translation invariance is the only symmetry and thus we de-symmetrise by restricting to zero-momentum operators. When the coupling g 2 is set to 0, the Hamiltonian is just that of the twisted XXX spin chain and the distribution is Poisson, cf. Fig. 2 (top left). For sufficiently large values of g 2 the theory becomes chaotic and the distribution is well described by the GOE Wigner-Dyson distribution as can be seen in Fig. 2 (top right). That the distribution is orthogonal is again consistent with the discrete symmetries as the previous notion of time-reversal symmetry extends to this case even though the rotational invariance is broken: For a spin chain with a Hamiltonian composed of terms of the form the previous combinations of complex conjugation and either spin-flip or parity is a symmetry of the Hamiltonian. This implies that the statistics of the twisted spin chain with non-nearest-neighbour interactions are also GOE. When β 2 = β there is presumably no completion into a long-range integrable theory and so the deformation should be strongly chaotic. In Fig. 2 (bottom left) we compare different values of β 2 by computing the level-spacing distribution for a range of values g 2 ∈ [0, 0.25] and then finding the value of ω for which P B with α = 1 best fits the data. It can be seen that the transition from ω = 0 to ω = 1 is slowest for β 2 = β with the transition being quicker for increasing values of β − β 2 . It is not shown in the figure but the curve for β 2 = 0 is very close to that of β 2 = 2/5 and so the dependence appears to be on |β − β 2 | mod π. Additional curves for β 2 increasing beyond those of β 2 = 1/5 + π/2 move closer to the original curve, and are coincident when β 2 = β + π, and so the . Bottom right: Brody-distribution parameter ω as a function of g 2 for spin chains of varying length, L ∈ [16,22], with β = β 2 = 1/5. dependence on g 2 appears to reach a limiting curve as β 2 − β = π/2. This provides an interesting example of a deformation which smoothly interpolates from weakly to strongly chaotic.
The authors of [52] also demonstrated that the cross-over behaviour of weak deformations has a different finite-size scaling than seen in strongly chaotic systems. As the size of the system increases, the value of the cross-over coupling, g c , at which the system transitions from integrable to chaotic level statistics generally decreases. In [94], the authors found numerical evidence that g c scales as a power of the system size, with the exponent being universal though possibly depending on the RMT ensemble relevant to the chaotic regime and so on the discrete symmetries of the model. For interacting models the transition from Poisson to GOE Wigner-Dyson distributions has a crossover coupling which scales as L −3 . The same result was found in [52] for the XXZ spin chain with a NNN-deformation which is not SU(2) invariant. In contrast, for the weakly-chaotic current deformation the cross-over was found to scale as L −2 . The scaling of the cross-over coupling thus provides another test for weakly chaotic systems.
Different authors use different definitions of the cross-over, but the results should be independent of specific definitions. In this work we define g 2 c to be the coupling for which the level spacing is best fit by the Brody distribution with ω = 1/2 and α = 1 4 . We numerically determine the crossover coupling for a given L and M by computing ω for a range of couplings and then finding the parameters a and b for which this data is best fit by the function 1 2 (tanh(ax − b) + 1) with the cross-over coupling being given by b/a. The system size, in terms of the dimension of the Hilbert space of each zero-momentum subsector, is a function of both L and M and in principle we might hold either fixed and vary the other but this may not be the most physical comparison. We compute the cross-over coupling for different values of L but hold S z /L as close to constant as possible. This is similar to [94] which considered chains at half filling, S z = 0, and [52] which also considered fixed total spin S z = 1, though we include chains of odd length in our computations. In Fig. 3 we plot the estimated cross-over coupling for the Hamiltonian (2.26) of the weakly, β 2 = β, and strongly chaotic, β 2 = β + π/2, cases for chains with L ∈ [16,22] with S z /L 7/16. For the strongly chaotic case we find that the dependence on L is g 2 while for the weak case we find We have quoted errors for the coefficients but these should not be taken too seriously -they are simply the standard errors of the fit to the conjectured form of g c and are likely a significant underestimate. The parameters for the strongly chaotic case are consistent with previous results of b c = 3. For the weak case we find b c 1 which differs from b c = 2 found in [52]. This may be due to the different global symmetries of the two models. While neither model has a global SU(2) symmetry, the anisotropy of XXZ model considered in [52] is a much more significant breaking than the twist we have used. Nonetheless, these results are compatible with the weakly chaotic case having a scaling exponent smaller than in the strongly chaotic case.

Marginally deformed SU(3) spin chain
We now extend our considerations to larger sectors of the field theory, and specifically the SU (3) subsector of N = 4 SYM theory. It consists of single-trace operators made from holomorphic combinations of three complex scalars φ i , with i = 0, 1, 2, which can be combined with fermionic partners into the components of three N = 1 chiral superfields Φ i . As before, we consider the effect of marginal deformations on the spectrum of anomalous dimensions for such operators and study the spectral statistics in the context of chaos and integrability, however restricting to one-loop order in this section. The general Leigh-Strassler deformation is given in terms of two complex parameters, h and q, which we parametrise as q = exp(iβ + κ) 5 , with β, κ ∈ R, and its superpotential is is closed and one can extract the corresponding planar dilatation operator. We follow the notation of [56] and write it as an SU (3) spin-chain Hamiltonian with nearest-neighbour interactions. The i-th spin-chain site, where i = 1, . . . , L, can be in one of three states and the Hamiltonian can be written as where the operators E m,n i act on the i-th site as The authors of [56] considered values of h and q for which this spin chain is integrable and found R-matrices satisfying the Yang-Baxter equation for the cases all of which consist of Hopf twists of the undeformed theory. We are interested in the general case of arbitrary complex parameters h and q, and in understanding the level statistics of the corresponding spin-chain spectra. As in the SU(2) sector we begin with a study of the model's symmetries to properly de-symmetrise the spectrum by restricting to subsectors of states with the same global quantum numbers. We consider closed spin chains of fixed length, L, with periodic boundary conditions. As the Hamiltonian and boundary conditions are then translation invariant, the states are labelled by the total momentum P and, motivated by the connection to single-trace operators, we focus on spin chains with P = 0.
As can be seen from the form of the individual terms in the Hamiltonian (3.3), the Hilbert space splits into three sectors, each having a fixed value of M = L i=1 M i counted modulo 3. The Hamiltonian additionally has a manifest symmetry under the cyclic flavour transformation and for L = 0 mod 3 the three sectors are mapped into each other by this transformation. By direct diagonalisation for small L = 0 mod 3 (we considered L ≤ 13), we find that the Hilbert space generically further splits into two subsectors which can be seen in the eigenvectors having many zeroes. This is due to the states possessing a charge associated with a combined spin-chain parity transformation P and a discrete flavour transformation C. The flavour transformation C sends all both having the eigenvalue 3hh. In these cases it is possible to define the generalised parity transformation using any of the three flavour transformations v ↔ v + 1 v = 0, 1, 2 and this symmetry separates the spectrum into two, with the degenerate states in different subspaces. Nevertheless, this does not fully de-symmetrise the spectrum and as a result we do not consider L = 0 (mod 3) M = 0 sectors in this work. Finally, we also study the system with h = 0, which includes the integrable case q = e iβ with β ∈ R but also non-integrable cases where q is not a pure phase. In this case the Hamiltonian not only commutes with the total value of M , but also the number of each individual flavours is conserved and so we must work with states of fixed length L and fixed number of sites, R, with |1 , and, S, with |2 .

Level-spacing statistics
The computation of the level-spacing distribution for spectra in the SU (3) sector is essentially the same as in the previously discussed SU(2) and twisted SU(2) cases: We numerically compute the energy eigenvalues by direct diagonalisation, then de-symmetrise the spectrum by focussing on specific sectors of states, and finally we carry out the appropriate unfolding and binning of the data. We first consider the spin chain with L = 11 and M = 0 and take the deformation parameters to have the numerical values h = 1 3 exp i 11 13 , q = exp i 7 10 + 5 13 . (3.9) These do not correspond to any of the integrable points (3.5) found by [56]. The corresponding spin-chain system has a total number of 5369 states which can be organised into 2806 states with even and 2563 states with odd generalised parity. The distribution of level spacings is shown in the top left panel of Fig. 4. Using the generalised Brody distribution (2.16) and assuming ω = 1 for a Wigner-Dyson distribution, we find the best fit is for α = 1.98 for the positive-parity states and α = 2.12 for the negative-parity states. If instead we assume α = 2, we find the best fit is for ω = 1.00 and ω = 1.05, respectively. We can attempt to fit both parameters and in this case the results are less convincing, in particular we find {α, ω} is {1.71, 1.10} for positive parities and {1.32, 1.37} for negative, but this seems to be a feature of the sensitivity of the fitting process. Thus this distribution appears to be well described by the GUE Wigner-Dyson distribution with α = 2. Similar results are found for L = 11 with M = 1 and M = 2. The case where L is a multiple of three has different symmetries but, once complete de-symmetrisation is carried out, the results are again similar, cf. bottom left panel of Fig. 4 for the case of (L = 12, M = 1) for the three values of Z.
Thus, in contrast to the chaotic SU(2) spin-chain models studied above, here we find that the level-spacing distribution is well described by the unitary instead of the orthogonal ensemble. This is explained by the fact that generally there does not seem to be a time-reversal operator for the deformed SU(3) Hamiltonian. Although there is the cyclic symmetry transforming the flavours at each site as in (3.6), this cannot be combined with complex conjugation to form a symmetry. Additionally, in the case of non-vanishing complex h the combination of spin-chain parity and complex conjugation is not a symmetry of the Hamiltonian (3.3) as  while for qq = 1 the term E +1, +1 i E , i+1 breaks parity symmetry as it is not mapped into the corresponding E , i E +1, +1 i+1 term. The distribution of spacings is largely independent of the values of (q, h), however we can see that for particular cases which are known to be integrable the distribution becomes Poisson. For example if we take i.e. β = 0 and κ = 0, we find the distribution is again well described by the GUE Wigner-Dyson distribution as shown in Fig. 5. If we take β = 0 but keep κ = 0, we find that the spectrum is chaotic but now GOE. In this last case the orthogonal ensemble's reappearance is natural as the Hamiltonian for q ∈ R is real, and so we expect the ensemble to be GOE. Alternatively if we take κ = 0 and β ∈ R, we find the Poisson distribution corresponding to an integrable Hamiltonian. In Tab. 1 we show the best fit values for the parameters ω and α of the Brody distribution. In Sec. 4 we will be interested in studying the dynamics of spin chains in the limit of infinite length and low energy. This is a difficult regime to access by direct diagonalisation methods due to the rapidly increasing size of the Hamiltonian mixing matrices. In general, see e.g. [45,90], one expects the low-energy states to be less chaotic than those from the bulk of the spectrum. Nonetheless, we can analyse the level-spacing distribution for the lowest 6% of energy states, see the right panel of Fig. 5 and the starred data in Tab. 1. Even for these low-energy states, the distributions are quite close to the Wigner-Dyson distribution although the case of β = 0 is further from GOE. One might suspect that for sufficiently large L even the low-energy states are chaotic and this is consistent with what we will find in the LL model, though the cross-over from integrable to chaotic is non-trivial.

Spectral rigidity
The different RMT Gaussian ensembles not only describe the nearest-neighbour correlations of the unfolded spin-chain spectra, but also long-range spectral properties and even the statistics of eigenvectors. In this section we consider the spectral rigidity which is a measure for the regularity, or rigidity, of a spectrum resulting from level repulsion. It can be measured by the Dyson-Mehta statistic ∆ 3 (l) [101] defined as  interval starting points ε 0 taken from a discretisation of the unfolded spectral range, excluding a number of points at the boundaries. Increasing the interval length l increases the number of probed levels and thus ∆ 3 is a measure for long-range correlations. The solid lines in Fig. 6 show the expected behaviour of the Dyson-Mehta statistic for matrices in the Gaussian orthogonal and unitary ensemble, as well as for an uncorrelated spectrum, see [88] for further details and e.g. [49] for an efficient way to perform the extreme value problem. The Dyson-Mehta statistic is related to other statistical measures like the number variance and spectral form factor, and these relations can be made explicit by realising that all of these measures essentially probe the two-level correlations in a spectrum. For example the spectral form factor g(t), which has been studied e.g. in the context of quantum chaos in the SYK model [102], is related to the Dyson-Mehta statistic via [103] ∆ 3 (l) = 1 π ∞ 0 dt t 2 g(t)G(lt/2) (3.14) with In Fig. 6 we plot spectral-rigidity data for the L = 11 and L = 17 deformed SU(3) spin chains for various parameter configurations (q, h). At integrable points the Dyson-Mehta statistic follows the linear behaviour of uncorrelated spectra, reflecting the absence of level repulsion in these spectra. At non-integrable points it follows the logarithmic growth of the GOE and GUE prediction, respectively, according to the discrete symmetry properties of the spin-chain model. Together with the results for the short-range nearest-neighbour statistics from the previous section, these findings for the Dyson-Mehta statistics imply that fluctuations in spectra of the considered planar supersymmetric theories are consistently described by the Poisson or Wigner-Dyson RMT distribution according to their symmetry structure.

Transition from integrable to chaotic dynamics
We now examine the transition between integrable and chaotic statistics in the SU(3) spin-chain model. In particular, it is integrable in the case where h = 0 and κ = 0, and transitions to a chaotic model as either h or κ become non-vanishing. We can thus view h = 0 and κ = 0 as perturbations about an integrable line parametrised by β. In the case of h = 0, the perturbation breaks the global symmetry which separately preserves the number of |0 , |1 and |2 sites, and so we cannot smoothly turn on this deformation in the level-statistics analysis as the size of de-symmetrised subsectors changes. However for the case of h = 0 and β = 0, the subspace remains the same as we turn on κ and so we can study the transition from Poisson to Wigner-Dyson more clearly.
We study chains of length L ∈ [13,17] and consider sectors with the excitation densities R/L ∼ 1/4 and S/L ∼ 1/6 held as close to constant as possible as we vary L. We take β = 1/ √ 5 so that the chaotic distribution is GUE. For a range of values of κ ∈ [0, 0.5] we compute the value of ω for which the Brody distribution best fits the data, cf. left panel in Fig. 7. For larger values of the length, and so larger values of the dimension of the subspace, the value of ω increases more rapidly toward the Wigner-Dyson value of 1 with increasing κ. Additionally, for larger L, i.e. L = 16 and L = 17, it can be seen that ω is greater than zero even for very small κ and so there seems to be no region of integrable dynamics at finite κ. This is consistent with the hypothesis that at infinite L the dynamics are chaotic for any finite value of the deformation.
The curves of increasing ω can be fitted by the function 1 2 (tanh(aκ − b) + 1) which can then be used to determine the cross-over parameter κ c for which ω = 0.5. The value of this parameter for different values of L is shown in the right panel of Fig. 7 where it can be seen to scale as ∼ 1/L 3 . Fitting the cross-over values to κ c = N c /L bc , we find the best fit is for b c = 3.0 ± 0.4 , N c = (6 ± 7) × 10 2 (3. 16) which suggests that the deformation is strongly chaotic and which is consistent with the fact that there is no known integrable completion of this deformed spin chain. For later comparison with the Landau-Lifshitz model, it is also of interest to compute the crossover for the case of β = 0 where the chaotic dynamics are of GOE type. In this case, when κ = 0, the global SU(3) symmetries are restored and the limit is no longer smooth as the subsector splits up into smaller sectors. As a consequence, for small values of κ, the best fit for ω in the Brody distribution becomes negative. Nonetheless, we can examine the cross-over behaviour for a range of values κ ∈ [0, 0.5], where we keep only those values for which ω > 0. With this caveat, one can again see, cf. Fig. 8, that the onset of chaotic dynamics occurs earlier for larger L. The cross-over parameter is, as previously, computed by fitting it to a tanh-function and the scaling with L is shown in the right panel of Fig. 8. The scaling is still 1/L 3 with the best fit values being indicating that the deformation is again strongly chaotic. These results imply that for any finite κ, or any κ = κ 0 /L b with κ 0 constant and b < 3, the system will be chaotic as L → ∞. In Sec. 4 we will study the large L limit with b = 1 which gives rise to a Landau-Lifshitz model. In this limit we additionally focus on the low-energy states and one can ask whether the cross-over behaviour is different for these states. Although this is difficult to analyse in the spin chain, if we restrict to the low-energy states there is a clear increase in the

Landau-Lifshitz Models
In the context of the AdS/CFT correspondence, studying the long-wavelength limit about the ferromagnetic vacuum of the spin chain proved a useful approach in understanding the connection between the spectrum of planar anomalous dimensions and the string energies. The generalised Landau-Lifshitz model, which emerges as an effective two-dimensional action in the limit of long chains, L → ∞ with λ/L 2 fixed, was shown to agree with the string sigma model action expanded in the same limit [63,70]. Let us briefly review the derivation of the Landau-Lifshitz action from the continuum limit for SU(2) spin-1/2 chains, for details see e.g. [62]. The starting point is the definition of the coherent state for a single spin-1 2 particle as | n(ψ, φ) = e iφ cos ψ|1 + e −iφ sin ψ|2 , (4.1) where the vacuum, spin-up, state is denoted by |1 and the spin-down state is |2 (see App. A for more details on the construction of coherent states). For the full spin chain we can then define a coherent state as the tensor product over each lattice site We now consider the time-evolution operator for a spin chain with Hamiltonian H = L i=1 H i,i+1 and the transition amplitude Ψ a |e −iHT |Ψ b between two coherent states Ψ a,b |. We begin by breaking the time interval T into N steps of duration δ and look at the limit N → ∞ and δ → 0, with N δ = T kept fixed. At each intermediate time step t j we insert the resolution of the identity (A.8) so that where we have made the identification |Ψ a = | n(t 0 ) and |Ψ b = | n(t N +1 ) . Expanding in δ and using the identity which follows from the definition of the scalar product (A.7) combined with the assumption that the trajectories are sufficiently smooth, we can write the transition function as In the usual fashion this can be considered as defining a path integral over all trajectories of the vector n on the unit two-sphere with appropriate boundary conditions where we have introduced the action The first term is the topological Wess-Zumino term that follows from the definition of the coherent states, while the second encodes the specific dynamics of the system. For the spin chain with corresponding to the undeformed planar SU(2) one-loop dilatation operator where we re-introduced the overall coupling λ 16π 2 , this second term is given by where n i = (cos 2φ i sin 2ψ i , sin 2φ i sin 2ψ i , cos 2ψ i ) so that n 2 i = 1. If we now make the additional assumption that the states are such that n i varies smoothly as a function of the lattice position, then we can replace the discrete index i = 1, . . . , L with a continuous parameter σ ∈ [0, 1] in the limit of the large L. Expanding the lattice fields so that where the derivative with respect to σ is denoted by the primed fields, and taking the same continuum limit for the first term in (4.7), we find the Landau-Lifshitz action The fields (φ(t, σ), ψ(t, σ)), or equivalently n(t, σ), can be viewed as describing the long wavelength modes of the spin chain about the ferromagnetic vacuum which corresponds to setting n(t, σ) to a constant vector. Holding λ/L 2 fixed while taking L → ∞, we see that this limit corresponds to a semi-classical limit of the spin chain where classical solutions describe long spin-chain energy eigenstates. Restricting to closed periodic spin chains corresponds to imposing the boundary conditions φ(σ = 1, t) = φ(σ = 0, t) and ψ(σ = 1, t) = ψ(σ = 0, t) (4.12) on the coherent state. Furthermore, to restrict to zero-momentum spin chains, we impose the condition P = 0 with the generator of world-sheet translations given by P = dσ(p ψψ + p φφ ) = L dσ cos(2ψ)φ .

Higher-derivative LL model
At higher loop orders one can map the planar dilatation operator to a spin chain involving longerrange interactions. A description in terms of long-wavelength modes is still possible, but the effective action includes terms with higher numbers of derivatives [70]. For a spin chain with SU(2) symmetry and next-to-nearest-neighbour interactions, the effective action is fixed by symmetries to be of the form where we have re-scaled σ by 2π compared to the previous section. The Wess-Zumino term is fixed by kinematic considerations and so is the same for all spin-1/2 chains, while the coefficients a 0 , a 1 , a 2 depend on the detailed structure of the Hamiltonian. The derivation of the effective action is in general more complicated than simply taking the continuum limit as above, and one must carefully include the quantum corrections before taking the continuum limit [70]. For the next-to-nearest neighbour Hamiltonian (2.6) with the choice of parameters corresponding to the two-loop planar dilatation operator, the coefficients were computed in [70] and found to be 6 Again we only keep the leading term in the limit L → ∞ while holding λ/L 2 fixed. The equations of motion following from (4.14) can be written as where we have also expressed the contribution from the Wess-Zumino term using n. We are interested in understanding the dynamics of this model and in particular the presence (or absence) of classical dynamics analogous to the quantum chaos found in the spin-chain system. To this end we wish to explore different regions of the space of solutions and, as a complete categorisation of solutions for this non-linear field theory is not feasible, we make a simplifying ansatz φ(σ, t) = ωt , ψ(σ, t) = 1 2 q(σ) , (4.17) where ω is a constant. Note that in this section q corresponds to a phase-space coordinate and not the Leigh-Strassler parameter. The ansatz (4.17) automatically satisfies the zero-momentum constraint P = 0. It also includes the constant vacuum solution ψ = const. and the LL limit of the folded spinning string solution considered in e.g. [63,66] and thus we refer to it as the spinning string ansatz. The equations of motion imply that where q (n) is the n-th derivative of q.
If we first consider the case of the undeformed LL-model, where λ = 0, the equation is essentially that of an exact pendulum and can be explicitly solved in terms of elliptic integrals. As described in [63], there are essentially two classes of solutions depending on the initial conditions: those where the angle ψ oscillates over some finite interval and those where the angle increases without limit. When we include the higher-derivative corrections, the equations are non-linear and we must generally resort to numerical methods. To gain some insight we can set ω to zero, integrate (4.18) and introduce the variable u = q (1) which satisfies u = a 0 a 1 u + 2 1 + a 2 a 1 u 3 + c (4.19) whereu denotes the derivative w.r.t. σ and c is an integration constant. This is simply the equation for a particle moving in a double-well potential combined with a constant force term given by c.
For the choice of parameters corresponding to the two-loop dilatation operator (4.15), the potential well is inverted and so one again finds both unbounded and bounded motion. To this we must add the sin q-driving term which, depending on the sign of ω, will either be a restoring force or tend to drive the particle away from the q-origin. If the initial conditions are such that the particle is close to the u-origin, with small initial u-velocity and acceleration, and the sign of ω is such that the driving sin q-term provides a restoring force, the motion will be bound. Alternatively for large initial values or when the driving force is sufficiently large, the motion appears unbounded though numerical issues make following the trajectory difficult. For initial conditions on q, q (1) , q (2) and q (3) giving bound motion, we plot the time evolution of q(t) in both the undeformed and deformed theories, cf. left panel in Fig. 9. This illustrates that in the theory deformed by the two-loop higher-derivative terms, the motion is somewhat more irregular than the undeformed case.
To obtain a more quantitative characterisation, it is useful to switch to a Hamiltonian formulation of the reduced dynamics. Substituting the ansatz (4.17) into the action yields the Lagrangian L = −a 1 (q (1) ) 2 − a 0 (q (1) ) 2 − (a 1 + a 2 )(q (1) ) 4 + ω cos q . (4.20) Following the general procedure of Ostrogradsky (briefly reviewed in App. B), the phase-space variables are q, q (1) , p 0 and p 1 , and we define the Hamiltonian It is easy to check that the resulting Hamiltonian equations of motion are equivalent to those of the original model (4.16) upon substitution of the ansatz (4.17). Note that while the reduced theory is one-dimensional, due to the higher derivatives the phase space is four-dimensional and so the theory is not trivially integrable and there is a possibility of chaotic motion. One quantitative approach to chaotic dynamics is to compute the Lyapunov characteristic exponents (LCEs) [104][105][106] (see [107] for a review and further references) which measure the rate of change of the separation between initially nearby trajectories. For chaotic motion at least one LCE is positive, while for regular motion all LCEs are zero and so the maximal LCE (mLCE) is an indicator of the chaotic nature of orbits in a dynamical system. More specifically for an autonomous Hamiltonian system the spectrum of LCEs of order 1, {λ i } with i = 1, . . . , 2N , where 2N is the dimension of phase space, and ordered so that λ 1 ≥ λ 2 ≥ · · · ≥ λ 2N , consist of pairs having opposite signs, λ j = −λ 2N −j+1 for j = 1, . . . , N , with at least two LCEs vanishing. For each additional independent integral of motion, one more pair of LCEs will vanish. There are well-established methods for numerically evaluating the LCEs of a given dynamical system including an implementation in Mathematica [108] which we followed closely (see App. C for additional details of our numerical methods for computing LCEs). Taking q (1) (0) ∈ [0, 1] and ω = −1/10 we find that the mLCE is less than 0.01, and similarly fixing q (1) (0) = 0.1 and letting ω ∈ [−0.2, 0.03] we find that the mLCE is less than 0.08. Such values are relatively small and suggest the motion of the higher-derivative theory corresponding to the two-loop dilatation operator is not chaotic.
One might worry that this contradicts the finite-L level statistics, however the lack of chaotic motion should perhaps not be too surprising. In similar systems, e.g. the Duffing oscillator or double pendulum, the motion is regular at low energies and only becomes chaotic when the energy is sufficiently large. In the LL model following from the spin chain, we cannot increase the energy too much without the motion becoming unbound/unstable. This is also not unexpected from the spin-chain description as the LL model corresponds to takingg 2 = g 2 /L 2 fixed while L is large. The ratio of the next-to-nearest-neighbour to nearest-neighbour terms (2.9) becomes in this limit a value for which the spectral statistics are intermediate between Wigner-Dyson and Poisson. Relatedly, the S-matrix of the deformed LL theory computed about the q = 0 vacuum is known to factorise [109]. However it should be noted that this calculation was perturbative in g and q, and so is strictly valid only near g 2 = 0 and q = 0.
For comparison, we can study the dynamics of the deformed LL model for values of the parameters a 0 , a 1 and a 2 where, in the ω = 0 limit, the double-well potential in the u-variable results in bound motion for arbitrary initial conditions. For example, we can repeat the phase-space analysis with a 0 = a 1 = −1 and a 2 = 3/2 . (4.23) In this case the dynamics appear to be chaotic as can be seen in the mLCE which, for ω = −10 and initial conditions q(0) = 0 = p 0 (0) = p 1 (0), as well as q (1) = 1 10 , is found to be λ 1 = 0.17. In Fig. 9 we show the numerical evaluation of the two largest LCEs for times up to T = 400. There is one non-zero value and one vanishing LCE as is expected for an autonomous Hamiltonian system with a four-dimensional phase space. The two remaining LCEs are simply, to within numerical accuracy, the negatives of those shown in Fig. 9. Thus we see that the higher-derivative model is chaotic for generic parameters a 0 , a 1 , a 2 , but not for the specific values (4.23) corresponding to the two-loop dilatation operator which is compatible with the weak chaos seen in the level statistics.

SU(3) LL model
We now turn to consider the deformed one-loop SU(3) sector described in Sec. 3, where the Hamiltonian density is given in (3.3) and where the coherent state is cf. (A.14). For simplicity we consider only the case with h = 0 and q = exp(iβ + κ), and when taking the L → ∞ limit we keepβ = βL andκ = κL fixed (see [76] for a discussion in the SU(2) case). Repeating the procedure described above we find for the action with H[ n] = (∂ σ θ −κ 2 cos 2ψ sin 2θ) 2 + sin 2 θ (∂ σ ψ +κ 4 (1 + 3 cos 2θ sin 2ψ) 2 + cos 2 θ(∂ σ φ − cos 2ψ(∂ σ ϕ +β)) 2 + sin 2 2ψ(∂ σ ϕ +β 4 (1 + 3 cos 2θ)) 2 + 9 4 (β 2 +κ 2 ) cos 2 θ sin 4 θ sin 2 2ψ . The undeformed version of this action was studied in [64,65], while the deformed version was studied in [84]. If we introduce the variables ρ 1 = cos θ , ρ 2 = sin θ sin ψ , ρ 3 = sin θ cos ψ we can write the Hamiltonian as This is not quite the deformed action found in [84] as it is missing a term which is a result of the coherent state missing the vacuum configuration consisting of totally symmetrised products of all three spin-chain states. As a consequence of this, the Landau-Lifshitz action does not match the "fast-motion" limit of the world-sheet theory for strings in the deformed AdS geometry. This limitation can in principle be remedied by either considering generalised coherent states or by making an appropriate change of basis [84].
In order to study the dynamics of the model, we again restrict to a specific class of solutions which effectively reduce the model to one dimension. We start from the ansatz which implies that the continuum limit of the spin-chain momentum vanishes cos 2 θ∂ σ φ + cos 2ψ sin 2 θ∂ σ ϕ = 0 (4.31) and so corresponds to a zero-momentum spin-chain state. Starting with the undeformed theory, β =κ = 0, and substituting the ansatz into the equations of motion, one sees that the φ and ϕ equations are automatically satisfied while those of θ and ψ follow from the reduced Hamiltonian More generally, in the caseβ = 0,κ = 0 the equations of motion follow from the reduced Hamiltonian with the additional terms Hκ =κ − 2p θ sin 2θ cos 2ψ + p ψ sin 2ψ(3 cos 2θ + 1) + 9κ(α − 1) sin 4 θ cos 2 θ sin 2 2ψ , (4.33) where setting the coefficient α = 0 corresponds to the naive coherent-state limit of the spin chain and α = 1 corresponds to adding the contribution H sextic . In general, the φ and ϕ equations are not satisfied by the ansatz (4.30) whenβ = 0 so the equations of motion do not reduce to a one-dimensional Hamiltonian problem in the same way and we will not consider this case further. The reduced dynamics, in both the undeformed and deformed cases considered, are thus given by a Hamiltonian system with a four-dimensional phase space. We can compute the Lyapunov exponents for different parameter values and initial conditions.  Table 2: mLCE forκ =β = 0 integrated for T = 1600 with initial conditions p θ (0) = 3/10, p ψ (0) = 11/10, θ(0) = 1, ψ(0) = 1.
Undeformed Model Starting with the undeformed theory and taking, for example, the initial conditions p θ (0) = 3/10, p ψ (0) = 11/10, θ(0) = 1, ψ(0) = 1 (4.34) with w φ = w ϕ = 1 and numerically integrating for T = 1600, we find that the mLCE is less than 0.005. Similar values were found for different initial conditions and values of w φ and w ϕ see Tab. 2. That all such calculations can be done for such long times at relatively low numerical precision (only 30 digits) is another signature of the non-chaotic nature of the system. Using the iterated method to calculate the full LCE spectrum, we similarly find that all exponents rapidly converge to zero, cf. Fig. 10.
Deformed Model Computations in the deformed theory are more difficult as, due to the chaotic nature of the system, greater numerical precision is required to compute an accurate solution. For example, with the same reference initial conditions (4.34), we can integrate the equations of motion withκ = w φ = w ϕ = α = 1 to T = 100 to an accuracy of 10 −9 at working precision 55. That is, integrating forward and back gives solutions which differ by no more than 10 −9 in any phasespace coordinate which are generally of order 1. However, by T = 110 the accuracy has decreased to 10 −3 and by T = 150 solution is completely unreliable. In order to get a reliable solution for this time interval it is necessary to increase the working precision to 75. Nevertheless, even after only T = 100 we have reasonably good results for the mLCE with the value 0.59 compared to the T = 150 value of 0.58. There is an inherent error in this method of computing the mLCE in that a random direction in phase space is used, but the use of finite-time intervals is the largest source of error in our estimates. In general the results seem to be reliable to within at least ±0.05 after T = 100. Using the iterated method to compute the full LCE spectrum, we can study the convergence over time. In Fig. 11 we plot the largest LCE for each of the considered configurations, two of the remaining LCEs are zero and the fourth is just the negative of the first to within numerical accuracy. It should be noted that in the iterated calculation the solution is reliable over each iteration, but generally not over the whole time interval. This potentially results in a certain averaging of the LCEs over different phase-space points in an uncontrolled fashion. Nonetheless, the results for the LCEs appear to be robust, the iterated method gives 0.54 for the mLCE with the same initial conditions and parameter values, (4.34),κ = w φ = w ϕ = α = 1. This is possibly due to the LCEs being constant over at least small regions of phase space.

Transition from integrable to chaotic dynamics
In the final part of this section, we study the transition from integrable to chaotic dynamics in the deformed LL model. Fig. 12 shows the variation of the mLCE with changingκ. One notable feature is that for small, but finite values ofκ the mLCE appears to vanish. For example, with the initial conditions (4.34) andκ = 1/4, w φ = w ϕ = α = 1, the mLCE is less than 5 × 10 −2 when we integrate out to T = 100, and 5 × 10 −3 when we integrate out to T = 1600. This suggests it converges to zero, see Tab. 3. By comparison, forκ = 9/32, we find the mLCE is 0.08 at T = 100 and 0.1 at T = 400, suggesting it will converge to a small but non-vanishing number. This implies that the dynamics at this phase space point only become chaotic when the deformation parameter crosses some critical value close toκ = 1/4. Further increasing the value ofκ generally, though not monotonically, results in an increasing value of the mLCE, cf. Fig. 12.
Note that different regions of the phase space generally give different Lyapunov exponents for fixed values ofκ. For example, the initial conditions p θ (0) = 3/10, p ψ (0) = 11/10, θ(0) = 1, ψ(0) = 1/5 are not chaotic forκ = 1/2 with the mLCE being less that 7 × 10 −3 , however they  become chaotic by κ = 3/4, see left panel in Fig. 12. Thus it seems that different points in phase space have different critical values, but for a sufficiently large value ofκ it seems reasonable to suppose that all regions of phase space are chaotic. We also studied the behaviour of the model for different values of the parameters w φ , w ϕ and α, see the right panel in Fig. 12. While the values of the LCEs differ, generally being lower than for w φ = w ϕ = α = 1, the behaviour is generically the same. In particular, the dependence onκ is the same. There does not seem to be a qualitative difference when we remove the additional sextic term and so chaotic behaviour for non-zero values ofκ is found in both the naive coherent-state limit and when the potential term needed to match the string model is included.

Conclusions
In this work, we have considered the statistical properties of anomalous dimensions in planar N = 4 SYM theory and its marginal deformations. We have found that the spectral statistics for the one-loop marginally-deformed theory, in an SU(3) sector with generic parameters, h = 0, β = 0, κ = 0, is that of GUE RMT. This implies that the model is chaotic according to the usual definition following from the BGS conjecture. Additionally, for the case of the imaginary-β theory, where h = β = 0 and κ = 0, the ensemble becomes GOE with the statistics still being those of a chaotic system. This provides a weak-coupling analogue of the classical chaos seen in the holographic string dual [39].
To further characterise the dynamics, we study the distribution of level spacings as the parameters transition from values where the system is known to be integrable to those where it is chaotic. For infinite-dimensional quantum many-body systems, it is expected that the transition to chaotic dynamics occurs for any non-zero value of the deformations, while for finite systems there is a smooth cross-over. We consider the scaling of the cross-over coupling, defined in terms of the Brody distribution, as a function of the spin-chain length L. For the case of the imaginary-β theory, the scaling is 1/L 3 which is consistent with generic quantum chaos, e.g. [94]. For comparison, we study the transition behaviour in N = 4 SYM theory and its twisted analogue, where the integrable one-loop nearest-neighbour XXX spin chain is deformed by a two-loop NNN-interaction. In this case the statistics at finite coupling are those of the GOE Wigner-Dyson distribution. The cross-over coupling at fixed L is higher than for generic deformations and, in the twisted case, scales as 1/L. We thus find that the two-loop theory is weakly chaotic in the sense of [52]. This is consistent with the fact that the two-loop deformation, while it breaks integrability at finite coupling, is perturbatively integrable and can be completed into an integrable long-range theory, namely quantum strings in AdS 5 ×S 5 . Viewed from a purely spin-chain perspective, this result is quite interesting as it suggests that the SU(2)-invariant NNN-deformation of the XXX model, and the twisted versions, are only weakly chaotic. This is a perhaps surprising result, however it is consistent with the realisation that the NNN deformation is weaker than an interchain perturbation [86], and that for the NNN deformation there are quasi-conserved quantities while for the next-to-NNN deformation there are no quasi-conserved charges at all [96,97].
Finally, in order to further our understanding of the relationship between the spin-chain statistics and string-theory dynamics, we studied the behaviour of generalised Landau-Lifshitz models. These models are known to capture the dynamics of both string solutions with large angular momenta and long-wavelength states of the spin chain. In the integrable case, they give a very explicit mapping between string solutions and spin-chain states, and thus provide a potential approach to understanding how the chaotic dynamics interpolate between weak and strong coupling. For the higher-derivative LL model corresponding to the two-loop dilatation operator, we find that the classical dynamics has vanishing Lyapunov exponents. This is consistent with the classical integrability of the string theory, but non-obvious from the spin-chain perspective. It is potentially explained by the weakly chaotic nature of the two-loop dilatation operator combined with the fact that the LL model captures only the low-energy states which are generally the least chaotic. For the imaginary-β theory, we find non-vanishing Lyapunov exponents consistent with classically chaotic dynamics. This provides further qualitative evidence for a matching between the string and spin-chain chaos. One notable feature of the imaginary-β LL dynamics is that the transition to non-vanishing Lyapunov exponents occurs at finite values of the deformation. This is familiar in classically non-integrable systems, but is different from that found in the spin-chain description. It suggests that there may be interesting transition dynamics in spin chains if one combines the thermodynamic limit with a low-energy limit.
A key assumption is that we can invert the expressions for p i,m i −1 to writeq i,m i −1 = x (m i ) i in terms of p i,m i −1 and q i,r i with r i < m i for each i = 1, . . . , n. However we do not express the remaining variables q i,r i with r i ≤ m i − 1 in terms of the momenta. We now define the Hamiltonian

C Numerical computation of LCEs
We briefly summarise our approach to the numerical computation of LCEs which follows closely the methods described in [108,107]. Given the n-dimensional phase-space variables x, to compute the LCEs we numerically solve the equations of motion with initial conditions x(t = 0) = x 0 , simultaneously with the variational equation where J is the matrix Jacobian and we take Φ(t = 0) to be the n × n identity matrix. The matrix Φ(t) provides a linear map on the space of tangent vectors which can be used to evolve the deviation vectors. That is, given some initial configuration x 0 which corresponds to a flow with coordinates F t ( x 0 ) at time t, we can consider a system with infinitesimally shifted initial conditions x 0 + w(0). The deviation vector at time t is then given by In practice, one numerically solves (C.1) and (C.2) for a large time T . One issue is that for chaotic systems it can be difficult to maintain sufficient numerical precision when integrating over large times as, almost by definition, small errors rapidly grow in the flow trajectory. To ensure that our numerical integration is sufficiently precise, we first integrate for time T , then take the end point as our new initial condition and integrate backward in time to check that we find the original trajectory. For example, if we consider the non-chaotic case of the higher-derivative LL model (4.21) with spinchain parameters (4.15) and ω = − 1 10 , we can consistently numerically integrate the equations of motion for T = 400 using the Mathematica function NDSolve with the "ExplicitRungeKutta" method with working precision of only 14 digits. By comparison, for the chaotic case where the parameters are as in (4.23) and ω = −10, we need to have a working precision of 60 digits so that the forward and backward integrations are everywhere consistent to ±0.003.
A second practical issue is that the vectors w(t) grow very quickly. To address this, one can split the integration range into k steps of length τ = T /k λ 1 ( x 0 )) = lim k→∞ 1 kτ k r=1 ln α r , (C. 6) with expansion coefficient α r = w(rτ ) w((r−1)τ ) . We can use the linearity of (C.4) to rescale w((r − 1)τ ) to have unit length at each step, and compute the expansion coefficient for the evolution of the unit vector over the interval [(r − 1)τ, rτ ].
More generally one can compute the full spectrum of LCEs. To do this one starts with the definition of the LCE of order p, which is given in terms of the volume spanned by p linearly independent vectors w 1 , . . . , w p whose evolution is given by the variational equation as in (C.4) w 1 (t), . . . , w p (t)) vol (p) ( w 1 (0), . . . , w p (0)) . (C.7) Further, one uses the known relation between the LCE of order p for a suitably chosen subspace and the p largest LCEs of order 1 In our numerical calculations, we start from a set of n orthonormal vectors w 1 (0), . . . , w n (0) and we again split the integration into k steps of length τ . After numerically integrating (C.1) and (C.2) we compute the evolved vectors w 1 (τ ), . . . , w n (τ ) and use the Gram-Schmidt procedure to compute an orthogonal, but not normalised, basis˜ w 1 (τ ), . . . ,˜ w n (τ ). The evolved volume is given by the product of norms vol (n) ( w 1 (τ ), . . . , w p (τ )) = ˜ w 1 (τ ) . . . ˜ w n (τ ) . (C.9) As above in the case of the mLCE we can rescale the n vectors so they are again orthonormal and then take these vectors to provide our initial conditions. Iterating the calculation we have that which can be approximated by computing for some large T = kτ until convergence is found.

C.1 Lorenz Model
A very well-studied system, originating in atmospheric science, is the Lorenz model described by the variables (x, y, z) with the equations of motioṅ x = σ(y − x) ,ẏ = x(R − z) − y ,ż = xy − bz .
(C. 13) We use this relatively simple system to check our numerical methods. For σ = 10 and b = 8/3 the model is known to be chaotic whenever the parameter R, the Rayleigh number, is greater than a critical value R 24.74. The trajectories of the system generally have one positive Lyapunov exponent and it follows from Liouville's formula that the sum of the Lyapunov exponents is −(σ + Step LCEs σ =10, b= 8 3 , R=28 Figure 13: Plot of the two largest LCEs for the Lorenz model. The initial conditions are x(0) = 29, y(0) = 2, z(0) = 20. The total evolution time is T = 400 which is broken into 5000 iterative steps.