One-Loop Non-Planar Anomalous Dimensions in Super Yang-Mills Theory

We consider non-planar one-loop anomalous dimensions in maximally supersymmetric Yang-Mills theory and its marginally deformed analogues. Using the basis of Bethe states, we compute matrix elements of the dilatation operator and find compact expressions in terms of off-shell scalar products and hexagon-like functions. We then use non-degenerate quantum-mechanical perturbation theory to compute the leading $1/N^2$ corrections to operator dimensions and as an example compute the large $R$-charge limit for two-excitation states through subleading order in the $R$-charge. Finally, we numerically study the distribution of level spacings for these theories and show that they transition from the Poisson distribution for integrable systems at infinite $N$ to the GOE Wigner-Dyson distribution for quantum chaotic systems at finite $N$.


Introduction
The eigenvalue problem for the dilatation operator, D, acting on the set of gauge-invariant local operators, O i , in N = 4 super Yang-Mills (SYM) theory, has been of continued interest due to its role as a proving ground for novel calculational techniques and because of its importance in the AdS/CFT correspondence. The operator dimensions, ∆ i = ∆ i (g YM , N ), are non-trivial functions of the coupling g YM and N , the rank of the gauge group.
In perturbation theory we can expand the dilatation operator in powers of the 't Hooft coupling λ = g 2 YM N D = k g 2k D 2k with g 2 = λ 16π 2 (1.2) and at each order in g 2 we can further consider the large-N expansion of the operator dimensions.
A key development [1] was the insight that for the so(6) sector of operators the one-loop, O(g 2 ), leading large-N anomalous dimensions can be computed by means of an integrable spin chain. Single-trace operators composed of L scalar fields were identified with closed spin chains of length L and the planar dilatation operator with a spin-chain Hamiltonian that can be diagonalised by use of the Bethe ansatz. This was subsequently extended to the full one-loop theory [2] and to higher orders in perturbation theory [3] as well as being observed at strong coupling [4]. This prompted a great deal of work and ultimately lead to non-perturbative results for the spectrum of planar anomalous dimensions first through the thermodynamic Bethe ansatz and subsequently by means of the Quantum Spectral Curve (QSC), see [5] for reviews. While less is known about non-planar anomalous dimensions there are a number of impressive perturbative results for specific operators. For example, twist-two operators at four loops were studied using standard Feynman diagramatics [6] as well as twistor methods [7], and the fourloop non-planar correction to the cusp anomalous dimension was computed in [8]. Additionally, the Hexagon formalism for correlation functions, [9][10][11], provides an integrability-based method to studying non-planar N = 4 SYM. The framework can be used to evaluate higher-point functions at higher genus and one can then in principle extract non-planar dimensions from OPE limits. A second approach, which was applied at tree-level in [12], is to directly compute the non-planar twopoint function by considering the four-point function with two operators taken to be the identity. Moreover, there are results for the complete dilatation operator including its non-planar parts which was found at one-loop for the so(6) sector in [13][14][15] and for the full theory in [16]. Also at one-loop, a spin-bit model that captures some features of the string-bit formalism for interacting strings [17] was used in [18] to compute non-planar corrections for operators in the su (2) and sl (2) sectors. The full two-loop dilatation operator in the su(2) sector was found in [3] and extended to the non-compact su(1, 1|2) sector in [19] by using the superconformal symmetry of the theory.
The problem of diagonalising the non-planar dilatation operator has itself been studied using a number of techniques. One approach makes use of group-theoretic insights, initially developed in the context of 1/2-BPS operators [20], to make an appropriate choice of basis operators diagonalising two-point functions. Alternatively, one may attempt to perturbatively compute 1/N corrections about the planar limit using quantum-mechanical perturbation theory. One splits the dilatation operator into a leading planar part and subleading off-diagonal terms, which mix single-and multitrace operators, and then computes matrix elements of the subleading terms in the basis of planar eigenstates. Such an approach was used to compute the large R-charge limit of non-planar dimensions of two-impurity BMN operators in the su(2) sector at one-loop and at two-loop level [13][14][15]. In this work we will generalise this approach for the one-loop dilatation operator in the su(2) sector by making use of Bethe states describing arbitrary numbers of excitations or magnons. In Sec. 2 we compute the action of the off-diagonal terms on Bethe states, schematically shown in Fig. 1(a), and write the overlap of the resulting state with the tensor product of two Bethe states, corresponding to a double-trace operator, in terms of off-shell scalar products. The latter can be computed efficiently using the algebraic Bethe ansatz [21,22] or, equivalently at this order, a Hexagon-like formalism [9]. We similarly find the action on tensor products of Bethe states, see Fig. 1(b), and compute the overlap with single-trace Bethe states in terms of off-shell scalar products. These overlaps can in principle be used to find corrected anomalous dimensions, or equivalently, spin-chain energies by performing the sum over all such overlaps. Naively, as the overlaps are of order 1/N , we expect the leading correction to be of order 1/N 2 . However, due to degeneracies in the planar spectrum, non-degenerate perturbation can fail and can result in order 1/N corrections to dimensions. The issue of planar degeneracies between operators with different numbers of traces was noted in [23,15,3], see also [24], and requires solving a non-trivial mixing problem. We will instead consider deformations of N = 4 SYM theory for which these degeneracies are lifted. In particular, we consider β-deformed N = 4 SYM which is a marginal deformation of the maximally supersymmetric theory [25] preserving N = 1 supersymmetry. We concentrate on the case β ∈ R for which the theory is exactly conformal to all loop orders in the planar limit [26]. Most importantly, the planar spectral problem for the β-deformed theory is described by an integrable twist of the undeformed spin chain. Twisted asymptotic Bethe equations at one and higher loops were found in [27], while the twisted QSC was found in [28] and used in [29] to study the anomalous dimensions of operators corresponding to the Konishi multiplet. In Sec. 3 we fix the form of the deformed non-planar dilatation operator, which includes additional double-trace terms, and then expand the action on operators in powers of 1/N . In addition to the leading planar piece and the subleading 1/N terms, which are similar to those occuring in the undeformed theory, there are 1/N 2 and 1/N 3 terms which contribute to non-planar dimensions. We compute the matrix elements of the subleading dilatation operator in the basis of Bethe states and study the BMN limit of large R-charge, J = L − 2 → ∞, where perturbation theory can be rewritten in terms of the effective loop-and genus-counting parameters [30] g = g 2 YM N 16π 2 J 2 and g 2 = The anomalous dimensions of two-impurity BMN operators, which additionally depend on an integer parameter n and rescaled deformation parameter b = βL/π, can be written in terms of rescaled energiesẼ (k) = J 2−2k E (k) ∆ n (g , g 2 , b, J) = L + g Ẽ (0) n (b, J) + g 2

2Ẽ
(2) n (b, J) + O(g 4 2 ) + O((g ) 2 ) (1.5) and we computeẼ (2) n (b, J) through O(J −1 ). Besides the lifting of degeneracies, there are a number of further advantages to considering the deformed theory. At a technical level it allows us to consider singular solutions to the undeformed Bethe equations. Such solutions correspond to finite energies but have singular wavefunctions and so matrix elements in the undeformed theory are not well defined. Instead they can be computed in the deformed theory, where the deformation parameter regularises singularities in the wavefunctions [31], and the limit of vanishing deformation parameter can be smoothly taken. More conceptually, the dependence of the spectrum on a continuous parameter allows for the phenomenon of level crossing whereby eigenvalues become degenerate at special values of the parameter. As was noted in the early days of quantum mechanics by von Neumann and Wigner [32], given a generic theory depending on a number of parameters it is necessary to tune at least two parameters to cause energy levels to cross and produce a degeneracy. Subsequently it was shown by Teller [33] that the surfaces E = E(β 1 , β 2 ) representing energy levels depending on two such parameters β 1 , β 2 are connected at points like the two sheets of a degenerate cone. In [34] the dimensions of local operators in N = 4 SYM were analysed as functions of the 't Hooft coupling and it was shown, by re-summing the large-N expansion, that when N is held fixed anomalous dimensions do not cross as λ is varied. Studying the one-loop spectrum as a function of the deformation parameter β, we similarly find that at finite N the anomalous dimensions repel and it is only at large N that they cross.
Level repulsion is characteristic of a chaotic system where energy levels are correlated and so avoid each other, while in an integrable system they are uncorrelated and move independently, crossing on occasion. The phenomenon of (non-)repulsing energy levels can be studied by looking at the distribution of spacings between neighbouring energy levels. If one computes the probability, P (s)ds, that the normalised spacing between adjacent levels lies in the interval between s and s+ds, one finds that for a generic, chaotic, quantum system P (s) → 0 as s → 0. For integrable systems it is generally the case that P (s) goes to a constant as s → 0 which reflects the presence of hidden symmetries in these models. In Sec. 4 we numerically study the spectrum of both the deformed and undeformed theories and show that in the planar limit the spectral distribution is Poisson, consistent with integrability, while at finite N the distribution is Wigner-Dyson and corresponds to that of the Gaussian Orthogonal Ensemble (GOE) random matrix theory. We are thus able to numerically study the transition from quantum-integrable to quantum-chaotic systems in the context of interacting four-dimensional gauge theory as we vary N . An analysis of the transition from Poisson to Wigner-Dyson statistics in a context similar to planar N = 4 SYM was performed in [35] which considered and integrability breaking deformation of the XXZ spin-chain. The appearance of quantum chaos in the spectrum is in fact quite natural if we view the dilatation operator as the Hamiltonian of the theory defined on R × S 3 and multi-trace operators as defining states somewhat analogous to large nuclei in QCD. Indeed it was the work of Wigner [36] and Porter & Rosenzweig [37] on the statistical properties of the energy levels of highly-excited nuclei that lead to much of the initial interest of physicists in random matrix theory [38].

