12 loops and triple wrapping in ABJM theory from integrability

Adapting a method recently proposed by C. Marboe and D. Volin for ${\cal N}$=4 super-Yang-Mills, we develop an algorithm for a systematic weak coupling expansion of the spectrum of anomalous dimensions in the $sl(2)$-like sector of planar $\mathcal{N}$=6 super-Chern-Simons. The method relies on the Quantum Spectral Curve formulation of the problem and the expansion is written in terms of the interpolating function $h(\lambda)$, with coefficients expressible as combinations of Euler-Zagier sums with alternating signs. We present explicit results up to 12 loops (six nontrivial orders) for various twist L=1 and L=2 operators, corresponding to triple and double wrapping terms, respectively, which are beyond the reach of the Asymptotic Bethe Ansatz as well as L\"uscher's corrections. The algorithm works for generic values of L and S and in principle can be used to compute arbitrary orders of the weak coupling expansion. For the simplest operator with L=1 and spin S=1, the Pad\'e extrapolation of the 12-loop result nicely agrees with the available Thermodynamic Bethe Ansatz data in a relatively wide range of values of the coupling. A Mathematica notebook with a selection of results is attached.


Introduction
The discovery of integrability in the AdS/CFT context [1] has led to important perturbative and nonperturbative results for a special selection of quantum gauge models in 2, 3 and 4 space-time dimensions [2]. The most studied examples are N =4 super-Yang-Mills (SYM), a natural supersymmetric generalization of Quantum Chromodynamics (QCD), and the N =6 super-Chern-Simons (SCS) theory in 3d, the so-called ABJM model [3]. Various aspects of these systems were scrutinized with the aid of integrable model techniques, such as the exact world-sheet and spin-chain S-matrix [4,5,6,7,8,9,10,11], the Asymptotic and the Thermodynamic Bethe Ansätze [1,12,8,13,14,15,16] for anomalous dimensions of local gauge-invariant operators, the quark-antiquark potential [17,18] and the expectation values of Wilson loops in the strong coupling limit [19,20]. Three-point functions [21,22] and the generalization, beyond the strong coupling limit, of polygonal Wilson loops [23,24,25] are currently the object of intense research.
At present, considering the results on perturbative [26,27], exact [28,29] and highaccuracy numerical methods [30], the spectral problem associated to the study of anomalous dimensions in planar N =4 SYM, can be considered virtually solved. One of the crucial final steps was recasting the infinite set of Thermodynamic Bethe Ansatz (TBA) equations into a much simpler matrix nonlinear Riemann-Hilbert problem, the so-called Quantum Spectral Curve (QSC) [31,32]. Starting from the TBA equations of [14,15,16], this simplification was highly nontrivial: it required a deep understanding of the branching properties of the TBA solutions [33,34,35] and their link with generalized T and Q Baxter's functions [36,35,37]. Finally, in the paper [26], C. Marboe and D. Volin proposed an iterative procedure for the exact determination of perturbative contributions, which was applied to compute up to 10 loops for many interesting operators. These results correspond to the sum of a huge number of Feynman diagrams and will probably remain -for a very long time -out of reach of the more standard approaches to quantum field theory.
The ABJM theory in 3d is a second interesting example of quantum gauge theory which can be considered to be exactly solvable in the planar limit. An exact description of the spectrum of anomalous dimensions was obtained by combining information from two-loop perturbation theory [38] and on the strong coupling limit, corresponding to the classical limit of type IIA superstring theory on AdS 4 × CP 3 [39,40,41]. This led to a conjecture for the Asymptotic Bethe Ansatz equations [42], describing operators with large quantum numbers, and to a proposal for the vacuum and the simplest excited state TBA equations [43,44]. The TBA formally encodes the full information on the anomalous dimensions in terms of the dressed coupling constant h, the so-called interpolating function which is an essential ingredient of the integrability approach to this model. As first noticed in [45], h(λ) is a nontrivial function of the t'Hooft coupling λ as, for example, it scales differently with λ at weak coupling (h ∼ λ) and at strong coupling (h ∼ λ/2). For the shortest unprotected operator, the so-called 20, the TBA equations were solved numerically at intermediate values of h in [46]. However, considering the parallel progress made for N =4 SYM, these remarkable achievements cannot be considered fully satisfactory and the purpose of [47] was to begin a reduction procedure which ultimately led to the Quantum Spectral Curve formulation of the spectral problem [48]. Starting from the latter set of equations, exact results on the slope function were obtained in [49], together with an important conjecture on the precise dependence of h on the t'Hooft coupling λ.
The main purpose of this article is to show that the perturbative scheme described in [26] can be adapted to the study of the ABJM theory. We will consider a set of sl(2)-like states corresponding to single-trace operators of the form [50] tr D S (1.1) OSp(2|2) subsector. In terms of the global OSp(6|4) symmetry of the theory, they are characterised by the Dynkin labels [L + S, S; L, 0, L]; in particular, the operator 20 carries the charges L = 1, S = 1. The spectrum of sl(2)-like states has been previously studied in [13,51,52,53,54,55,56,46], and explicit results were known up to four loops. Although the current work parallels in many ways the N =4 SYM case, the difference in the analytic properties between the two models makes the study of the ABJM theory nontrivial and still very challenging, as both the iterative procedure of [26] and the associated publicly available Mathematica program had to be revised step by step. The algorithm shows that the answer, at a generic loop order, is expressed as a linear combination, with algebraic 1 coefficients, of alternating Euler-Zagier sums [57,58,59,60]: where a k = 0 to avoid divergence. In contrast, the results of [61,26] show that anomalous dimensions in the sl(2) sector of N =4 SYM can contain only Euler-Zagier sums with all signs strictly positive, or Multiple Zeta Values (MZV). Examples of results can be found in Appendix A. We see the appearance of the simplest alternating sum, ζ −1 = − ln(2), at 6 loops for twist-1 operators; irreducible multiple sums such as e.g. 2 ζ −1,−3 start to appear from 8 loops, corresponding to double wrapping. We point out that the weak coupling expansion of the exact form of h(λ) conjectured in [49] yields only ζ 2n ∝ π 2n terms (see equation (A.2)); the general transcendental structure of the results appears therefore very similar when they are expressed in terms of the t'Hooft coupling. Nonetheless, it is certainly worth studying whether this rearrangement of terms can reveal some indirect evidence in support of the conjecture of [49]. We leave this for future investigations. We performed explicit calculations up to 12 loops for the simplest twist-1 and twist-2 operators, however the procedure is valid for any values of the twist and spin, and can be pushed in principle to an arbitrary order. Let us collect some comments and observations. We find that the anomalous dimensions do not contain terms of transcendentality one. It would be interesting to see whether there are further selection rules for the answers. In the case of N =4 SYM, it was observed in [61] and [26] that, up to 10 loops, only a subclass of the possible combinations of MZVs appear, related to single-valued polylogarithms [63,64], and it is thus natural to wonder whether the same happens for ABJM (either in terms of h or of λ). Another preliminary observation concerns the shift symmetry: it was noticed in [51] that the anomalous dimensions for the operators (L, S) = (1, 2n) and (L, S) = (1, 2n − 1), n ∈ N + are the same, at the level of the ABA, up to four loops. Our results suggest that this may be the manifestation of an exact symmetry, holding for the maximum and next- 1 In this paper, we restrict to states characterised by a rational Baxter polynomial. In this case, the coefficients are actually rational.
2 As a historical remark, this alternating double sum first appeared in the computation of the anomalous magnetic moment of the electron at three loops in QED [62].
to-maximum transcendentality part 3 of the full anomalous dimension, including wrapping corrections. So far we checked this explicitly up to 12 loops for 1 ≤ S ≤ 6. Finally, the analysis of the highest transcendentality part of the anomalous dimensions reveals a significant difference between ABJM and N =4 SYM. In the latter case, such terms are completely determined by the ABA and single-wrapping corrections. This property makes it possible to compute this part of the answer as an exact all-loop expression [61]. In the ABJM case, we discover that the anomalous dimensions are also affected, at all levels of transcendentality, by double and perhaps even higher wrapping corrections 4 . The difference between the two theories in this respect reminds us of the situation with the slope function [68], which is a purely ABA quantity in N =4 SYM but is affected by wrapping corrections in the case of ABJM theory [69,49]. It would be nice to see how the computation of [61] could be generalised in ABJM. The rest of the paper is organized as follows. In Section 2, we review the QSC construction of [48] and setup some useful notation. In Section 3, we describe the logic of the iterative algorithm. In Section 4, we present some checks of the results, against Lüscher corrections and by comparison with the TBA data of [46]. Section 5 contains the conclusions and some comments on future directions of research. In Appendix A we list the results for a selection of operators. Appendices C and D contain technical details completing our description of the algorithm.