Non-Planar Dilatation Operator
In order to fix our conventions and notations we briefly review the one-loop dilatation operator of N = 4 SYM. We follow closely the work [3] where more details regarding the calculations and generalisations to higher loops can be found. N = 4 SYM theory contains six scalar fields, (φ I ) a b , I = 1, . . . 6, a, b = 1, . . . N which transform as a vector of the SO(6) SU (4) R-symmetry and in the adjoint representation of the SU (N ) gauge group. We will restrict ourselves to the su(2) sector comprising operators made of products of traces of two complex scalar fields, Z = 1 √ 2 (φ 5 + iφ 6 ) and X = 1 √ 2 (φ 1 + iφ 2 ), and so we consider operators such as Tr(Z 1 ) , Tr(XZ 1 XZ 2 ) , Tr(XZ 1 XZ 2 Z 3 )Tr(XZ 4 )Tr(XX) . This sector is known to be closed under the action of the dilatation operator and does not mix with operators containing other scalars, field strengths or fermions. To describe the action of the dilatation operator it is useful to make use of the notation for functional derivatives of fields, for example where it is assumed that A and B do not contain any Z's. The N −1 terms are due to the fact that we are considering the SU (N ) gauge theory. This is not particularly important for N = 4 SYM and we could equally well consider a U (N ) gauge group, however it will become relevant when we subsequently consider the β-deformed theory. Using this notation, the tree-level dilatation operator in the su(2) sector can be written as D 0 = Tr(ZŽ) + Tr(XX) (1.10) and simply counts the number of fields present in a given operator. The one-loop correction to the dilatation operator is then given by [13] where the normal-ordering markers : : indicate that the functional derivatives do not act on the fields in D 2 itself. We can find the action on multi-trace operators by repeated use of identities (1.9). For example on the length-six single-trace operator, Tr(X 2 Z 4 ), we have where we see that the leading term in a large-N expansion corresponds to a superposition of singletrace operators and the subleading term is a double-trace contribution. In general, we can decompose the action of the one-loop dilatation operator on multi-trace operators into planar and non-planar pieces The planar piece, H (0) , leaves the number of traces in an operator unchanged, while the non-planar corrections, H ± , which are suppressed by a factor of 1/N , increase/reduce the number of traces in a given operator. In order to find the eigenvalues of D 2 , one can first solve the planar problem using integrability and then attempt to use perturbation theory to find the 1/N k corrections.

Planar Theory and Integrability
Mapping the problem of computing anomalous dimensions to that of computing integrable spinchain energies proved to be an important step in solving the planar spectral problem. We will make use of the spin-chain notation and the results from integrability to organise the computation of non-planar corrections. We thus review the coordinate Bethe-ansatz approach to integrable spin chains here. Single-trace operators with M insertions of X fields in a background of (L − M ) Z's will be denoted as Tr( (1.14) Multi-trace operators with K traces and M insertions of X fields, where M = K k=1 M k , can be denoted by products of such states which is an element of the symmetrized tensor product. It will often be convenient to use the compressed notation |{n} and k |{n (k) } .
We can write general single-trace states as linear combinations of this basis in terms of a wavefunction ψ {n} which can depend on the quantum numbers describing the particular state. For example we will consider states with M impurities characterised by the momenta {p} = {p 1 , p 2 , . . . , p M } of M excitations, or magnons (1.16) The sum is over the positions of excitations ranging over the nested values 1 ≤ n 1 < n 2 < · · · < n M ≤ L. We will compute overlaps of such spin-chain states and so we define the natural dual basis 1 so that the scalar product for states |ψ and |ψ is given by The action of the planar dilatation operator can be defined in this basis and is given by the well-known formula where S(p j , p k ) is the two-magnon S-matrix In this definition we have made a particular choice for the overall, non-physical, phase of the wavefunction which is convenient for our subsequent purposes. For these states to satisfy periodic boundary conditions the momenta must satisfy the Bethe equations, i.e. for each j = 1, . . . , M Each eigenstate corresponds to a solution of the algebraic equations (1.23) and the energy eigenvalue is given as a sum over individual magnon energies The cyclicity of the trace for gauge-theory operators becomes the condition that the spin chain is invariant under the shift n j → n j + 1 and so we consider only states which satisfy the condition M j=1 e ip j = 1 . (1.26) It is useful to introduce rapidity variables u j = 1 2 cot for each excitation, which we can use to rewrite the S-matrix and the individual magnon energies as (1.28) It will be useful to define the quantity so that the S-matrix is given as Additionally, it is convenient to define the normalisation factors for states involving momenta p 1 , p 2 , . . . corresponding to rapidities u 1 , u 2 , . . . as and generalisations such as N (p, q) = N (p)N (q). Finally, similar to [39], we will use the following short-hand notation for products

Perturbative Non-Planar Anomalous Dimensions
In this section we study the action of the non-planar dilatation operator on the planar eigenstates, the Bethe vectors, and subsequent overlaps with other Bethe states. We consider first the action of H − on a double-trace operator corresponding to the product of a length L q Bethe state |{q} with Q excitations and a length L r state |{r} with R excitations. There are L q L r terms corresponding to the action of the dilatation operator on each pair of sites of the two spin chains. However, the terms where it acts on a Z field at the i-th site of the spin chain can be rewritten using the Bethe equations, so that they become equivalent to the action on a Z at the first site. This can be seen by gathering all the terms in the state (1.16) with a Z field at the i-th site and using the cyclicity and on-shell condition to write them with the Z field at the first position: Analogously, the action on an X field at the i-th site can be rewritten as the action on the same chain with the X field placed at the first site where the terms on the right hand side all correspond to single-trace operators. The overlap with a dual state {p}|, of length L p = L q + L r and with P = Q + R excitations, can then be computed where we define the sets {t λ } b a = {t λ(a) , . . . , t λ(b) } and denote products of exponentials over such sets using the notation e iL(t λ ) b a = b i=a e iLt λ(i) . The factors of P L (z) in (2.4) correspond to the geometric sums of exponentials in the wavefunctions which can be rewritten as sums over ordered partitions (2.5) Using these notations we can write a similar expression for overlaps of H + as sums over ordered partitions: The normalisation in this case is defined slightly differently with N + = N /S, where S is a symmetry factor that equals 2 when the states in the double trace are equal and 1 otherwise. Carrying out the geometric sums via (2.5) makes these formulae useful for analysing states of arbitrary lengths. However, while these expressions are reasonably compact, they involve sums over permutations for each of the sets of external momenta and so they rapidly become impractical as the number of excitations grows. The same growth is known from the computation of spin-chain scalar products in the coordinate Bethe ansatz and by making use of known results in this case we can find further simplifications.

Matrix Elements from Spin-Chain Scalar Products
The scalar product of two Bethe states involves double-sums over permutations and so is generally complicated to evaluate. Fortunately, there are well-known formulae for such scalar products which were developed in the Algebraic Bethe Ansatz (ABA) approach to integrable spin chains (see Sec. A for a brief review). In the case where both sets of momenta {k} and {l} do not satisfy the Bethe equations, i.e. they are off-shell, the scalar product can be written as a sum over partitions of the sets of momenta into subsets of equal cardinality [40], see (A.14). Similar simplifications can be used to rewrite the expressions of the overlaps (2.4) and (2.6). Each term in the formulae for the overlaps non-trivially involves one momentum of an excitation from the single-trace operator, which we label p j , and one excitation momentum from the double-trace operator, i.e. from either {q} or {r}, which we label as q i or r i . The remaining momenta are simply contracted using a rescaled spin-chain scalar product We can thus write the overlaps (2.4) and (2.6) in terms of the off-shell scalar products by splitting the single-trace excitation momenta into three subsets, {p} = s ∪ t ∪ {p j }, with the cardinality of s equal to that of {q}î (or {r}î) and the cardinality of t equal to that of {r} (resp. {q}). In terms of off-shell scalar products, the overlap of H − can then be written as a sum over all such splittings and that of H + as In addition to the scalar products of Bethe states these expressions involve (e ip − 1) factors, which are essentially the same as arise in the planar dilatation operator, and ordering factors for which we introduced the notations (2.11) These terms account for the phase acquired by the p j magnon as it is shifted around the chain before being contracted with a magnon on the double-trace operator. For each configuration there are two different ways to carry out this reordering and the overlap is a linear combination of both. Spin-chain scalar products have previously appeared in the context of N = 4 SYM in the computation of structure constants. In the all-order hexagon approach, [9], structure constants are written as sums over partitions of the magnon excitations and it was noted that this formulation is related to the scalar-product formula of Korepin [40]. It is therefore convenient to use a tree-level version of the hexagon formulation of scalar products , (2.13) in order to rewrite the overlaps (2.9,2.10).

A Hexagon-like Formulation
In the previous section we obtained the non-planar dilatation operator overlaps as sums over partitions of the rapidities, in a way that is reminiscent of the Hexagon formulation of three-point correlation functions [9]. In that context, the partitions of the rapidities arise naturally in the largevolume regime where the correlation function is broken down to its simplest building blocks, the Hexagon form factors. Crucially, these form factors satisfy a set of axioms which, together with the diagonal symmetries and some educated guesswork, can be used to obtain an all-loop description of structure constants. A particular feature of the Hexagon is its conical defect which is associated with the existence of three asymptotic regions and corresponds to a monodromy composed by three crossing operations. In the context of non-planar overlaps between a single-trace and a double-trace operator, a similar role seems to be played by the three distinct traces. In this section we investigate the properties of the objects arising from the action of H + and H − and find that they satisfy the same form factor axioms as appear in the context of correlation functions. The sum over determinants occuring in our rewriting of off-shell scalar products (2.12) can be found in a straighforward way from the hexagonalization of three-point functions [9]. To be precise we consider the three-point function of two unprotected operators in the SU (2) sector, one with X excitations and the other withX, and one rotated half-BPS operator. The X andX fields must be Wick contracted at tree level in order to produce a non-vanishing contribution. If there are l Wick contraction between the excited operators, then the structure constant is with the splitting factor defined as The Hexagon function H in this particular configuration is simply related to our determinant expression (2.13) by The Hexagon description of three-point functions allows the evaluation of general configurations where all three operators are excited. If we now let two of the operators haveX excitations, while the other is composed of X fields, then the structure constant becomes where l ij denote the number of Wick contraction between operators i and j at tree level and the sum over partitions is further restricted by the fact that H(α|β|γ) is non-vanishing only when the cardinality of β matches that of α ∪ γ. It is interesting to note that while (2.17) is given as a sum over partitions of three sets of rapidities, a naive tree-level evaluation would give rise to geometric sums naturally yielding a sum over five partitions The equivalence of these descriptions follows from the following tree-level relation We should stress that, computationally speaking, (2.17) is not necessarily a more efficient version of (2.18) as the objects H(α|β|γ) do not have a known compact determinant description. Nevertheless, while the computational gain might not be considerable, there is a conceptual advantage due to the fact that the Hexagon functions can be bootstrapped. First, they obey the Watson equation Together with the diagonal symmetries of three-point functions, these form-factor axioms allowed the determination of the hexagon functions at any value of the coupling [9]. With this in mind, we can attempt a similar rewriting of the dilatation operator overlaps. It is useful to work with normalised spin-chain states where we divide by the norms of on-shell Bethe states {p} = {p}|{p} . These can be conveniently calculated using the Gaudin formula (A.16) which for the coordinate Bethe-ansatz normalisation is with φ defined in (1.23). This can be combined with the normalisation factors of the overlaps to define a new normalisation factor where we used the fact that solutions of the Bethe equations are invariant under complex conjugation, e.g. {u * } = {u}, [41], and the cyclicity condition to simplify the expressions 2 . The overlap with normalised external states can then be written as where H is the same as in (2.19), and we define the new functions H − i as We remind the reader that µk denotes the set of rapidities µ without µ k , following the notation introduced earlier in (1.33), while the short-hand notations for products are defined in (1.32). This decomposition of the overlap in (2.24) seems to fit the splitting of Figure 1(b) particularly well. By cutting the pair of pants depicted in that figure, one would naively expect the side facing away to be represented by the original hexagon H of (2.19), while the side facing forward should lead to something new as it contains the action of the commutator from the dilatation operator. Similarly, the overlap V + can be rewritten as where the normalization is now N + = N /S, with S a symmetry factor that equals 2 when the states in the double-trace are the same, and 1 otherwise, and we have further defined the functions Unfortunately, in this case we are not able to write any of the new objects in terms of the original hexagon function H, since the partitions of the rapidities β into µ and ν appear with a distinct structure. The decomposition is however very similar to that of V − , and seems to match once again the intuition derived from Figure 1(a), with the cutting producing a product between a simpler structure with a more complex ones. We wish to emphasize that while these formulae appear quite involved, they can be straightforwardly evaluated once the rapidities are known, by using, for example, Mathematica. While these expressions are a post hoc massaging of the expressions in (2.9,2.10), when written in this form they clearly resemble the formulas for structure constants. Importantly, the new objects H + i and H − j also obey the Watson equations (2.20) and decoupling conditions (2.21). Note that these properties of H + i and H − i follow from those of the object H(α|β) defined in (2.16). This is non-trivial nevertheless, as it occurs only for certain functions of the rapidities in the summands of (2. 25,2.27). This hints at the possibility that the non-planar dilatation operator overlaps can be written in terms of Hexagon-like objects and potentially determined even at higher orders in perturbation theory. The corrections to the energies of single-trace operators are obtained through the action of H + and H − and a sum over intermediate double-trace operators, which has the natural representation of cutting the torus into two pairs of pants. The fact that the overlaps themselves have a decomposition into hexagon-like objects therefore seems to indicate a possible tesselation of the torus as depicted in Figure 2.
There is an implicit notion of crossing that comes with the decoupling condition. It is natural to imagine that, once such an operation is defined, the excitations can be moved around, so that we relate the hexagon-like objects to a single function where all rapidities are on the same edge. It is upon crossing of the excitations in (2.21) to the same edge that a particle-antiparticle pair X(u 2γ )X(u) can form in a manifest way and decouple from the corresponding form factor. Such a formulation of H + i and H − i with all excitations on the same edge would also be the ideal setup for implementing a bootstrap of those objects. Unfortunately, crossing operations do not commute with the perturbative expansion, and since our one-loop analysis gives access only to the more complicated form of these objects, we were not able to explore further the possibility of such a bootstrap programme. Figure 2: Both V + and V − can be seen as a pair of pants where the asymptotic regions correspond to the three distinct traces involved in the overlap. We have found that each of them can be decomposed into hexagon-like objects satisfying the Watson and decoupling conditions. By glueing them together one can reconstruct the torus, thus finding the non-planar corrections to two-point functions.