The Quantum Spectral Curve equations
In this Section, we will recall the QSC formulation introduced in [48], valid non-perturbatively for finite values of the coupling h, which will be the basis of the algorithm presented in Section 3.
The unknown functions appearing in the problem are defined on the complex domain of the rapidity u. At finite h, they in general exhibit an infinite ladder of square-root branch cuts, which are described more precisely below. We shall denote byf the second sheet evaluation of a function f , obtained by crossing the cut on the real axis, and by f [n] the shifted value f (u + i n 2 ) obtained by avoiding all cuts. The functions {P i } 4 i=0 are particularly simple since, on the principal Riemann sheet, they have a single branch cut running along u ∈ (−2h, 2h). They satisfy the constraint The functions {ν i } 4 i=1 instead display infinitely many cuts for u ∈ (−2h, 2h) + in, n ∈ Z, but have a simple periodicity/anti-periodicity propertỹ

2)
3 An analogous phenomenon may occur for an asymptotic degeneracy of N = 4 SYM considered in [65,66]. The lifting of the degeneracy due to wrapping corrections was calculated in [66], and appears not to affect the highest transcendentality part of the result. 4 We are currently unable to disentangle double from higher orders of wrapping. One could try to tackle this problem using the generalised Lüscher approach described in [67].
where σ = ±1. The necessity to include anti-periodic solutions (with σ = −1) was not discussed in [48], where the QSC equations were first proposed. Notice however that the introduction of this sign is fully consistent with the construction of [48]: in fact, the ν a functions were introduced to parametrize quadratically a periodic matrix µ ab which was proved to satisfyμ ab = µ [2] ab . It will appear very clearly from the weak coupling analysis that, in order to capture all states, one needs to allow both signs when taking the square root of this relation. The value of σ in (2.2) depends on the operator under consideration. For instance, for L = 1 one has σ = (−1) S ; more generally, σ depends on the mode numbers of the state and can be defined in terms of the associated Baxter polynomial, as will be discussed in Section 3.
The complete analytic structure is described by a set of Riemann-Hilbert type relations In addition to the previous constraints, we impose some regularity requirements on the solutions: the P and the ν functions should have no poles on any sheet, and should remain bounded at the branch points (locally they behave like √ u ± 2h). The global quantum numbers, L ∈ N + (twist), S ∈ N + (spin) and ∆ (conformal dimension) for the states described in this paper are encoded in the asymptotic behaviour of the solution at large u, through the relations while the four components of ν are required to have the leading power-like asymptotics The anomalous dimension γ ∼ O(h 2 ) is defined as in ∆ = L + S + γ. For the following, it is important to underline some further algebraic consequences of the equations presented above. First, using (2.3) and (2.1), one finds and, subtracting (2.10), divided by P 12 , to (2.3), shifted by +i and divided by P [2] 12 , we obtain the following matrix TQ Baxter equation: ak (u) P [1] 12 (u) χ kb . (2.12)