Anomalous Dimensions from Overlaps
Our main goal in calculating the above overlaps is to perturbatively compute the leading non-planar correction to operator anomalous dimensions. The general idea is to apply first-order quantummechanical perturbation theory. We denote the planar energies as E (0) and their non-planar corrections at order N −k as E (k) . Given a single-trace operator characterised by momenta {p(u)} solving the Bethe equations and with planar energy E (0) ({p}), the non-planar correction is . (2.28) The sum over I is taken over all intermediate double-trace states where we must sum over all lengths 1 < L q < L p − 1 and for each length sum also over all solutions {q}, {r} of the Bethe equations corresponding to operators of lengths L q and L r with planar energy As a simple example let us consider the unprotected operators of length six in the [2, 2, 2] SO(6) representation. There are two single-trace operators with planar energies and rapidities given by both of which mix with the double-trace operator with rapidities u (4,2),1 = −u (4,2),2 = 1/2 √ 3 and planar energy E (0) (4,2) = 12. The overlaps can be simply found from the general formulae (2.24) and (2.26) and V + (u (4,2) , ∅; u (6,2a) ) = V + (u (4,2) , ∅; u (6,2b) ) = 6 √ 2. The resulting non-planar corrections are These results are in agreement with those found by direct calculation [42], and follow from diagonalising the dilatation operator [3].
In the case where there are only two magnons we can in fact solve the Bethe equation for any length, L, if we consider only cyclic solutions with j e ip j = 1. In terms of the momenta such solutions are given by with n ∈ Z and 0 < n < L−1 2 . Given such a complete set of solutions it is possible to numerically carry out the sum over intermediate states so that we can compute quite efficiently the corrections to energies even for states with quite long lengths, e.g. n = 1 for L = 100, 250, 400, which to six digits gives (2.34) From this and similar numerical examples it can be seen that the corrections to the energies of long operators scale as L 2 /N 2 . This is essentially the well-known BMN limit [30] where one considers operators with large R-charge, J. The non-planar corrections to two-magnon states in the BMN limit were computed in [13,23], see also [43,15] and shown to be It is straightforward to check that our general expressions reproduce this result by substituting the two-magnon rapidities, u n,1 = −u n,2 = L−1 2πn , into (2.24) and (2.26) and taking the large L limit. We must consider the overlaps with all double-trace operators consisting of a vacuum state of length (1 − r)L and two-magnon states with rapidities u m,1 = −u m,2 = rL−1 2πm . Following [15], we then expand in L, sum over m = 0, . . . , ∞ and approximate the sum over intermediate lengths by an integral over r from 0 to 1. At leading order in J = L − 2 this reproduces (2.35), while at subleading orders we find the same result but with J replaced with L−1 = J +1 which is the natural parameter from the perspective of the Bethe equations.
It is naturally interesting to consider higher numbers of excitations. For example at L = 7 with three excitations, i.e. for states in the [3, 1, 3] representation, we have two single-trace operators with planar dimensions E (0) (7,3a/b) = 10. Due to the degeneracy of the states, a naive application of relation (2.28) will fail as it is not clear which linear combination of the Bethe states to use as planar eigenstates. We may use the fact that the two degenerate states are distinguished by their transformation under the parity operation [3]. This operator, P, reverses the order of fields within each trace, for example and commutes with the complete non-planar dilatation operator. Thus the non-planar eigenoperators must have definite parity, and consequently also their planar limits. The rapidities for the two L = 7 and M = 3 Bethe solutions u (7,3a) and u (7,3b) can be easily found using the method (and Mathematica programme) of [44]. They can be seen to transform into each other under parity which acts on finite rapidities by u i → −u i while rapidities at infinity are left invariant. The two parity eigenstates can then be formed from the corresponding Bethe eigenstates as |± = 1 √ 2 (|u (7,3a) ±|u (7,3b) ). Having identified the proper planar linear combinations, we can proceed by computing the mixing with double-trace operators. We choose as our basis of double-trace operators where we have labelled the Bethe states by the magnon rapidities rather than the momenta and Both of these operators have positive parity and the linear combination 2 3 |u (5,3) 5 |∅ 2 − 1 3 |u (5,2) 5 |∞ 2 is the remaining non-protected operator in the [3,1,3] representation. The other linear combination is a descendant of a two-excitation double-trace operator from [2,3,2]. The non-vanishing overlaps following from (2.24) and (2.26) are while the overlaps involving |u (4,2) 4 |∞ 3 are all zero, which is expected since primary operators cannot mix with descendants. Now by applying (2.28) for the parity eigenstates, we find that the non-planar corrections arise from the mixing of the positive-parity eigenstate |+ with the doubletrace state (which has planar energy E (0) = 8) and are given by which agrees with [45,3]. The occurrence of degenerate parity pairs in the planar limit is quite general and so to use non-degenerate perturbation theory we must work within sectors of definite parity. Unfortunately, as has been noted by several authors, for example in [23], [15], [3], the energies following from the Bethe equations demonstrate an additional degeneracy which is relevant to the mixing problem between multi-trace operators. For example, if we consider the two-excitation states with Bethe solution (2.33), states with different lengths, L a , L b , and mode numbers, n a , n b , but equal ratios na La−1 = n b L b −1 will have equal energies. Correspondingly, in the planar limit, the single-trace state corresponding to the spin-chain state |{ 2πm While it is less straightforward to show, analogous degeneracies generally also occur for higher excitation numbers. As just one example, if we consider the L = 8 operators with M = 3, we see that there are three solutions to the Bethe equations. Two of which, whose rapidities we denote u (8,3a) and u (8,3b) , are degenerate parity pairs with energy E (0) (8,3a/b) = 8, while the third is a singular solution with energy E (0) (8,3s) = 12. There is in this case a positive parity double-trace state which is degenerate and which mixes with the positive-parity linear combination of single traces. The mixing matrix can be computed from the overlaps and is from which we can compute the leading corrections to the energies E (1) = ±8 √ 2. We can now proceed to use the corresponding eigenstates to find the subleading 1/N 2 corrections. As we proceed to longer lengths and more impurities, the need to diagonalise the mixing matrix will rapidly become difficult. One way to avoid this problem is to deform the theory to remove such degeneracies. In principle, if we can completely solve the deformed problem, one can then hope to remove the deformation parameter however as this requires resumming the 1/N corrections before removing the deformation we will only be able to make preliminary steps in this direction.
There is another reason for considering the deformed theory which has to do with the singular solutions of the Bethe equations. Already at L = 6 and M = 3 there is a solution u 1 = i/2, u 2 = −i/2 and u 3 = 0 for which the Bethe wavefunction is singular and naive application of the above formulae will lead to unphysical infinities. It is possible to regularize the Bethe equations by the introduction of a twist, see [31] for a useful discussion and further references, which is equivalent to the deformation parameter we introduce below. We can use the solutions of the twisted Bethe equations and the overlaps of the deformed theory to compute non-planar energies which reproduce the undeformed results in the limit of vanishing deformation.

β-deformed SYM Theory
We now turn to the β-deformed N = 4 SYM theory preserving N = 1 supersymmetry. The theory's Lagrangian may be obtained from the undeformed Lagrangian by replacing all products of fields by a Moyal-like -product where the non-commutativity occurs in the internal SU (4) R-symmetry directions [46]. Using N = 1 superspace formalism this corresponds to adding a singletrace deformation to the superpotential, however when written in terms of the component fields this results in both single-trace and double-trace deformations of the Lagrangian [47,48]. For this theory the U (N ) gauge group is no longer conformal at the quantum level due to the couplings of U (1) scalars. These degrees of freedom decouple at the infrared fixed point corresponding to the SU (N ) theory and we will thus consider only the SU (N ) gauge group.
In the remainder of this section we use blue colour to denote how the β-deformation changes terms existing in the undeformed theory, and use purple to emphasize terms that are new.

β-deformed Dilatation Operator
The planar dilatation operator for the deformed theory has been previously studied using both integrable methods [27] and direct field theory computations [49]. The non-planar dilatation operator can in principle be directly computed from the deformed Lagrangian using standard Feynman diagrammatics or perhaps more efficiently using on-shell methods [50]. We instead fix its form by using symmetries and known one-loop results. The form of the single-trace part of the dilatation operator is simply inherited from the undeformed theory and is fixed by the planar theory. It is found by replacing the commutators in (1.11) by the β-deformed commutator [., .] β defined via the R-charges of the fields. In the su(2) sector spanned by X = φ 14 and Z = φ 12 the only relevant commutator is This is supplemented by a double-trace term which is necessary to make the theory exactly conformal [48]. The form of this term follows from the deformed action [47,48] and in the su(2) sector it becomes We fix the coefficient of this term by imposing that the operator Tr(XZ) is a protected operator This has been shown perturbatively at one-and two-loop level by direct calculation [51]. Using these conditions we find that the deformation leaves the tree-level dilatation operator (1.10) unchanged, while the one-loop correction (1.11) gets deformed to It is important to note that although the double-trace term is suppressed by 1/N , it can be relevant at leading order when acting on short operators and results in the vanishing anomalous dimension of Tr(XZ). For longer operators in the planar limit this term is not relevant, however it is essential in understanding the non-planar corrections. The fusion and splitting formulas (1.9) imply for the action of the one-loop dilatation operator on single-trace states where the double-trace part of the dilatation operator contributes the triple-trace term at order 1/N 2 . For the action of the dilatation operator on double-trace states we find Relations (3.5) and (3.6) suggest that the deformed one-loop dilatation operator can be decomposed into planar and non-planar pieces similar to the undeformed case (1.13), however we now find subleading contributions and so we decompose (3.4) as As for the undeformed case H β has a contribution which leaves the number of traces unchanged and so we have diagonal overlaps which we consider in Sec. 3.3.

Deformed Planar Theory
The action of the planar dilatation operator on single-trace operators of length L > 2 is quite similar to the undeformed action and is given by It can be related to the integrable deformation of the Heisenberg XXX-Hamiltonian [52] so that the planar spectrum can still be solved using integrability. As the Hamiltonian H D commutes with i σ z i we can still consider sectors with fixed excitation number M = L − i σ z i . The vacuum corresponding to M = 0 is the same as in the undeformed theory, i.e. (1.20), and has energy E (0) (∅) = 0. Similarly, the one-excitation eigenstate is given by the usual Bethe state, but its energy becomes which we see is no longer degenerate with E 0 when p = 0. For more excitations the spectrum and wavefunctions are still given by the Bethe ansatz (1.34) but now with the deformed S-matrix The Bethe equations determining the momenta are as in the undeformed theory (1.23) but with the S-matrix replaced with S β and the trace cyclicity condition still requires that the momenta satisfy exp(i j p j ) = 1. The dependence of the S-matrix on the deformation parameter can be removed by defining the shifted momentap However, this new parametrisation makes the parameter β manifest in the Bethe equations and cyclicity condition. Introducing the rapidity variable u = 1 2 cotp 2 they are given by thus in terms of the rapidity variable both S β and the corresponding function h β can be defined as in the undeformed case, i.e. via (1.28) and (1.29). One consequence of the deformation is that the degeneracy occurring in the undeformed theory between single-trace and double-trace operators (2.40) is lifted. This can be seen directly in the case of single-trace operators corresponding to two-magnon states, |{k, −k} L , by solving the Bethe equations to the first non-vanishing order in the deformation parameter For generic real values of β there will be no integersm and L 1 such that k(m, L−L 1 ) = k(m, L) and hence no double-trace operator will be degenerate with the single-trace operator. While we do not have a similar proof for states with more excitations, direct diagonalisation of the dilatation matrix for operators with short lengths shows that the degeneracy of excited states is lifted in all cases of operators which were unprotected in the undeformed theory. This reduced degeneracy increases the number of operators for which we can compute the non-planar corrections to the energies by using non-degenerate perturbation theory.

Matrix Elements and Dimensions
The action of the non-planar dilatation operator on Bethe states and the corresponding overlaps can be computed by essentially the same method as for the undeformed theory. For the deformed theory, if we wish to compute the corrections to the energies to order O(1/N 2 ) we must consider not only the off-diagonal contributions from H ± β but also the diagonal contributions from H β .
Off-diagonal Overlaps. We can write the overlaps of H ± β using the notation of Sec. 2.1. As for the undeformed theory, the solutions of the deformed Bethe equations are invariant under complex conjugation which can be used to simplify the expressions. For H − β the overlaps are almost identical to (2.9) The function h in this formula has exactly the same form as in the undeformed theory (1.29) when written in terms of the on-shell rapidities, and similarly for the S-matrices implicit in the scalar products. We should remember however that the rapidities themselves depend on the deformation parameter through the Bethe equations. The overlaps of H + β involve additional contributions in the case where one of the traces has length two and are given by As in the undeformed theory, dividing by the norms of the external states we can define the normalised overlaps V ± β .
Diagonal Overlaps. The contribution H β , which does not occur in the undeformed theory, to the dilatation operator (3.7) contains both length-preserving and -changing parts. Here we are interested in the former, since the computation of non-planar corrections at order 1/N 2 to the anomalous dimensions requires solely the diagonal overlap {p}|H where the δ Lp,2 -term arises from the enhanced contribution of the last double-trace term in (3.5).
The diagonal overlap can then be written in terms of ordered partitions (2.5) as Note that P L (z) vanishes for |z| ≥ L, cf. (2.5). In terms of scalar products of normalised Bethe states (2.7) the overlap can be written as where κ and λ are of the same cardinality, which must be smaller than x − 1. Finally we can divide by the square of the norm of the external state, p 2 , to defined normalised overlaps V β ({p}) which can be then used to compute the energy corrections of single trace states.
Anomalous Dimensions. As the deformation lifts many of the degeneracies present in N = 4 SYM we can use the overlaps in the deformed theory to compute the corrections to energies for a wide range of states by using the deformed analogue of (2.28) (3.21) The additional input to such a calculation are the solutions to the deformed Bethe equations. Solving the deformed Bethe equations is generally a non-trivial task, however for short lengths it can be done either for specific numerical values of β or by starting with the undeformed result and perturbatively solving for β 1. The latter is particularly useful when we wish to use the deformation as a regulator of singular solutions of the undeformed Bethe equations. One must be careful with the order of limits as the one-loop anomalous dimensions are functions of both β and N and we may choose to first expand in large N and then small β, E(β N −1 ), or alternatively first in small β and then large N , E(β N −1 ). In general these expansions will not commute. For example, let us consider the L = 6, M = 3 single-trace operator described in the planar undeformed theory by the roots {u 1 = 0, u 2 = −i/2, u 3 = i/2} with planar energy E (0) = 12. This solution is singular as it has rapidities separated by i. It has a vanishing 1/N 2 correction to the energy due to the su(2) symmetry which ensures there is no other operator with which it can mix. We can study the same operator in the deformed theory where the mixing problem is non-trivial and we find from direct diagonalisation that through O(N −4 ), and keeping only the leading terms in the β-expansion, we have In this expression we can see that the leading non-planar term does not reduce to the vanishing 1/N 2 undeformed answer in the β → 0 limit and in fact the 1/N 4 term is singular. There will be additional singular terms at subsequent powers in the 1/N expansion that would need to be resummed to recover the smooth limit. As the wavefunction in the deformed theory is perfectly regular we can use (3.16,3.17) and (3.20) to compute the overlaps, then take the β → 0 limit and use these expressions to perturbatively compute the undeformed non-planar correction. To be explicit, we need to solve for the deformed rapidities to O(β 6 ) to find a non-singular wavefunction since where we have kept only the first leading term in the β-expansion. The deformation is particularly important when calculating the norm of the singular state, which diverges in the β → 0 limit. However the overlaps themselves are smoothly vanishing in this limit and so, consistent with the symmetries, there is no mixing in the undeformed theory. If we instead use these overlaps in the perturbative formula (2.28) we find a cancellation between the powers of β in the overlaps and the energy differences. As the diagonal contribution is of order β 2 , it gives no leading contribution and we find E (2) in agreement with the result from direct diagonalisation. Let us now turn to the L = 8, M = 3 singular Bethe state |u (8,3s) with planar energy E (0) (8,3s) = 12. This operator is not protected by symmetry in the undeformed theory. Instead it mixes with the double-trace operator |u (6,3s) 6 |∅ 2 made up of the length-6 singular state and length-2 vacuum. From directly diagonalising the corrected energies are where we see that due to the degeneracy the correction is O(N −1 ). Due to the singular nature of the Bethe solution we cannot directly use the overlap formulae of N = 4 SYM but again we can compute the mixing matrix using the regularised singular states in the deformed theory. We find that in this case they are non-vanishing in the β → 0 limit and give the correct 1/N corrections. This corresponds to the case where β N −1 . Alternatively, in the deformed theory as the degeneracy between the two states is lifted, with the planar energy of the single-trace state becoming E (0) (8,3s) = 12 − 36β 2 + O(β 4 ), we can use the same overlaps in non-degenerate perturbation theory in the small β-limit with β N −1 to find the 1/N 2 corrections in the deformed theory. The contribution of the overlaps between the two regularised singular states to E (2) (8,3s) is +4/β 2 , i.e. it is singular in the limit β → 0. There are additional overlaps with other double trace states however they are subleading as is the diagonal contribution which is O(β 4 ). Thus for β N −1 the non-planar corrections start at order N −2 but are singular in the β → 0 limit. This demonstrates that, in general, the two limits β → 0 and N −1 → 0 do not commute.

BMN Limit
We will now look at the two-excitation single-trace solutions and their non-planar corrections in the BMN limit [30] of the deformed theory. First of all, we need to analyse carefully the rapidities that solve the Bethe equations. The solutions are periodic in the deformation parameter with period π and they are symmetric around β = π/2. In general the solutions are parametrized by an integer n which is given as 2πn = Lp 1 − i log(S β (p 1 , p 2 )) . that the rapidity for a strictly positive mode number n is given by with momentum conservation requiring the other excitation in the solution to be u −n . Meanwhile, the zero-mode solutions have a distinct expression where the expansion parameter becomes the square root of the length The planar energies in the BMN limit can then be computed through (3.10), yielding Note that despite the unusual expansion of the zero-mode rapidities u ± 0 , the expansion of the corresponding energy is free of any square roots. Finally, while at the leading order the rapidities u ± 0 seem to be only a particular case of u ±n , it is important to note that the expression for the Gaudin norm, denoted here by N ψ = ||ψ|| 2 , differs already at the leading order by a factor of two Now that we understand the behaviour of the Bethe solutions, we can study the non-planar corrections to the energies in the BMN limit. The strategy is to expand the dilatation operator overlaps obtained in section 3.3 and plug them into (3.21) written explicitly as (2) n from the dilatation operator overlaps. We focus here on mode number n = 1, 2 for singletrace operators up to length 100. We observe that for n = 2 the large L limit is approached differently for even-and odd-length operators. However, fitting the two curves with polynomials in 1/L we find that the mismatch in the coefficients starts only at the subsubleading order.
In the deformed theory there are three contributions that we need to consider: thus leading to a further factor of L. Therefore, the one-loop non-planar energies E (2) ψ scale at most as L 2 which, combined with the colour factor 1/N 2 , produces at leading order in the BMN expansion a factor of g g 2 2 , with the relevant expansion parameters introduced in (1.4). While in principle we can find the BMN corrections to the overlaps at any subleading order, the expansion of the energies eventually breaks down due to the approximation of the summation over intermediate states by an integration. More precisely, we find that for mode number n > 1 the sub k -leading BMN correction to the integrand with k > 1 has simple poles at lengths L /L = n /n with n = 1, . . . , n − 1. As explained in Figure 4 this failure of the BMN expansion is in fact expected and agrees with the numerical experiments performed.
Let us then start with the configuration of uneven splitting 1(A). We consider first an external state with positive mode number n while the intermediate double-trace operator has two excitations on the trace of size L = rL and is described by another positive mode number n . While L is smaller than L, the deformation parameter for the double-trace solution is still expanded in terms of the length L of the single-trace operator as in (3.28), so the rapidities, energies and norms of the double-trace states are written as and analogously for the zero-mode expressions. The overlaps in the BMN limit then become H − nn = 32(r − 1)r 3 n 2 sin 2 (πrn)L 2 (n 2 − r 2 n 2 ) 1 + 2 (r − 1)n 2 n 2 − r(n 2 − rn 2 )b 2 rn 2 (n 2 − r 2 n 2 )L + π cot(πrn) (2r − 1)n 2 + (3 − 2r)b 2 nL + . . . , H + n n = 32n 2 sin 2 (πrn)L (n 2 − r 2 n 2 ) 1 + 2r (r − 1)n 2 n 2 − r(n 2 − rn 2 )b 2 n 2 (n 2 − r 2 n 2 )L + π cot(πrn) (2r − 3)n 2 + (1 − 2r)b 2 nL + . . . . The contribution to the non-planar energies is the combination of those two cases, leading to Note that in this particular case the subleading correction of the Euler-McLaurin formula is vanishing. Furthermore, the contribution of the intermediate zero-mode is crucial for the simplicity of this formula, which would otherwise be plagued by more complex functions such as z 0 dt sin(t)/t. Finally, taking b = 0 yields the non-planar correction to the two-excitation energies of N = 4 SYM. In that theory the rapidities are more naturally written as an expansion in even powers of 1/(L − 1), thus matching the result obtained here. As explained above, the integral approximation at the subsubleading order is not well defined for n > 1, and that is manifested here by the presence of poles in the integrand. However, the expression for n = 1 appears to be well defined at any order, thus allowing us to obtain  n is once again O(L 2 ) and so we can safely ignore the contribution of these modes at the order we wish to consider. The non-planar correction to the energy of the zero-mode solution is then given by where the second term of the first line corresponds to the subleading correction in the Euler-MacLaurin formula (3.34). As expected E (2) 0,A vanishes in the undeformed theory, where the operator becomes a descendant of the chiral primary.
In order to study the second splitting configuration 1(B), where H + β leads to double-trace operators with an excitation in each of the traces, we need to consider the single-excitation solution in more detail. The rapidity in that case is given by Note that L here is again the length of the external single-trace operator as this solution turns out to be independent of the length of the trace it describes and is defined solely in terms of the scaled deformation parameter from (3.28). The energy of this state is (3.43) and the norm is simply the length of the operator. It is important to remember that the singleexcitation solution of length 2, Tr(ZX), is an exception to this formula. The operator is in fact protected due to the contribution of the double-trace term to the planar dilatation operator. When the single-trace operator corresponds to a non-zero mode, the overlaps are We can already see that this will contribute at the subsubleading order, so it suffices to consider the leading integral approximation which gives Meanwhile, for an external zero mode, the overlaps are 12r(r − 1) + (1 + 6r(r − 1) + 4r 2 (r − 1) 2 )b 2 π 2 6L + . . . , We now wish to perform the sum over operators as in equation (3.33). Looking at the expressions for the planar energies (3.31,3.43), we see that E (0) 0 − 2E (0) vanishes at leading order, which effectively enhances the leading order non-planar correction of the energy to O(L 1 ) (note that both H + 0 and H − 0 are subleading). However, this reasoning does not apply when one of the traces has length 2 due to the double-trace term of the dilatation operator, and so that particular double-trace operator contributes at a further subleading order. Therefore the non-planar correction to the energy from this splitting configuration becomes where, as before, the second term of the first line corresponds to a non-vanishing subleading contribution in the Euler-McLaurin approximation. Taking b to zero we see that both E (2) n,B and E (2) 0,B vanish as expected. In that limit the single-excitation solutions correspond to descendants of the undeformed theory, and so these splitting configurations are not expected to contribute to the non-planar corrections of energies in N = 4 SYM. Finally, in the deformed theory we should also consider the diagonal contribution of the dilatation operator 2. However, the overlap grows linearly with L, and so its normalized contribution to the non-planar energy starts only at O(1/L), which goes beyond the order we wish to consider here. Therefore, the non-planar correction to the energy of two-excitation single-trace operators is, at O(L 0 ) in the BMN limit, given by the sum of the off-diagonal uneven splittings from (3.38,3.41) with the symmetric ones in (3.45,3.47) (3.48) We can then write the scaling dimensions of two-excitation states in the BMN limit of large R-charge J = L − 2, which for the non-zero modes yields , (3.49) reproducing the result of [13] at leading order, while for the zero mode we have . Berry and Wilkinson [53]). The situation is quite different for systems with additional symmetries, such as integrable systems, where degeneracies occur even as only a single parameter is varied. For the spectrum of N = 4 SYM this implies that operator dimensions depending on parameters λ and N avoid crossing for generic fixed values of N as λ is varied, as was borne out in [34]. In our case, being at one-loop, the λ-dependence is trivial however we can study the spectrum as a function of both the deformation parameter, β, and the rank, N , of the gauge group. By numerically solving for the eigenvalues of specific families of operators we can see the behaviour of the scaling dimensions as we vary β and as an example we consider the length-six, three impurity states in Fig. 5. For finite values of N the energy levels mostly repel and even at points where they appear to come close they do not in fact cross, maintaining a separation of ∼ 1/N 2 . There is one obvious exception which clearly does cross other levels at finite N . This is a double trace state that does not mix with other operators, receives no 1/N corrections and so is effectively uncorrelated with the other states. This is due to the fact that at half-filling the charge conjugation transformation Z ↔ X combined with the parity transformation (2.36) is a symmetry which commutes with both the impurity number and the one-loop non-planar dilation operator. The double-trace operator is the only L = 6, M = 3 state with negative charge with respect to this transformation. This points to the fact that in order to avoid trivial crossings we must consider operators which have the same quantum numbers. At large values of N we can see the appearance of crossings which occur at special values where β/π ∈ Q; for example in Fig. 5 there are crossings in the planar limit at β = π/4 and β = π/6. These points correspond to values where the β-deformed theory becomes equivalent to an orbifold of N = 4 SYM, e.g. [54], which are known to have enhanced structures such as additional regions on the Coulomb branch.
The level repulsion at finite-N suggests the transistion from a quantum integrable model to a chaotic system. To further explore this it is interesting to compute the distribution of level spacings. Given a spectrum of energy levels one can easily show that if we assume they are uncorrelated the spacing of successive levels satisfies Poisson statistics, P P (s) = e −s . That this distribution is a good description of integrable sytems has been numerically shown in a range of models including manybody systems such as the Heisenberg spin chain, the t-J model at its integrable supersymmetric point and the Hubbard model [55]. There are significantly fewer analytical results, however one important result by Berry and Tabour [56] showed that for a "generic" integrable, semi-classical system the distribution of energy levels is indeed Poisson. There are known examples of integrable models for which this is not the case, such as [57] where by considering finely tuned multi-parameter Richardson-Gaudin models, integrable models with non-integrable statistics were found. However in that case even small changes in the parameters resulted in a restoration of the integrable distribution.
It is well known that in Random Matrix Theory (RMT) the joint probability distribution for the eigenvalues, x 1 , x 2 , . . . , x S , of S × S Hermitian random matrices is given by with α = 1, 2, 4 corresponding to orthogonal, unitary and symplectic ensembles, respectively. Furthermore, the distribution of spacings between eigenvalues, normalised so that the mean spacing is unity, can be well approximated by the Wigner-Dyson distribution