Iteration scheme
In this Section we describe the iterative procedure for computing the weak coupling expansion of anomalous dimensions, where, in contrast with the N =4 SYM setup, each term in (3.1) corresponds to two loop orders in ABJM theory (odd loops are vanishing in this model).
The procedure takes as input the integer quantum numbers L and S. It should be noted that there are in general different states associated to a given choice of these charges. Each state is expected to be associated unambiguously to a solution of the 2-loop Bethe Ansatz [1], and is thus identified by the Baxter polynomial Q(u). As will be explained below, the latter is selected at the first iteration of the procedure, from the solutions of the 2-loop Baxter equation (3.20).

Analytic assumptions
The weak coupling expansion is based, as in [26], on a specific ansatz for the form of the P functions at finite values of h. It is enough to specify the ansatz for {P} 3 i=0 , since we will always consider P 4 to be defined by (2.1). Without loss of generality, we shall set A 1 = 1 and A 2 = h 2 . This is possible thanks to the symmetries of the QSC equations, analogous to the H-symmetry described in [31] (see Appendix B for more details). The choice of scaling for A 2 is important, and comes from the observation that, due to the classical value of the conformal dimensions, Introducing the Zhukovsky variable we can always assume that, on the first Riemann sheet characterised by |x(u)| > 1, the P's can be expanded as Above, M L (u) and K 2L+1 (u) are polynomials in u of degree L and 2L + 1 respectively, with coefficients depending on h, where A 0 (h) and A 3 (h) are the same as in (2.6), and we are simply introducing the dependence on h for clarity. It is by determining these coefficients order by order in h that one is able to compute the anomalous dimensions using (2.7),(2.8).
We shall make a crucial assumption on the behaviour of the polynomials, as well as all the coefficients of the expansion at weak coupling: This is the same scaling that turned out to be appropriate for the sl(2) sector of N =4 SYM theory described in [26].
To setup an iterative algorithm, we consider the perturbative expansions of the P and ν functions for h ∼ 0, keeping fixed the value of the rapidity u. These expansions are called normal scaling expansions [26] and read more explicitly where p i is defined by P i = (hx) −L p i . Naturally, the normal scaling expansion of the p's is simply related to (3.4)-(3.7) by expanding at weak coupling both the coefficients and the Zhukovsky variable ). The procedure also computes the normal scaling expansions of p i : It is important to describe how the analytic properties are modified in the perturbative limit. The branch cuts at u ∈ (−2h, 2h) + iZ shrink to zero, and, at every order in the expansion, poles normally appear, correspondingly, at positions u ∈ iZ. However, the underlying analytic structure of the QSC imposes some very nontrivial constraints. In particular, we shall use the fact that at finite coupling the combinations are free of cuts and analytic on the whole real axis. This means that the normal scaling reexpansion of these combinations should not develop a pole around u = 0 at any perturbative order. This set of constraints is used in an essential way in the procedure. Besides, with the disappearance of the cuts, there is no direct relation between the values of p i,ns,k (u) and p i,ns,k (u) at any given order. Nonetheless, the two expansions are nontrivially related through the finite-coupling ansatz (3.4)-(3.7), since, after the replacement x → x = 1/x, the latter yields an expansion for P i , which is expected to converge at least in an open neighbourhood of |x| = 1 [26]. As we shall see, this will provide precious information for the perturbative method.