(4.2)
A particularly important feature is that when approximating a physical system by random matrices the appropriate ensemble can be determined from the symmetries of the Hamiltonian, such as rotational or time-reversal symmetry, and does not depend on the specifics of the interactions. In particular, for a Hamiltonian with time-reversal and rotational symmetry it is appropriate to choose the Gaussian Orthogonal Ensemble (GOE) corresponding to α = 1. In order to similarly analyse the spectrum of one-loop anomalous dimensions we must first focus on a specific sector comprising states which have the same quantum numbers for any operators which commute with the dilatation operator. For the β-deformed theory we thus consider states with fixed bare dimension, L, and excitation number, M . Additionally we remove all zero energies and we do not consider the sector with M = L/2. The former correspond to protected states whose dimensions are fixed by supersymmetry and, as previously mentioned, the latter contains states which are related by symmetry and so are degenerate. We then numerically compute the spectrum and order the dimensions in this sector E 1 , E 2 , . . . , E S so that E 1 ≤ E 2 ≤ · · · ≤ E S . One important technicality is that the spectra of RMT are normalised so that the mean level density is unity. So in comparing with physical systems the spectrum must be rescaled to remove the overall dependence on the energy. This procedure is called unfolding and as we do not know the mean level distribution we find a way to approximate it. We describe our procedure in App. B, and find that the final result is relatively insensitive to the specific details of the unfolding. After carrying out this step we label the unfolded spectrum of anomalous dimensions from smallest to largest: x 1 ≤ x 2 ≤ · · · ≤ x S . From the unfolded spectrum we estimate the distribution of level spacings between consecutive levels by computing s i = x i+1 − x i , then binning the data and calculating the fraction that occur in each bin. The results naturally depend on the bin size and a choice is made such that small changes do not significantly affect the overall results. The estimate of the probability distribution naturally improves with larger numbers of states and so one must compute the dimensions for relatively long operators. In Fig. 6 and Tab. 1 we present the results for L = 16 states with N = 17 and β = 0.9. By visual inspection it is apparent that the GOE Wigner-Dyson distribution (α = 1) closely matches the data for most values of M . To be more quantitative, one can fit the data to the Brody distribution which is a one-parameter generalisation smoothly interpolating between the GOE Wigner-Dyson distribution (ω = 1) and the Poisson distribution (ω = 0). In Tab. 1 we show the best fit values of ω for different values of M which are generally close to one suggesting that this is the appropriate value for the distribution at relatively small values of N . This fit captures the Gaussian behaviour of the exponential decay of the tail and the fact that the distribution goes to zero as s → 0. If we assume the distribution is of the Wigner-Dyson form we can perform a fit to the general form (4.2) and find the best fit value of α which from Tab. 1 can again be seen to be approximately one. It is clear that the fit is better for higher excitation number as the values for M = 2 are furthest from those of the GOE. The values for M = 0, which are protected operators, and M = 1, which are protected in the undeformed theory, clearly do not fit the Wigner-Dyson distribution but we do not have a clear explanation for why the M = 2 fit is so poor. In general however we find that the Gaussian Orthogonal Ensemble describes the non-planar distribution of energy levels in the su(2) sector of deformed N = 4 SYM. We can repeat the computation for the strict planar spectrum, however in this case there are additional symmetries that we must account for. In particular the number of traces in a given operator is conserved under the action of the dilatation operator and so we must work at fixed number of traces. In the single-trace sector this reduces the problem to essentially that of the integrable twisted XXX spin chain which is known to satisfy Poisson statistics while the multi-trace sectors are uncorrelated tensor products and so also have the same distribution, see Fig. 7. We can again compare the planar spectrum to the Brody distribution (4.3) and we find that for most impurity numbers the fit is best for a value of ω 0, though there are a small number of cases where the value is larger. If we combine the distributions of single-trace states with M ≥ 3 together we find −0.4 ≤ ω ≤ 0.2 depending on how we bin the data while for the double-trace operators we find −0.4 ≤ ω ≤ 0.4 with the results generally close to zero. Thus the spectrum appears to be well described by the Poisson distribution.
One can see how the distributions change as the spectrum transitions from chaotic to integrable by considering large a sequences of values of N . In Fig. 8 we plot such a sequence of distributions of spacings for L = 15 states. For N = 15 we find the expected Wigner-Dyson distribution while for N = 50 and N = 100 we find distributions between Wigner-Dyson and Poisson with ω 0.71 and ω 0.39 respectively while for N = 200 we find ω 0.05 and the distribution appears to be approaching Poisson. However this is not quite the case with the value of ω further decreasing as we increase N giving ω −0.61 for N = 10 6 . This is due to an excess of points occurring toward s = 0 due to the decoupling of sectors with different numbers of traces which, as explained above, should be considered separately.
We can of course ask about the statistics of the spectrum for the undeformed theory. In this case there are additional symmetries even in the non-planar theory and, in order to find the Wigner-Dyson distribution, we must carefully desymmetrize the spectrum. As mentioned previously, the parity operation, (2.36), commutes with the full dilation operator in the undeformed theory and so non-planar eigenstates with different parity are uncorrelated. Additionally the full global su(2) symmetry is present in the undeformed theory and so it is necessary to work with only highestweight states. This means that at fixed length and excitation number there are fewer available states and consequently the statistics are of poorer quality. For example, if we consider positive parity L = 16 states with M = 3 we find only 315 distinct states. Nonetheless, the general features of the Wigner-Dyson distribution can be seen in a plot of the level spacings, Fig. 9, and the best fit to the Brody distribution occurs for ω 0.9. This suggests that the statistics of the undeformed spectrum are similarly described by GOE random matrix theory.

Conclusions
We have considered the problem of computing non-planar anomalous dimensions in N = 4 SYM and its deformations. The approach we have followed involves two steps: first one must obtain the mixing matrix, and then find its eigenvalues with the method of quantum-mechanical perturbation theory.
In this work we have mostly focused on the first half of the problem by finding the matrix elements of the one-loop dilatation operator in terms of the Bethe rapidities. While direct application of the dilatation operator can in many cases yield the mixing matrix in a similarly efficient fashion, our formulas are given in terms of partitions of the Bethe rapidities, and therefore they are especially advantageous when the number of excitations is small. In those cases we are able to easily evaluate the overlaps, even for long operators where direct diagonalisation would be infeasible, and the bottleneck in computing anomalous dimensions is the determination of the Bethe rapidities. While there are tools for efficiently computing such rapidities, most notably the Baxter Q-function method of [44], carrying out the sums over solutions is still non-trivial. At a more conceptual level, we found that the matrix elements can be written in terms of Hexagon-like objects satisfying both the Watson and decoupling conditions. While our methods are not obviously related to the hexagonalization of the torus, this decomposition hints at the possibility of an approach similar to [58], where four-point functions are built through the OPE, but the OPE data itself is computed within an integrable framework. Similarly, the matrix elements of the dilatation operator might have a more general description which determines their form at higher orders in the perturbative expansion. In order to study this further it would be useful to determine the overlaps at higher loops and to investigate if hexagon-like objects can be found for other sectors of the theory.
One issue with our approach to the diagonalization of the mixing matrix is that it assumes a non-degenerate spectrum of excited states. There are however many degeneracies in the planar spectrum of N = 4 SYM, and so we also considered the β-deformed theory where these degeneracies are lifted. A second advantage of the β-deformation is that it provides a useful regularization for the singular solutions occuring in the su(2) sector of the N = 4 spin chain. The action of the dilatation operator in the deformed theory yields several new structures and for the purpose of evaluating 1/N 2 corrections to the spectrum it is necessary to include an additional diagonal overlap and the contribution of the double-trace term. As an application of our method we computed the anomalous dimension of two-excitation states in the BMN limit through subleading order. We extracted the corresponding coefficients from fits to numerical data at lower lengths and found the results agreed with at least 8 digits of precision. As the problem of degeneracies occurs in other sectors of the theory additional twists will be needed. For example to study the sl(2) sector it may be useful to consider the integrable dipole deformation [59].
The problem of summing over intermediate states increases with the excitation number of the operators under consideration and will rapidly become unfeasible. To compute such sums in the deformed theory it would be of great advantage to generalise the Baxter Q-function method for determining rapidities to the twisted case. It may also be possible that the sum over solutions is simpler than the individual terms and that the computational methods based on algebraic geometry discussed in [60] can be fruitfully applied. It will likely be of interest to study the semi-classical limit of the non-planar corrections where both the number of excitations and the spin-chain lengths are taken to be large. For the planar theory this limit proved to be of great use in making contact with the strong-coupling classical string description. One tool to carry out the sum over intermediate states in this thermodynamic limit is the Quench Action [61], where the sum over Bethe solutions is replaced by a functional integral over root densities which can then be evaluated by saddle-point approximation.
There are of course alternative methods for studying non-planar anomalous dimensions such as Hexagonalisation. There are in principle two approaches that can be taken within this formalism. On the one hand, hexagonalisation can be used to compute four-point functions on the torus [11] and OPE limits then used to extract anomalous dimensions. While in this method one can restrict to correlation functions of protected operators, it gives access only to sum rules of OPE data. On the other hand, it is possible to consider the two-point function on the torus [12] by taking the four-point function with two identity operators, and while this approach does not solve the problem of diagonalizing the mixing matrix, pursuing it beyond tree-level might provide an alternative way of finding the matrix elements of the dilatation operator.
In addition to computing specific operator dimensions it is also of interest to understand the general properties of the spectrum. To this end we analysed the distribution of level spacings and found that at infinite N the one-loop spectra of both N = 4 SYM and its deformation were well described by the Poisson distribution which is characterisic of integrable systems. This could be seen to transition at finite-N to the Wigner-Dyson distribution of chaotic quantum many-body systems which suggests that the statistical properties of the finite-N spectrum can be well described by a GOE random matrix model. Quantum chaos has in recent years been studied extensively in the context of the holographic duality between the SYK-model of N (0+1)-dimensional Majorana fermions and Jackiw-Teitelboim gravity on AdS 2 [62]. The distribution of the level spacings for the SYK model was numerically computed in [63], see also [64], and it was shown that it is Wigner-Dyson with all three ensembles, GOE, GUE and GSE, occuring depending on the value of N . It would be naturally interesting to study this chaotic behaviour at higher loop-orders in N = 4 SYM and whether, by the holographic correspondence, we can describe the properties of interacting quantum strings on anti-de Sitter space by RMT.