Summary of a single iteration
Each iteration of the algorithm starts from the knowledge of p 1,ns,n and p 2,ns,n , complemented with limited information on the normal scaling expansion of p 1 and p 2 , and allows to compute every other quantity at the same perturbative level. The input needed on p 1,ns,n and p 2,ns,n consists only in the coefficients of their Laurent expansions around u = 0, up to the O(1) or O(u) term, respectively. For n = 0, all this information can be easily extracted from the the ansatz (3.4)-(3.5), and we find p 1,ns,0 = 1, p 2,ns,0 = 0, (3.15) For the reader's convenience, we list below the main steps of the algorithm: Step 1: Fix the singular parts of p 0 and p 3 At a generic iterative order, the first step is to determine the singular part of p 0,ns,n and p 3,ns,n . Since the P's have a single cut on the principal Riemann sheet, in the normal scaling expansion they can develop pole singularities only at u = 0. From the weak coupling expansion of the ansatz (3.6), (3.7), we see that, once the pole part is fixed, the regular part around u = 0 depends only on the coefficients of the polynomials (3.9). This is very important in order to keep the number of unknowns in the algorithm under complete control at every iteration.
To fix the pole parts, we consider the sum of equations (2.3) and (2.10) for a = 1, 2, taking into account (2.2). Then we obtain Notice that the combinations in brackets are always of the form (3.13), and therefore free of poles at u = 0 at any order. The order-O(h 0 ) values of the p's and hx are also regular at u = 0. Therefore, the pole part of p 0,ns,n (u) and p 3,ns,n (u) can be determined without requiring the knowledge of ν a,ns,n . At leading order, the above computation simply fixes the singular part to be zero.
Step 2: Solve the Baxter equation for ν 1 , with unfixed polynomial part of p 0 ; then fix the polynomial part using analyticity Then we consider the first of the Baxter equations (2.11), reading explicitly At leading order in the weak coupling expansion, the rhs can be dropped since P 2 ∼ O(h 2 ), and (3.19) reduces to the well-known 2-loop Baxter equation : where ν 1,ns,0 (u) ∝ Q [−1] (u) and where M (0) L (u) is the coefficients of h 0 in the weak coupling expansion of M L (u). It is well-known that the (3.20) admits a discrete set of solutions determining simultaneously the transfer matrix eigenvalue T 0 (u) and the Baxter polynomial Q(u) = S j=1 (u − u j ). Physical solutions are also required to satisfy the zero-momentum condition (ZMC) [38], which can be written as The sign ambiguity in (3.22) is an important difference as compared to the N =4 SYM case. Notice that the sign on the rhs of (3.22) must be identified with the value of σ, since requiring the absence of poles in (3.14) at the leading order we obtain precisely In summary, picking a Baxter polynomial corresponds to choosing a state in the theory and unambiguously defines the value of σ to be used in all subsequent iterations. We will discuss shortly how to fix the precise proportionality coefficient α in At higher orders, plugging the previously computed values of p 1,ns,k , p 2,ns,k , k ≤ n, ν 2,ns,l , l ≤ n − 1, together with the pole part of p 0,ns,n into (3.19), we are required to solve an inhomogeneous Baxter equation of increasing complexity. The technique for its solution, coming from [26], is explained in Appendix C. The output of the procedure returns ν 1 as a combination of rational functions of u and generalized Hurwitz zeta functions. The solution depends on the still undetermined coefficients of the polynomial M L,n (u), plus an additional set of integration constants φ per a,k , φ anti a,k (see Appendix C.3 for details). To fix all these parameters, we use the analyticity requirements coming from the regularity of (3.13),(3.14) for a = 1, which can be seen as the higher loops analogue of imposing the ZMC and polynomiality for Q and T 0 in (3.20). This yields a linear system of equations, which determines all coefficients, apart from m (2n) 0 (which is unfixed due to the symmetries of the system, see Appendix B) and another coefficient, which is the analogue of α at leading order, and corresponds to the freedom to re-define the solution as ν 1,ns,n (u) → ν 1,ns,n + α n Q(u). We will show in the next step how to resolve this ambiguity.
Notice that the size and complexity of this linear system increases with the order of iteration. This is one of the delicate parts of the algorithm, since the solvability of the linear system requires to use nontrivial relations between the Euler-Zagier sums that are the numbers generated by the procedure and make up the matrix of coefficients. Luckily, a full set of relations -reducing any multiple sum to a combination of lower weight zetas from a minimal basis -have been tabulated up to transcendentality twelve in the datamine of [60] 5 . For the 12-loop computation, we used these results up to weight nine (which also critically improves the program's speed and memory usage); at least one more iteration should be relatively simple to perform making full use of these decomposition rules.
Step 3: Compute ν 3 , p 2 and p 4 ; fix the constants α n We can now easily determine ν 3 from the first equation in (2.3): where we know every other term at the relevant perturbative order. The second equation from the system (2.4): can then be used to determine p 2 . Considering the small-u behaviour of p 2,ns,n (u), and matching it with the small-u expansion predicted from the previous iteration (see Step 6), it is possible to fix the coefficient α n . In particular, at leading order, the constraint (3.16) allows us to fix 6 Finally, we compute p 4 using the quadratic constraint (2.1). Notice that, since p 2,ns,0 ∼ O(h 2 ), this only requires to know p 3 at the previous perturbative order. 5 In the case of non-alternating sums, these relations are provided up to transcendentality twenty-two. 6 The sign of the square root in (3.25) is unimportant as the QSC is symmetric under νa → −νa.
Step 4: Solve the Baxter equation for ν 2 , with unfixed coefficients for the polynomial part of p 3 ; impose the analyticity constraints We then consider the second of the Baxter equations (2.11): Notice that, even at the first iteration, the equation is inhomogeneous. Besides, its homogeneous part differs from the Baxter equation (3.20) for a sign in front of the transfer matrix eigenvalue. As explained in Appendix C.2, we can solve this equation and determine ν 2,ns,n as a function of the yet unfixed coefficients of the polynomial part of p 3 . By requiring the regularity of the combinations (3.13) and (3.14) for ν 2,ns,n , one again finds a system of linear equations whose solution fixes all of the coefficients of the polynomial K 2L+1 at order h 2n , apart from the three coefficients k 2L , which are left free due to the symmetries of the system (see Appendix B), and the coefficient of the highest power A (2n) 3 . The latter contains the anomalous dimension, and will be determined in the next step.
Step 5: Compute p 1,ns,n and fix the anomalous dimension; compute ν 4 In terms of the previously found solution for ν 2 , we can compute p 1,ns,n using the first equation of the system (2.4): (3.27) By matching the expected leading behaviour at u ∼ 0, we finally fix the last remaining coefficient A (2n) 3 , and determine the relevant term of the anomalous dimension γ 2n+2 from (2.8).
At this stage we can also compute ν 4,ns,n by inverting the second equation of (2.3): σ ν [2] 2 = −P 0 ν 2 + P 3 ν 1 + P 1 ν 4 . (3.28) Step 6: Compute the seed for the next iteration Finally, starting from the knowledge of { p a,ns,n } n j=0 , we shall show how to determine the quantities needed for the next iteration, namely p a,ns,n+1 and the small-u behaviour of p a,ns,n+1 up to the same powers as in (3.16) (a =1,2).
To achieve this, we start by noticing that the ansatz (3.4)-(3.7), which underlines the whole method, is naturally organised in terms of the parameter y = h/x. Besides the normal scaling expansion, it is very convenient to consider also the double scaling expansion, obtained by expanding p 1 and p 2 at fixed value of y. We can summarise the difference between the two expansions by writing normal scaling: double scaling: The seed quantities for the n + 1-th iteration will be computed from the truncated double scaling series: S a,n (y) = n j=0 h 2j p a,ds,j (y), a = 1, 2. (3.31) Let us then show how to obtain (3.31) order by order. The main observation is that, setting y → h/ x(u) = hx(u) in (3.31) and re-expanding at normal scaling hx(u) ∼ u + O(h 2 ), one should find precisely the truncated series for p a,ns (u), up to the same order. The precise relation between the two expansions is 7 [26]: p a,ds,n (y) = p a,ns,n (u)| u→y − p a,ns,n (u) u→y , a = 1, 2, where p a,ns,n (u) is the coefficient of h 2n in the normal scaling expansion of S n−1 (hx(u)).
Relation (3.32) is used in the program to compute S a,n (y) = S a,n−1 (y) + p a,ds,n (y).
, it is enough to re-expand S a,n (h/x(u)) at fixed u to obtain p a,ns,n+1 (u).
Finally, notice that the ansatz (3.4)-(3.7) implies that the small-y behaviour of the double scaling expansion is, for every l > 0, where p a,ns,n+1 (u) can be computed from the expansion of S a,n (hx(u)). This relation enables us to fix the small-u behaviour of p 1,ns,n+1 (u), p 2,ns,n+1 (u), up to and including the O(u 0 ), O(u) terms, respectively, as needed to restart the algorithm.

Testing the results
The simplest unprotected operator belonging to the symmetric sector of ABJM is the 20 with su(2) and sl (2) representatives with charges L = 2, S = 1 and L = 1, S = 1 respectively. By iterating the procedure illustrated in Section 3, we obtained the 12-loop result reported in Appendix A. In the same Appendix and in an attached notebook we also list the dimensions of a few other operators.
Up to 4 loops, we can directy compare these results with the literature [51,52,53,55,56]. In particular, we tested our method for L = 1 and several values of S (S = 1, ..., 20) against the formula of [51,56], where we have used the shorthand notation H a 1 ,...,a k ≡ H a 1 ,...,a k (S) for the generalized harmonic numbers, and the wrapping term W(L, S) was calculated in [51]: The Baxter function to be used in the formula above for L=1 states is The 6-loop contribution, instead, is a new prediction. However, since it is still not affected by double-wrapping (that sets in starting from eight loops for twist-1 operators 8 ), we can test it against the NLO expansion of Lüscher's formula, which is discussed below.
In the case of twist-2 operators, we checked our results, up to four loops, against the ABA prediction [51,56]: (4.4) At 6 loops we have the first appearance of single-wrapping effects, calculated in [56] for S = 2, 4, ..., 30 by using the formula (4.2) with L=2 and Q(u) = 3 F 2 (−S/2, iu + 1/2, −iu + 1/2; 1, 1; 1), (4.5) while the 6-loop ABA result is too long to be written here, but it can be found in [51]. In order to check the 8-loop result, we need to expand further the Lüscher terms, as for the twist-1 6-loop term, involving corrections to the rapidities and the first contribution from the transcendental part of the dressing factor: this will be done in Section 4.1. We found perfect agreement with all these results.

More checks with Lüscher corrections
Single-wrapping effects can be calculated at all orders using Lüscher's method [70], first generalized in the case of N =4 SYM in [71,72]. We shall use the following Lüscher-like formula, introduced for ABJM by [73,74]: where S AA and S AB are the SU (2|2)-invariant S-matrices found in [75], F b = 0, 1 for bosonic or fermionic indexes, respectively, andẼ M (q) is the mirror energy [76,77] E M (q) = 2 arcsinh Using the techniques developed in [72] for the N =4 SYM case, it is possible to calculate the integrand appearing in (4.6) at various orders in the weak coupling expansion. In particular, in the symmetric case A = B in (4.6), we can use S-matrix elements squared with respect to those used in [72] and the same (fused) scalar factor S M,1 sl(2) (q, p i ). At the leading order at weak coupling, (4.6) reduces to (4.2) for given L and S. More interestingly, expanding (4.6) to higher powers of h, we performed further nontrivial checks of the QSCbased computation. As discussed in [72], to go beyond the leading order it is essential to consider also the finite-size corrections to the physical rapidities or momenta, which yield an additional contribution to the anomalous dimension (4.8) The asymptotic momenta p i appearing in (4.6) and (4.8) are fixed by the Asymptotic Bethe Ansatz equations where Θ(x ± k , x ± j ) is the BES dressing phase [8], the corrections δp k 's are determined by the system [72] Lδp where S(p k , p j ) is the scattering phase involved in the corresponding ABA equation, and In the case L = 1, S = 1, the corrections in (4.8) vanish, since the rapidity is fixed at u = 0 by the trivial Bethe equation (x [1] /x [−1] ) L = −1, and the contribution to the energy determined by the ABA is exactly γ L=1,S=1 ABA = √ 1 + 16h 2 − 1 at any order. The only complication at 6 loops comes from the transcendental part of the dressing factor. The latter however can be treated analytically using the techniques explained in [78]. Finally, one gets γ L=1,S=1 6,wrapping = 288 ζ 3 − 768 + 320 ζ 2 + 504 5 ζ 2 2 + 384 ζ 2 ζ −1 , (4.12) while γ L=1,S=1 6,ABA = 256, matching exactly the result reported in Appendix A. We also checked the 6-loop term in the L = 1, S = 2 case, for which we had to calculate the corrections (4.8)-(4.11) to the asymptotic momenta Once the asymptotic energy γ L=1,S=2

6,ABA
= 224 is added, this agrees with the analytic result reported in Appendix A up to 18 digits. To improve the accuracy by reducing the errors coming from the truncation of the sums, we employed the technique described in Section 5 of [79]. First, we derived a large-M approximation for the (symmetrised) integrands appearing in (4.6) and (4.11). After rescaling q → sM and keeping the first six orders, one finds The contribution of the M > 2000 terms can then be estimated by integrating the functions a n (s) and b n (s) and performing the sum from M = 2001 to ∞. Adding this correction, the agreement with the analytic result coming from the QSC method reaches 48 digits. Such a high precision is sufficient, using EZ-Face [80], to recover the rational coefficients of the zeta functions appearing in the QSC result. Starting at order h 8 , Lüscher's formula (4.6) fails to give the complete answer for twist-1 operators, due to the appearance of double-wrapping effects. However, in N =4 SYM one has that the highest transcendentality part of the result is affected only by single wrapping at any coupling [79,61]. Based on this expectation, we computed the ABA and single-wrapping contribution to the maximally transcendental part of the anomalous dimension, at 8 loops and beyond. In the present case, this does not match the result of the perturbative solution of the QSC, showing, as anticipated in the Introduction, that these terms also include double wrapping corrections. To gain more insights, it would be interesting to extend the Lüscher-like techniques to the double-wrapping order, using the method developed in [67].
For twist-2 operators, instead, double wrapping effects start only at 12 loops, and we can still test the 8-and 10-loop results. In particular, for S=2 we obtained (4.16) from the perturbative expansion of the Bethe equations (4.9). Then the ABA contribution to the energy is that, added to the corresponding terms of (4.17), agrees with the result reported in Appendix A up to 25 digits. A 49 digits agreement is obtained by adopting the asymptotic resummation method described above.
Since the complexity and the size of the integrands involved in these calculations increase drastically with the value of S and the loop order, we managed to perform just some lower-accuracy checks for the L = 2, S = 2 operator at 10 loops and the L = 1, S = 3, 4 operators at 6 loops. Keeping the first 200 terms in the sums, we found

Comparison with TBA
The exact anomalous dimension of the 20 operator at finite values of h was studied by F. Levkovich-Maslyuk in [46], up to values of h ∼ 1, by solving numerically the TBA equations proposed in [44]. In Figure 1 is shown the interpolation of the TBA data of [46], together with various truncations of the perturbative result,  , (4.21) which is in remarkably good agreement with the TBA data up to h ∼ 0.4, as can be seen in the figure. This goes beyond the radius of convergence of the perturbative series, whichby analyticity arguments supported by a fair amount of numerical evidence in N =4 SYM [8,81,26] -is expected to be |h| < 1 4 . The approximation is still rather good decreasing the number of loops, e.g., with the optimal choice of the order of the Padé approximant, the relative error in the estimate of |γ(h) − γ ABA (h)|, compared to the TBA value at h ∼ 0.3, is 2.4 × 10 −2 , 1.2 × 10 −2 , 2.1 × 10 −3 for the Padé extrapolation of the 8-, 10-, and 12-loop results, respectively.

Conclusions
The discovery of integrability in the context of the AdS/CFT correspondence has opened the exciting opportunity to solve certain interacting gauge theories in 3d and 4d. Within the powerful integrability setup, exact expressions for particular physical observables are often fairly easy to obtain, revealing the beautiful simplicity of the final outcomes. One of the recent achievements in this research field is related to the study of anomalous dimensions and their weak coupling expansion through the so-called Quantum Spectral Curve.
In this article, relying on the success obtained for N =4 SYM, we have developed an algorithm to solve perturbatively the QSC for the ABJM model, finding a perfect agreement both with the TBA numerical results and with earlier perturbative computations and predicting several new terms. We are planning to optimise the Mathematica code and to render it publicly available.
Along the way, we have recorded some important formal differences between the present case and N =4 SYM. These include the appearance, for the ABJM model, of Euler-Zagier sums with both positive and negative signs, MZVs of even arguments ζ 2n and the experimental observation that the highest transcendentality part of the anomalous dimensions is not completely determined by the ABA and single wrapping corrections but also contains double (and possibly higher) wrapping contributions.
There are a few interesting open problems related to the current work. First, as already mentioned in the Introduction, it would be interesting to understand better the structure of the perturbative outcomes, and investigate whether some patterns emerge. One could also try to use our results to deduce expressions for the anomalous dimensions for all values of the spin, as was done in [27] in the case of N = 4 SYM. This would be particularly interesting for the purpose of exploring various conjectures on the large-spin behaviour [51,82]. The generalization of this perturbative approach to other sectors of the theory is another obviously important problem. Although we do not have yet concrete results to discuss, the generalization to the full OSp(2|2) sector should be fairly straightforward, while setting up the iterative procedure for generic operators appears to be more involved. Finally, it would be very nice to complement the current results by adapting to ABJM the numerical technique developed in [30] for the non-perturbative solution of the QSC, or to try to transfer some of these powerful methods to the study of the QSC for the Hubbard model [83].

A. Sample results
In this Appendix we present some explicit results, which can also be found (including three more operators) in the Mathematica notebook Results.nb attached to the present paper.
Using reduction formulas such as the ones in [84], the result can be written in terms of a small number of non-reducible sums. We will use the basis of [59,84], which is conjectured to be minimal. We find empirically that the l-loop result involves only the basis elements ζ a 1 ,...,a k with weight w = k i=1 |a i | ≤ l − 2 and depth k ≤ l − 1 − w. In particular, the 12-loop result can be written in terms of the following sums: We report the results in terms of the interpolating function h(λ). Assuming the validity of the conjecture of [49] for this quantity, namely, the anomalous dimensions can be rewritten in terms of the true coupling constant λ using the expansion: the first two orders of which have been verified by direct calculations [52].

Twist-1
For the 20 operator with L=1, S=1, Q(u) = u, we found: where for brevity we have grouped together a number of terms that are common to the S=1 and S=2 cases, at orders h 8 , h 10 and h 12 : We shall give the following results up to double wrapping. For L=1 and S=3, Q = u 3 − 5 4 u:

B. Symmetry
In this Appendix we shall discuss a simple symmetry of the QSC equations, which allows us to choose freely the constants A 1 and A 2 and four of the coefficients of the polynomials (3.8). The simplest symmetry of the QSC equations (see [32] and Section 2.3 of [28] for the very similar N =4 case) is the transformation where R is any 4 × 4 constant matrix satisfying To preserve the structure of the QSC, we should also require that the transformation does not change the ordering of the magnitudes of the P functions at large u, The most general form of R compatible with these constraints has six degrees of freedom and can be written as: and the transformation acts on the P functions as From (B.4), it is simple to see that the constants A 1 and A 2 characterising the leading asymptotics P 1 ∼ A 1 u −L , P 2 ∼ A 2 u −L−1 can be fixed to arbitrary values by an appropriate choice of α 1 and α 3 . Similarly, relations (B.5),(B.6) show that the coefficients m 0 (h), k 0 (h), k L (h) and k 2L (h) defined in (3.8) -which correspond to certain coefficients of the large-u expansion of P 0 and P 3 -can be chosen freely by tuning α 5 , α 6 , α 2 and α 4 , respectively. Accordingly, these numbers are not fixed by the algorithm at any order in h.

C. Solving inhomogeneous Baxter equations
In this Appendix, we shall present the basic method to solve the inhomogeneous Baxter equations (3.19) and (3.26) encountered in the procedure.

C.1 Solving the inhomogeneous Baxter equation for ν 1
At a generic perturbative order O(h 2n ), equation (3.19) reduces to where q 1 = ν [1] 1,ns,n , U 1 a source term of increasing complexity, and T 0 is the zero-order transfer matrix. At leading order, we know that the source term is zero, and the regular solution is the Baxter polynomial Q. To solve the generic case, we follow the method of [26]. Considering the ansatz q 1 = Q f [1] 1 , it is simple to see that, in order for (C.1) to be fulfilled, f must satisfy where we have denoted ∇ + g = g − g [2] and ∇ − g = g + g [2] . We introduce the inverse operators Ψ + and Ψ − , such that ∇ ± Ψ ± g = g. We shall give a precise operative definition of the operators Ψ ± in Section D.2, by explaining how they act on the functions generated by the algorithm. Using these operators, an inhomogeneous solution of (C.1) can be found as .
Moreover, setting U 1 = 0 and Ψ − (0) = Φ 1,anti , where Φ 1,anti = −Φ [2] 1,anti is a generic antisymmetric function, we find a second independent solution of the homogeneous equation, (C.4) Introducing for convenience the notation and putting all pieces together, the general solution of (C.1) can be written as , where Φ 1,per = Φ [2] 1,per denotes a generic i-periodic function. Following [26], it is convenient to rewrite this expression in order to cancel its apparent poles at the Bethe roots where Q = 0. To achieve this, we introduce two polynomials A and B, of degree S − 1, defined by the condition The Baxter equation then implies that with R a polynomial of degree L − 2. Introducing the constants r ±,k and the polynomial C through we see that (C.5) can be written as where η −k (u), k ∈ N + are certain polygamma functions (see equation (D.4) below), solutions of Likewise, the inhomogeneous solution of (C.3) can be written in the following form, which only involves poles at positions u ∈ iZ: (C.13)

C.2 Solving the inhomogeneous Baxter equation for ν 2
The equation arising from (3.26) at the (n + 1)-th iteration of the algorithm has the form (C.14) with q 2 = ν [1] 2,ns,n . One can proceed in a similar way as for (C.1), by paying attention to the different sign in front of the T 0 q 2 term. This implies that a simple, homogeneous solution of the equation for U 2 = 0 is simply given by q 2 = Q Φ 2,anti , where Φ 2,anti = −Φ [2] 2,anti is any anti-periodic function. Likewise, we see that an independent family of solutions to the homogeneous equation is described by where Φ 2,per = Φ [2] 2,per is i-periodic and Z is defined in (C.11). Making the ansatz q 2 = Q f [1] 2 , we find (C. 16) so that a solution to the inhomogeneous equation can be found as q 2 = Q f [1] 2,inhomo , with , (C. 17) and the most general solution can be written as .
(C. 18) More explicitly, in terms of the same polynomial C and constants r k,± defined in (C.9), (C.17) can be expressed as

C.3 Periodic coefficient functions
To construct periodic/anti-periodic functions without introducing unphysical poles or an unphysical exponential growth at infinity, we consider the following combinations where η k is defined in (D. for a = 1, 2, where φ per a,j , φ anti a,j are free parameters. The number of terms included in the sum increases with the perturbative order: at the n-th iteration, we may take Λ = 2(n − 1).

D. Functions generated by the algorithm
In this Appendix, we give a precise definition of the operations Ψ ± . This will show that the algorithm always produces answers in a specific algebra of functions comprising: • rational functions of u, with poles only in the set u ∈ iZ ; • the functions η A , with A a multi-index A = a 1 , a 2 , . . . , a k , a i ∈ Z \ {0}, which are a generalisation of Hurwitz zeta functions [61] , • periodic/anti-periodic functions constructed as explained in Section C.3.
The presence of alternating signs is an important difference as compared to the N =4 SYM case, leading to the appearance of alternating Euler-Zagier sums in the results.

D.1 The η functions
When the sums are convergent, the operators Ψ + and Ψ − can be implemented as Ψ + (g) = which is simply expressible in terms of polygammas: It should be noted that the sum (D.2) is divergent for a = 1, but we may use (D.3) as regularised definition of η 1 . The ambiguity in this choice does not affect any physical results. In fact, notice that the η functions enter the algorithm through the solution of inhomogeneous Baxter equations and this ambiguity only goes into a redefinition of the integration constants described in Section C. Another important point to underline is that the particular regularisation of η 1 introduces iη 1 (i) ≡ ζ reg 1 ≡ γ Euler-Mascheroni in the solution of the QSC. While this number enters the expansions of various P and ν functions (which are not directly physical due to the symmetry described in Appendix B), it always cancels out of the anomalous dimensions.
We then define η a,B , with B any multi-index, as follows where by convention we set n 0 = −1. The sum is convergent for a k = 1. One can define the marginally divergent cases where a k = 1 through the stuffle algebra (see Appendix A of [61]), which allows one to express them in terms of convergent η functions and η 1 . Relation (D.6), shows that the Laurent expansion of the η functions around u = i naturally produces multiple Euler-Zagier sums, which therefore enter the algorithm. The precise relation is Each rational function of u is then broken into a polynomial part and a sum of inverse powers, leading to the appearance of η functions Ψ ± 1 u |a| = η ±|a| . (D.12) As in [26], we then define recursively the action of Ψ ± on the products of rational functions and η's. One starts form noticing the following relation (which is a consequence of the nested definition (D.5)) η a,A = (sgn(a)) η [2] a,A + η [2] A u |a| , (D. 13) for any multi-index A, which leads immediately to ±|a|,A . (D.14) Besides, when more general combinations like η B /(u+i k) |a| are encountered, one can use (D.13) until either a term of the form (D.14) is met, or the η runs out of indices (η ∅ = 1). To deal with the polynomial-times-η part, we use the relation b,A ± (sgn(b)) Ψ ±ign(b) (u |a| ) [2] η b,A ∓ (sgn(b))Ψ ± Ψ ±sgn(b) (u |a| ) [2] η [2+2n] A (u + in) |b| .
(D. 17) In conclusion, the operators Ψ ± are fully defined on the algebra of trilinear combinations of rational, η and P k functions. We notice that, since some of the equations adopted in the algorithm are quadratic, products of two η functions may also be generated along the way. However, by using stuffle algebra relations such as the ones in [61], these can be converted into expressions that are linear in all the η functions. Up to the order we reached, it actually turned out that this step was unnecessary, since all the source terms in the inhomogeneous Baxter equations were already linear in the η's.