A Overlaps from the Algebraic Bethe Ansatz
The algebraic Bethe ansatz (ABA), see [65] for introductions, provides a powerful framework for studying integrable systems such as the spin chains arising in the one-loop planar dilatation operator. Of particular interest in this work are the computationally efficient formulae for scalar products of Bethe states [21,22]. These scalar products have previously appeared in the context of N = 4 SYM structure constants and we will mostly follow the conventions of [39].
Central to the ABA approach is the monodromy matrix,T a (u), which is an operator depending on the spectral parameter, u ∈ C, and acting on the tensor product of the L spin-chain sites, (C 2 ) ⊗L , and an extra auxiliary space, V C 2 , labelled by the index a. ConsideringT a (u) as a 2 × 2 matrix, whose entries are operators acting on the spin chain, we can writê The commutation relations of these entries can be found from the relations where the R-matrix, R a 1 a 2 (u − v), is an operator acting on the two auxiliary spaces labelled by a 1 and a 2 and, for the theories we consider, depending on the difference of the spectral parameters u and v. Viewed as a 4 × 4 matrix, mapping (C 2 ) ⊗2 → (C 2 ) ⊗2 , we can write where we have introduced the functions The trace of the monodromy matrix over the auxiliary space defines the transfer matrix,T (u) = TrT a (u), and it follows from (A.2) that transfer matrices with different spectral parameters commute. The Hamiltonian of the spin chain is given by the log derivative of the transfer matrix evaluated at u = i/2 while the higher conserved charges can be found by further expanding the logarithm of the transfer matrix near u = i/2. The eigenstates of the transfer matrix thus simultaneously diagonalise the Hamiltonian and all higher charges. One can define Bethe states as where the pseudovacuum is defined by C(u)|0 = 0 and satisfies The operators B(u i ) can thus be viewed as creating excited states whose relative normalisation is given by, see [39], where we use the product notation (1.32). The dual states in the ABA are defined by where the dual vacuum satisfies 0|B(u) = 0 and 0|A(u) = 0|a(u) and 0|D(u) = 0|d(u) . and, following [40], can be written as a sum over partitions of the excitations. The partitions are defined by splitting each set of excitations, {u} and {v}, into subsets, α ∪ᾱ = {u} and β ∪β = {v}, with the cardinality of α is equal to that of β. The scalar product is then given as sgn(α)sgn(β)d α aᾱa β dβk αβ kβᾱk αᾱ kβ β det t αβ det tβᾱ (A.14) where and sgn(α) is the signature of the permutation required to put α ∪ᾱ into the canonical order {u}. This formula is valid for arbitrary Bethe states, even those whose rapidities do not satisfy the Bethe equations and which are thus said to be "off-shell". In the case where one set of rapidities does satisfy the Bethe equations, they are said to be "on-shell", the formula can be dramatically simplified to the calculation of a single determinant [22]. There is a further simplification when both sets of rapidities are on-shell and equal. In this case, as the set of rapidities is invariant under complex conjugation, the quantity I M is related to the norm of the Bethe state and is given by Gaudin's formula: where φ k is defined in (1.23).

B Unfolding Procedure
For an ordered spectrum E 1 ≤ E 2 ≤ · · · ≤ E N we define the level density function and the cumulative spectral function, or staircase function, We now separate these spectral functions into smooth and fluctuating parts so that for small separations where D = 1/n av (E i ) is the local mean spacing. These new variables thus capture the nature of the spectral fluctuations about the smoothed behaviour. However, without a priori knowledge of the smooth or mean level density for a physical system we must use approximate methods to compute the unfolded spectrum. There does not appear to be an optimal procedure and so we use a relatively straightforward method. We select each n-th energy from the spectrum {E i } and then perform a piece-wise linear interpolation to define I av .
Fortunately the final result does not appear to be particularly sensitive to the choice of method. For example, we took n = 10 but alternative choices such as n = 2, 200, 400 all give similar results, and so the values for the unfolded spectrum are likely reasonably robust, see Fig. 10. The procedure does cut-off the first and last n-elements and so has edge effects, however as we are interested in differences of energies the overall shift has no effect and the differences in the tails of the unfolded distribution do not modify the final results significantly. In fact for the anomalous dimensions the unfolding process has only a very minor effect and could have been neglected.