Exploring the spectrum of planar $AdS_4/CFT_3$ at finite coupling

The Quantum Spectral Curve (QSC) equations for planar $\mathcal{N}=6$ super-conformal Chern-Simons (SCS) are solved numerically at finite values of the coupling constant for states in the $\mathfrak{sl}(2|1)$ sector. New weak coupling results for conformal dimensions of operators outside the $\mathfrak{sl}(2)$-like sector are obtained by adapting a recently proposed algorithm for the QSC perturbative solution. Besides being interesting in their own right, these perturbative results are necessary initial inputs for the numerical algorithm to converge on the correct solution. The non-perturbative numerical outcomes nicely interpolate between the weak coupling and the known semiclassical expansions, and novel strong coupling exact results are deduced from the numerics. Finally, the existence of contour crossing singularities in the TBA equations for the operator $\textbf{20}$ is ruled out by our analysis. The results of this paper are an important test of the QSC formalism for this model, open the way to new quantitative studies and provide further evidence in favour of the conjectured weak/strong coupling duality between $\mathcal{N}=6$ SCS and type IIA superstring theory on $AdS_4 \times CP^3$. Attached to the arXiv submission, a Mathematica implementation of the numerical method and ancillary files containing the numerical results are provided.

was the discovery of a link with integrability, at both weak and strong coupling [4,5], in the planar ('t Hooft) limit of the duality.
The QSC is probably the ultimate simplification of the spectral problem, as it allows to compute numerically the spectrum at finite coupling with high precision [19,20] and to inspect analytically interesting regimes such as the BFKL limit [21,22] or the weak coupling expansions [23,24]. It has also been generalised 1 to the γ and η deformations [26,27], to the so-called fishnet theory [28], to the quark-antiquark potential [29,30] and very recently also to the calculation of correlators of three cusps in a special limit of N = 4 SYM [31].
Despite the considerable progress made in this research field, there are still many interesting open problems and possible generalisations, see e.g. [17]. A practical problem is the fact that, while the QSC potentially allows to study the anomalous dimension of arbitrary operators, it is rather difficult to find starting points ensuring the convergence of the iterative algorithm on a given chosen operator. Usually weak coupling data can be used efficiently as an initial seed for the numerics. An initial step towards covering of the full spectrum of SYM was taken in [24], solving the QSC at one loop for a wide set of states. However, while at weak coupling there exists an efficient procedure to solve the QSC even beyond 10 loops [22,23], at strong coupling a systematic perturbative approach is still missing (see however [20] for progress in this direction).
Planar N = 4 SYM is not the only interesting AdS/CF T -related integrable theory. Further examples, that for different reasons can be considered equally important, exist and concern supersymmetric conformal gauge theories in lower space-time dimensions. As a matter of fact, these models are intrinsically more complicated, are not maximally supersymmetric, and are currently much less understood compared to N = 4 SYM.
The AdS 4 /CF T 3 case -the main subject of the current paper -was introduced by Aharony, Bergman, Jafferis and Maldacena (ABJM) in [32] and is potentially very important since it involves, on the AdS part of the correspondence, a 4D quantum theory of gravity.
In the integrable planar limit the gauge side of the duality corresponds to the N = 6 superconformal Chern-Simons theory, while the gravity side is described by the type IIA superstring theory on AdS 4 × CP 3 [33][34][35][36][37] (see also the review [38]).
In contrast to N = 4 SYM, in the ABJM theory integrability leaves unfixed the interpolating function h(λ) [34,39], which parametrises the dispersion relation of elementary spin chain/worldsheet excitations and enters as an effective coupling constant in the integrability-based approach. An important conjecture for the exact form of this function was made in [40] by a comparison with the structure of localization results. The conjecture was extended in [41] to the more general ABJ theory [42]. The proposal of [41] suggests that the only difference between ABJM and ABJ corresponds to the replacement of h(λ) with an explicitly defined h ABJ (λ 1 , λ 2 ), where λ 1 , λ 2 are the two (apparently) independent couplings of the ABJ theory. Therefore the analysis performed in this paper is also potentially relevant to the more general ABJ model.
Anomalous dimensions of single trace operators with asymptotically large quantum numbers are described, at all loops, by the Asymptotic Bethe Ansatz (ABA) equations, conjectured in [43] and derived from the exact worldsheet S-matrix of [44]. From the S-matrix, it was possible to obtain the leading order finite-size corrections (see for example [45][46][47][48][49]). The exact result, including all finite-size corrections for short operators, is formally described by the infinite set of TBA equations and corresponding functional relations, proposed in [50][51][52]. The latter equations were solved numerically in [53] and the anomalous dimension of the operator 20 [33] computed up to h = 1. The results of [53] were so far the only finite coupling data for the spectrum of this model.
The QSC equations of the ABJM model [54,55] have been already used to compute the so-called slope function in a near-BPS finite coupling regime [40] and to develop an efficient algorithm for the weak coupling expansion in the sl(2)-like sector [56].
The main purpose of this work is to compute the finite coupling spectrum of a set of short operators by solving numerically the QSC equations. The TBA results of [53] for the operator 20 are extended far beyond h = 1, until the dual string description starts to emerge clearly, and a finite coupling analysis of other states in the sl(2|1) sector is performed. This work can be considered as the first step toward a more systematic study of the ABJM finite coupling spectrum [57]. As an attachment to the arXiv submission, we provide a simple Mathematica implementation of the numerical method for the subsector of parity-even operators; more general versions of the code are available upon request.
The rest of the paper is organized as follows. Section 2 contains a review of the basic QSC equations derived in [54,55]. The structure of the numerical algorithm is schematically described in section 3 where the differences with respect to the N = 4 SYM case are underlined. In section 4 the numerical results obtained for two of the simplest and most studied operators in the sl(2)-like sector are reported and compared with existing finite coupling results [53] and strong coupling predictions [40,58]. Furthermore, the analytic structure of the function Y 1,0 associated to the operator 20 is carefully investigated, extending the computation performed with TBA techniques in [53] to larger values of the coupling constant and showing the absence of critical values of h for this state. See for example [59,60] for a discussion concerning the possible emergence of critical values of the coupling constant in N = 4 SYM.
In section 5, the numerical analysis is extended to a couple of operators that do not belong to the sl(2)-like sector. One of the novel feature here is the emergence of a nontrivial h-dependent phase in the QSC equations, i.e. P(h), which is explicitly computed both numerically at finite coupling and perturbatively at weak coupling. The results of this paper are an important test of the self-consistency of the QSC formulation of [55] even for these more complicated operators.
The paper ends with a series of concluding remarks and four appendices which contain technical details about the symmetries of the QSC equations, the reconstruction of the TBA solution from the Q functions, analytic weak coupling expansions and the numerical results for ∆(h).

Review of useful equations
In this section we review the basics of the QSC formulation presented in [54,55]. It involves a large number of Q functions, depending on the spectral parameter which we denote as u. Among them a primary role is played by the Q functions denoted as P and Q, which can be viewed as a quantum version of the quasi-momenta parametrising classical string solutions in AdS 4 × CP 3 . P and Q functions roughly correspond to CP 3 and AdS 4 degrees of freedom, respectively.
The P functions enter a self-consistent formulation of the spectral problem, the Pνsystem, which is a closed set of discontinuity relations -a non-linear Riemann-Hilbert problem -involving a finite number of unknown functions defined on a Riemann surface with an infinite number of sheets. On the reference Riemann sheet, the P's have a single, square-root type branch cut running from −2h to 2h. The ν's instead are required to fulfil the following quasi-periodicity condition: where f (u) denotes the analytic continuation of f (u) to the next sheet through the cut u ∈ (−2h, 2h). In (2.2), P(h) is a state-dependent phase that may be, in general, a nontrivial function of the coupling constant h, as will be discussed in more detail in section 5.1.1. Setting with (P 1 , P 2 , P 3 , P 4 , P 5 , P 6 ) = (−P 4 , P 3 , P 2 , −P 1 , −P 6 , −P 5 ) and the constraints P ac P cb = δ ab ←→ P 5 P 6 − P 2 P 3 + P 1 P 4 = 1 , ν a ν a = 0 , (2.4) the Pν-system is written as Notice that, as a consequence of these equations, the P functions exhibit an infinite sequence of evenly-spaced short cuts from the second sheet onward with branch points located at u = ±2h + iZ (see figure 1). The ν functions possess evenly-spaced short cuts with branch points at u = ±2h + iZ on all the Riemann sheets. The P's and the ν's are required to be bounded and free of singularities, other than the branch points at u = ±2h + iZ, on every sheet of the Riemann surface. In addition, equations (2.4)-(2.6) need to be supplemented with the following large-u asymptotics of the P functions: where 8) with no contraction over the index B. In (2.7) and (2.8), the charges M andM corresponding to a given state can be identified as (see [55] for more details) where L is the spin chain length, γ is the anomalous dimension and K i the excitation numbers in the ABA description of the states [43,44], in sl(2) grading. In particular, the following ordering holds: which implements unitarity of the representation of the superconformal algebra. An important consequence of the analytic properties of the P functions, which is fundamental for the numerical algorithm, is that they admit a convergent series representation: where x(u) is the Zhukovsky variable defined as . (2.13) In principle, the set of equations (2.4)-(2.6) already contains all the information necessary to compute the planar AdS 4 /CF T 3 spectrum at fully non-perturbative level. However the currently available algorithms for the numerical solution at finite coupling of the AdS 5 /CF T 4 spectrum [19,20] c|i . (2.14) The Q a|i 's are required to be analytic in the whole upper half plane (they turn out to have an infinite ladder of short branch cuts in the lower half plane starting from Im(u) = −1/2), to have power-like asymptotics at large u and can be normalised as where κ ij is an anti-symmetric matrix, independent of u, whose only nonzero entries are κ 14 = κ 32 = −κ 23 = −κ 41 = 1. They can be used to construct the Q and τ functions -the AdS 4 counterpart of the P and ν functions -as The corresponding Qτ -system is with τ i ≡ e −iP κ ij τ ++ j , and κ ij is the inverse of κ ij . The cut structure of the Q functions is represented in figure 1, while the τ 's inherit from the ν's the infinite set of evenly-spaced short cuts at u = ±2h + iZ and have the 2i-periodicity property τ [+4] i = τ i . Below, we will use a more convenient vector notation for the Q functions: To summarise the large-u asymptotics of the main Q functions it is convenient to introduce the combination of charges 19) which allow us to write Finally, upon a specific choice of basis for the solutions of the system (2.14), the Q functions and their analytic continuations fulfil a further set of constraining equations, the so-called gluing conditions, 3 which were derived in [55] for the ABJM model in the case of halfinteger spin and real values of h. These equations are the main ingredient of the numerical algorithm. For the validity of the gluing conditions, we choose a solution of (2.14) with a particular large-u asymptotic expansion of the "pure" form [22,55]: where, for the purposes of this paper, all coefficients B (a|i),n are real. Then the gluing conditions are (see [55]): where the coefficients δ 1 , δ 2 are constrained by

The algorithm
This section contains the description of the algorithm implemented for the numerical computation of the AdS 4 /CF T 3 spectrum, using the equations reviewed in section 2.

General setup
Following [23], we expand the P functions in powers of the Zhukovsky variable x(u) around x = ∞ with coefficients {c A,n } n≥0 as in (2.12). In analogy with the N = 4 SYM case [19], it is simple to deduce from the analytic properties of the P functions that (2.12) converges everywhere on the first sheet of the uplane, which corresponds to |x(u)| > 1 in figure 3, and also in an elliptic region around the cut u ∈ (−2h, 2h) on the second sheet (see figure 2). The convergence region is indeed bounded by the position of the nearest branch points, which lie at u = ±2h ± i on the second sheet. On the other hand, an analogous Laurent expansion in x(u) for the Q functions would not even converge on all points of the cut u ∈ (−2h, 2h). This is the reason why one starts with a series representation for the P's and not for the Q's.
As will be discussed in the following sections and in [57], the convergence of (2.12) is a particularly delicate issue at strong coupling and when analytic continuation to complex values of h is numerically implemented. Since in the current context we will mainly adopt the x-plane perspective, it is useful to introduce here some general concepts. Formula (2.13) maps the complex plane of u into the complex plane of x according to the following rules and the whole second sheet in the u-plane is mapped into the interior of the unit circle. Thus the series (2.12) converges everywhere in the x-plane except for the inner green disk in figure  3 containing all the branch points belonging to the second sheet. Since the gluing conditions are evaluated on the cut between −2h and 2h (unit circle in figure 3), it is crucial that the unit circle falls completely inside the convergence region of the P functions. As h is increased along the real axis, the branch points approach x = ±1 and correspondingly the numerical algorithm, which is based on a truncation of the series (2.12) becomes less and less efficient, since the cutoff N A should be accordingly increased with h. The final objective is to determine -with very high numerical accuracy -the set of coefficients {c A,n } N A n=0 solving the QSC equations together with the gluing conditions. Since the conformal dimension ∆ is related to c A,0 via (2.7)-(2.10), in the following we will work with the equivalent set of unknowns Schematically, the algorithm consists of two main blocks. In the first block Q I (u) ≡ Q I (u; X) and Q I (u) ≡ Q I (u; X) are formally evaluated on the cut u ∈ (−2h, +2h) in terms of the parameters X. This enters as a subroutine in the main part of the program, in which X is fixed by imposing the gluing conditions on the cut; for this purpose, the iterative Levenberg-Marquardt procedure is used.
Part 1: computing Q I and Q I . The strategy to compute Q I (u) and Q I (u) on the cut u ∈ (−2h, +2h) is the following. Starting from the ansatz (3.1) for the P's, first we compute Q + a|i (u) on the cut by solving (2.14). Below, the method to solve such a finite difference equation is explained in detail. Then, using formula (2.16) and its analytic continuation we find Q ij (u) and Q ij (u). The computation of the Q functions is therefore reduced to the solution of equation (2.14) for Q + a|i (u) on the segment u ∈ (−2h, +2h); this is done in a two step calculation. First we find an approximate solution to (2.14) at some u with large integer imaginary part Im(u) = N s . For this purpose, we truncate the large-u asymptotic series representation (2.21): which allows to reduce the finite difference equation (2.14) to a much simpler linear system 4 for the unknowns B a|i,n N a|i n=1 , where the leading order coefficients B a|i,0 are known up to a gauge choice. Then, iterating N s times equation (2.14), the large-u solution can be shifted down to u ∈ (−2h, +2h): Part 2: fixing the coefficients X. Here the strategy is to fix X by imposing the gluing conditions on the cut u ∈ (−2h, +2h). A suitable functional F( X) is built out of the gluing conditions, in such a way that it vanishes when the latter are satisfied. Then X is obtained looking for a root of F( X). We begin by rewriting (2.22)-(2.24) in a discretised form together with their analytic continuation f i (u k ; X) obtained by replacing Q I → Q I and Q I → Q I . In (3.6)-(3.11) the coefficients are defined as give a discretisation of the interval (−2h, 2h). For numerical convergence the optimal choice of discretisation points corresponds to the zeros of the Chebyshev polynomials of the first kind adapted to the interval, i.e. u k = 2h cos( 2k−1 2Np ) , (k = 1, . . . , N p ). Notice that the gluing conditions (2.22)-(2.24) are defined up to the parameter δ 2 , which cannot be fixed a priori. This ambiguity is lifted by treating it as a genuine additional unknown: leaving equations (3.6)-(3.11) formally unchanged 5 . The next step is to arrange (3.6)-(3.11) into a (12 N p )-dimensional vector where the generic element is labelled by the multi-index I = (i, k). It is then natural to define the functional F( X) as the squared norm of the vector f ( X), i.e.
is a real and positive defined quantity which vanishes when the gluing conditions are fulfilled.
The Levenberg-Marquardt method appears to be the right choice for a minimisation problem of this kind, as observed in [19,20]. To implement this iterative procedure efficiently, an initial guess -close enough to the solution -for the parameters X (0) is needed. For small values of h up to h 0.30, analytic data from the weak coupling expansions 6 provide a good starting point for ∆ and {c A,n } N A n=1 . Whereas to move outside the weak coupling regime in practice it is necessary to change the coupling in small steps, using an extrapolation to obtain the initial configuration for a given value h.
Similarly to [19], for most of the operators discussed in the following, we find that for the convergence of the algorithm it is sufficient to impose the validity of a subset of the gluing conditions: in particular, f i (u k ; X) with i = 3, 5, 6 are found to be sufficient in most cases. Instead, an additional gluing condition, i.e. f 2 (u k , X) or equivalently f 4 (u k , X), is necessary for the numerical convergence of the method in the non-symmetric case discussed in section 5.2. This is so far an experimental observation and it would be interesting to clarify why this is the case.
Concerning the convergence of the iterative algorithm, it turns out to be very important to reduce the space of the parameters to a submanifold where the solution of the QSC is non-degenerate. For example, the quadratic constraint (2.4) was imposed at each iteration of (2.14) by considering P 4 as a function of the other five P functions throughout, therefore reducing the set of independent parameters to {c A,n } N A n=1 with A = 4. Furthermore, the QSC admits a continuous family of symmetries, which implies that infinitely many different sequences of coefficients {c A,n } may be used for the description of the same physical state.
The gauge fixing of these extra symmetries is discussed in detail in appendix A. 5 In [19] a slightly different method is used to remove a similar ambiguity in the N = 4 SYM gluing conditions, namely the unfixed parameters are evaluated by constructing normalisation independent combinations of Q functions. While this is perfectly equivalent, we found the method described here to be numerically more stable in some challenging regimes, such as at strong coupling and in the proximity of the branch point at h ∼ i/4 [57]. 6 For the symmetric operators discussed in section 4.1 and 4.2 we used the analytic results of [56] up to 8 loops, whereas for the operators studied in section 5.1 and 5.2 weak coupling data are computed by generalising the same method, based on the Pν-system, to non-symmetric sectors (see appendix C for the explicit results).  Table 1: Set of parameters used to get some of the conformal dimensions of the sl(2)-like operator with L = 1, S = 2.
Finally, before discussing the applications of the algorithm, let us focus briefly on the precision of the numerical results, which is mainly affected by the cut-off parameters: • the truncation order in the power series (3.1), N A ; • the truncation order in the asymptotic series (3.4), N a|i ; • the imaginary part of the large u approximation, N s ; • the number of sampling points, N p . Table 1 displays specific values of the cut-off parameters and the computing times corresponding to some of the results listed in table 7 of appendix D.1. The algorithm is implemented in a Mathematica notebook using a processor with 16 cores at 2.10 GHz each and 32 GB of RAM. The precision of the results was estimated by considering the number of stable digits, written in between brackets, as the truncation parameters are slightly increased.

Spectrum at finite coupling: sl(2)-like operators
The Bethe Ansatz Equations (BAEs) describing the asymptotic spectrum of the ABJM sl(2)like sector are [43] x + 4,k where x 4,j = x(u 4,j ) and σ is the BES dressing factor [7]. The BAEs (4.1) have to be supplemented by the zero momentum condition (ZMC) The states described by equations (4.1) and (4.2) correspond to single-trace operators of the form and do not form a proper closed sector of the theory, but rather a collection of states within the wider sl(2|1) sector 7 . The ABA predictions for the conformal dimension of these operators is The large-u asymptotics of the P functions, that are one of the main initial inputs of the algorithm, can be easily selected, for this sector, by setting K 1 = K 2 = K 3 = 0 and K 4 = K4 = S in (2.9). The result is [54] and besides the symmetry P 6 = P 5 can be imposed.
The first non-perturbative numerical study was performed in [53] by solving numerically the TBA [50,51] up to h = 1. The results obtained using the QSC-based algorithm described in section 3 are reported in appendix D.1, table 6. As shown in figure 4, the TBA data are consistent with our results and give an important independent test of the correctness of the method. However, the numerical precision of the data obtained in this paper with the QSC is much higher.
Since the main motivation for the study of this model resides in the weak/strong coupling AdS/CF T duality, it is particularly important to explore the large h behaviour of the spectrum. Strong coupling predictions for operators in the sl(2) sector are based on analytic continuation from the classical folded spinning string solution [40,58], and is expected to be applicable only to operators with even S. Therefore, we found no independent strong coupling predictions for the operator 20. However, the high precision results in table 6 allow us to extract a numerical prediction for the first few strong coupling coefficients of ∆ L=1,S=1 . By analogy with the S even case, we shall assume the following ansatz [19] where g = 2πh + ln 2. This is a natural parameter for the strong coupling expansion (4.8), since g ∼ (λ − 1/24)/2 at large λ, where the shift −1/24 is expected from string theory considerations [40,66]. Table 2 contains the numerical strong coupling coefficients ∆ (n) fit obtained by fitting the results in table 6 with the ansatz (4.8). Following the method of [20], the coefficients ∆ guess , corresponding to the following analytic expression To guess the fourth coefficient ∆ guess we assumed a certain similarity with the known result for S even [40,58], see also section 4.2. We stress that already the leading order coefficient in (4.9) deviates from the ones naively obtained interpolating the even-S results, confirming that this operator belongs to a different trajectory. It would be interesting to reproduce (4.9) by an analytic computation, identifying the appropriate family of classical solutions. It is important to remark here that, while the strong coupling limit of the ABA formula (4.7) gives the leading order of (4.9) matches instead the expectations of [53], and the known λ 1/4 strong coupling behaviour of analogous anomalous dimensions in N = 4 SYM [2]. In figure 5, a nice overlap of the numerical data with both (4.9) and a diagonal [6/6] Padé approximant of the weak coupling expansion up to 12 loops [56], is observed. The Padé prediction is: As noted in [53], the solution of the TBA becomes numerically unstable beyond h = 1. In the TBA setup, the unknowns are the so-called Y functions and the origin of the instability was identified by the author of [53] in the apparently divergent behaviour of ln (1 + Y 1,0 (u)) around h = 1. The aim of this section is to investigate numerically the behaviour of the function Y 1,0 for the operator 20 at strong coupling using the QSC-based algorithm. In doing so we will answer an interesting question on the behaviour of this function raised in [53]. As was noted in [53], Y 1,0 (0) gets closer and closer to the value −1 as the coupling is increased toward h 1. In view of the structure of the TBA equations, which involve convolution integrals over ln (1 + Y 1,0 (u)) , it is quite interesting to determine if Y 1,0 actually crosses the value −1 as the coupling increases further. In fact, in the presence of singularities crossing the integration contour, the TBA should be modified with the inclusion of extra residue terms [67] or, equivalently, by implementing a "desingularisation" procedure that guarantees the correct analytic continuation of the physical state [68,69]. Speculation on the existence of such critical values in the context of the TBA for N = 4 SYM was presented in [59,60].
The function Y 1,0 (u) can be reconstructed efficiently starting from the numerical solution of the QSC, which in fact allows to compute all Y functions [55]. For technical details see appendix B. It is noteworthy that, while the approach of the value −1 causes a severe instability in the TBA equations, it is harmless from the point of view of the QSC, as it simply corresponds to the zero of a Q function approaching the cut.
From the numerical outcomes displayed in figure 6, we see that Y 1,0 develops a wider and wider plateau as h is increased while (1 + Y 1,0 ) remains positive. These numerical results strongly suggest that there are no contour-crossing singularities of the kind Y 1,0 (u) = −1 for any finite value of h.
Finally, it is also interesting to investigate the analytic structure of Y 1,0 (u) as u is continued to the complex plane. Since Y 1,0 (u) is in general a complex-valued function for u ∈ C, we studied the real ratio  Since an analytic strong coupling expansion is available for all the sl(2) states with even S [40,58], it is interesting to analyse also the operator with length L = 1 and spin S = 2. In this case the BAEs (4.1)-(4.2) reduce to u 4,2 = −u 4,1 and (4.13) The numerical results obtained for the conformal dimension are reported in (4.14) In figure 8, we also plotted the strong coupling predictions obtained from the ansatz (4.8) in [40,58]:  be considered as a further strong evidence of the gauge/string duality involving the ABJM theory. A numerical prediction for the unknown coefficient at order g −5/2 is also reported: unfortunately this could not be fixed neither from the exact slope derived in [40] nor from the 1-loop results of [58].
Finally, we would like to remark that also in this case not even the leading order of (4.15) can be predicted correctly by solving the BAEs (4.13) in the strong coupling limit:  Indeed, the solution of (4.16) is which differs from the correct result (4.15) by a factor 2.
To get (4.16), only the so-called AFS phase [70] is needed as leading order contribution of the dressing factor σ: a similar computation was performed for the first time in [70] for the su(2) sector of N = 4 SYM, where also the O(1) term was correctly predicted by the leading order BAEs, as noted also in [71,72].

Spectrum at finite coupling: non-symmetric sl(2|1) operators
The QSC equations reported in section 2 allow to explore also other sectors, less studied compared to sl (2). In general, the single-trace sl(2|1) operators contain fermionic fields D + , ψ 4+ and ψ 1 † + acting on the vacuum (Y 1 Y † 4 ) L . The corresponding BAEs are obtained by setting K 1 = K 2 = K 3 = 0 in the sl(2) grading of the full Osp(2, 2|6) BAEs of [43]: Besides the intrinsic physical interest, the QSC equations for non-symmetric states exhibit novel features which are worth to be investigated, for instance the appearance of the nontrivial phase P(h) in equation (2.2). In addition, this computation can also be considered as a strong test for the consistency of the QSC in its general form. Indeed, the derivation of the gluing conditions for non symmetric operators in [55] was based on an unproven conjecture for the asymptotics of the functions τ i . Furthermore, our findings confirm that the equations of [55] form a closed system even without the exact knowledge of the state-dependent function P(h).  Table 4: Coefficients of (4.8) for the non-symmetric state with L = 2, S = 1, u 4 = u4.
as a particular case of (5.1)-(5.2). The resulting conformal dimension in the ABA approximation is For this state P 4 = −P4 = 2π 3 + O(h 2 ), therefore it is one of the simplest example where the total momentum for particles of kind 4 and4 is different from 0 or π, corresponding to a non-trivial weak coupling value for P(h): P(0) = 2π 3 . As in the previous cases, weak coupling expansions of ∆ and the P's are necessary as initial input for the iterative procedure. Since they are not available in the literature, they are computed from scratch adapting the algorithm developed in [56] and using the symmetric large-u asymptotics (4.5) for the P's but different ansatzs for the large-x expansions of P 5 and P 6 . The corresponding 8-loop perturbative results are reported in appendix C. Besides serving as initial input of the numerical algorithm, they may be considered as original findings interesting by their own: in particular, it turns out that P 6 (u) = P 5 (−u).
In contrast to the cases discussed previously, in the current case not all the P's have a definite parity in x(u). As a direct consequence of this fact a resonance problem appears when solving (2.14) in the large-u limit, due to the overlap between some of the exponents of (2.20). This problem is overcome following the strategy described in appendix A.3.
The numerical data in appendix D.2, table 8, are used to predict the strong coupling coefficients reported in table 4, and allow us to conjecture the following strong coupling expansion for the spectrum of this operator: Notice that the leading order coefficient in (5.8) can be also obtained from the large h limit of (5. showing nice agreement.

Computation of P
The generalisation of the algorithm developed in [56] to non-symmetric sectors allows us to compute analytically the first five non-trivial coefficients of P(h) = −i ln νa(0) νa(i) , (a = 1, . . . , 4) for the L = 2, S = 1 operator: In general, the exact expression for P(h) was first derived in [55]. Taking into account the non-trivial monodromy of the logarithm the complete result is , (n ∈ Z) , with F (z|k 2 ) and K(k 2 ) being the incomplete and complete elliptic integral of first kind with modulus k, respectively.
The sign ambiguity in e iP can be lifted by comparison with the leading order of the weak coupling expansion, indeed it corresponds to the ±1 ambiguity observed in the symmetric sector in [56]. In order to compute (5.13), we need to evaluate Using the Qτ -system (2.17), the quantity e F can be written in terms of the output of our numerical algorithm as where the matrix of functions f j i (u) = δ j i − τ i (u) τ j (u) can be computed as a particular combination of Q functions, as explained in appendix B. This allows us to compute P L=2,S=1 (h) non-perturbatively.
The numerical finite coupling results are displayed in figure 10 and compared with the ABA expression of the total momentum P 4 18) and the ABA approximation [55] of P(h) where Q α (u) = Kα j=1 (u−u α,j ) , (α = 4,4). As one can see in figure 10, while apparently there is no match between P L=2,S=1 (h) and P ABA 4,L=2,S=1 (h) except for h ∼ 0, interestingly the exact and ABA results converge both at small and large h. Notice also that, using the solution of (5.9), it is easy to check that P ABA L=2,S=1 tends to zero at strong coupling. In general, for any solution of the ABA (5.1)-(5.3) with Bethe roots scaling as h at strong coupling, the leading order of P ABA turns out to be quantized in integer units of π. A natural question is whether the exact formula (5.13) always reduces to (5.19) at strong coupling, and therefore if P is quantized at strong coupling for a generic state. In general, for non-symmetric sl(2|1) operators the large-u asymptotics of the P functions generalise to Differently from the symmetric case, where the following ansatz [56] for the large-x expansions of P 5 and P 6 was used, with p L (u) being a polynomial of degree L, here we use (5.23) This is the starting point to generalise the analytic algorithm of [56] to a state with L 0 = 0, and in particular to compute the 10-loop weak coupling expansions of ∆(h) and the P's needed as initial input of the program for the state with L = 4, S = 1 and L 0 = 2. The resulting expressions are reported in appendix C, with ∆(h) matching the ABA result up to h 8 .
In contrast to the previous case, here all the P functions turn out to have a definite parity in x(u), then no resonance occurs. Moreover, P L=4,S=1 (h) is found analytically to vanish up to O(h 10 ), and this behaviour is confirmed at non-perturbative level by evaluating numerically formula (5.13) up to h = 3.2. Finally, as already mentioned in section 3, in this case the inclusion of an extra gluing condition turns out to be necessary to guarantee the convergence of the algorithm.  Table 5: Strong coupling coefficients for the non-symmetric state with L = 4, S = 1, The numerical results are reported in table 9  and with the conjectured strong coupling asymptotics (see table 5): Again, the ABA predicts the same leading term as in (5.25) It is quite surprising that, for the non-symmetric sl(2|1) states discussed in this paper, the strong coupling limit of ∆ ABA matches the corresponding numerical predictions at leading order. We do not have a physical explanation for this fact. It may be just an accident, or a specific property that distinguishes between non-symmetric and symmetric operators of the sl(2|1) sector in ABJM.

Conclusions and outlook
In this paper we presented a numerical method allowing the study of the spectrum of planar ABJ(M) theory at finite coupling, in principle for any operator. The method is based on the Quantum Spectral Curve formulation obtained in [54,55] and on the numerical algorithm proposed for N = 4 SYM in [19]. Our results are first of all an important test of the QSC formulation of [54,55], which itself was based on a long chain of conjectures. 8 Besides, attached to this article we provide a simple implementation of the algorithm in Mathematica, 9 which we hope will facilitate future studies of the model. 8 In fact, this project was carried out in parallel with [55] and early numerical results were very important for developing the QSC formulation for the ABJM model. 9 For N = 4 SYM the algorithm is attached to the arXiv version of [19], see also [17].
The numerical method gives access not only to the spectrum itself but to the full set of Q functions. In the spirit of the Separation of Variables method, the Q functions may play a role not only in encoding the spectrum but also in the description of structure constants and more general observables. Encouraging results in this direction were obtained recently in [31] and [73].
An interesting generalisation of the QSC equations would be to allow for analytic continuation in the spin, similar to what done in N = 4 SYM in [21,22,40,74], where this allows to reach a BFKL regime, relevant for high-energy scattering, where the theory is similar to QCD. It would be interesting to investigate whether a similar regime exists also for the ABJM model or whether there are qualitative differences. This in turn could help reveal new properties of the spectrum and amplitudes, see [63]. Allowing for complex spin would require a modification of the algorithm presented here, in particular a change in the gluing conditions. For ABJM theory the first steps in this direction were taken in [40] and [75].
Another interesting, almost completely unexplored problem, is the study of the analytic dependence of the spectrum on the coupling constant. The spectrum has branch points in the complex domain, whose nature can be investigated efficiently numerically. We plan to report on this problem soon [57].
It would also be interesting to extend the ABJM QSC to the twisted case, in particular in view of the interest of the 3D integrable fishnet model obtained as a double scaling limit of twisted ABJM theory [76,77]. Together with the recently much studied 4D fishnet model [18,28,76,78], this non-supersymmetric model, which has an explicit Lagrangian description, allows for a direct all-loop connection between integrable spin chains and Feynman diagrams and could be very useful to develop the integrability approach for observables beyond the spectrum.
The integrable description of cusped Wilson lines is also naturally a very important open problem for ABJM theory, which could be solved adapting the QSC method as done in [29]. This in particular would allow for many strong tests of the conjecture for h(λ) due to the existence of independent localisation results (see e.g. [79][80][81] for recent results) Finally, in order to make the numerical QSC method universally applicable to any operator, it would be important to develop a systematic weak coupling algorithm covering the full Osp(2, 2|6) spectrum, in analogy to the work started in [24] where a fully automatic method of weak coupling expansion for any operator in N = 4 SYM was described. This would be very useful since the numerical algorithm requires rather precise initial data for the iterative procedure in order to converge on a given operator. It is also worth noticing that the study of the finite coupling regime with numerical methods requires extensive computational time and power. The situation worsens as the most interesting regimes -the strong coupling limit and close to branch points for complex coupling [57] -are approached. Therefore, an optimisation appears to be very desirable to pursue the numerical analysis at a more satisfactory level.

Acknowledgments
We thank Mikhail Alfimov, Lorenzo Bianchi, Sergei Frolov, David Grabner, Nikolay Gromov, Julius, Christian Marboe, Fedor Levkovich-Maslyuk, Stefano Negro and Dmytro Volin for interesting discussions and suggestions. We are especially grateful to Nikolay Gromov, Fedor Levkovich-Maslyuk and Grigory Sizov for sharing the Mathematica code for N = 4 SYM at an early stage of the project.
This project was partially supported by the INFN project SFT and the EU network GATIS+. For most of the duration of this project AC was supported by a postdoctoral fellowship from the University of Turin, and is presently supported by King's College London and the STFC grant (ST/P000258/1).

A Symmetries of QSC equations and gauge fixing
A.1 Symmetry acting on the coefficients {c A,n } n≥0 and its gauge fixing In this appendix we describe the symmetry of the QSC equations which acts directly on the coefficients {c A,n } n≥0 . It is easy to verify that the transformation where R b a is a constant unit-determinant matrix, leaves invariant the algebraic form of all the equations presented in section 1, and does not alter the analytic structure of any function involved. In addition, to preserve the asymptotics of P ab and Q a|i , R b a should have the following form where we used the fact that the charges are ordered as in (2.11), and r i ∈ R for the reality of the coefficients {c A,n } n≥0 . Choosing suitably r 2 1 , r 2 , r 5 , the transformation (A.1) can be used to fix the values of A 1 , A 2 , A 5 , 10 in (2.7): where a i are arbitrary constants. The remaining six-parameter freedom can be used to enforce the following additional gauge fixing conditions:

A.2 Residual symmetry acting on Q ij
There exists a further algebraic symmetry of the QSC equations: where G j i is a constant matrix satisfying G j i κ jk G k l = κ il . The large-u asymptotics of the QSC is preserved if G j i is lower-triangular. The condition of pure asymptotics, which is always enforced in our algorithm, breaks this symmetry almost completely, since requiring that Q a|i (u) has large-u expansion of the form (2.21) forbids a generic mixing between different columns of this matrix. However, for physical operators it is always true that the asymptotics of the second and third column of Q a|i (u) differ by an integer (see [55]): therefore, a mixing between these two particular columns does not spoil the assumption of a pure asymptotic expansion. This implies that the equations used in the numerical algorithm are invariant under a one-parameter family of symmetries, given by (A.5) with G j i taking the form: where f ∈ R is an arbitrary real parameter 12 . As we discuss below, the presence of this zero mode in our equations produces a singularity in the linear system used to compute the large-u expansion of Q a|i , due to the "resonance" between two solutions. The careful treatment of this problem is discussed in detail in appendix A.3. However, this issue does not arise for all states corresponding to P functions with definite parity in u, and additionally have integer (as opposed to only half-integer) spin. For these states, we can consistently demand that the large-u expansion (2.21) goes in even powers of 1/u only. This requirement forbids the mixing described by the transformation (A.7), so that in this case the resonance problem discussed below does not occur. 11 The particular choice of ai is irrelevant in principle. However, practically this choice is important for the convergence of the algorithm. As a rule of thumb, we should choose ai, in such a way that the coefficients {cA,n} n≥0 are roughly of the same size for any A at any fixed value of n. 12 Notice also that this transformation acts on the Q functions as follows: leaving the other Q functions unchanged. This map preserves the form of the gluing conditions, with exactly the same coefficients.

A.3 The resonance problem
Let us introduce the following notation The presence of the symmetry (A.7) implies that the computation of the coefficients B (a|i),n n≥0 is ambiguous for the column i = 2 and n ≥ N 0 . To be more explicit, let us consider the large u expansion of (2.14). At order 1/u n , we obtain a linear system of equations of the form: where M (i),n and V (i),n are respectively a 4 × 4 matrix and a 4-vector which contain the coefficients {c A,n } n≥0 and B (a|i),m m≥0 for m < n. Solving the system order by order, the matrix of coefficients appearing at order n can be fully determined in terms of the coefficients {c A,j } N 0 −1 j=1 and the charges, while the vector V n also depends linearly on the six coefficients c A,N 0 with A = 1, . . . , 6.
The concrete manifestation of the zero mode is that at the critical level we have This means that system has either zero or a one-parameter family of solutions, depending on whether a certain constraint is met on the vector of coefficients V (i=2),(n=N 0 ) . This constraint gives a linear equation to be satisfied by the coefficients {c A,N 0 } 6 A=1 . A convenient way to impose this condition is to rewrite the set of equations appearing at n = N 0 as a linear system for a different set of unknowns. Indeed, the linear system: can always be rewritten as 13 where A 0 is an arbitrarily chosen index 1 ≤ A 0 ≤ 6. Above,M (i=2),N 0 is a new 4 × 4 matrix of coefficients andV N 0 a new 4-vector. In particular,M (i=2),N 0 depends on the coefficients Choosing A 0 appropriately, the matrix of coefficients in (A.11) will have non-zero determinant, so that the system can be solved unambiguously. For the operator considered in section 5.1, the choice A 0 = 2 ensured that the system was non-singular. 14 It is also a particularly convenient choice since the coefficient c 2,n is not involved in any of the gauge fixing conditions (A.4), and it is therefore clear that the latter do not clash with equation (A.11). Notice that, as expected, we still have a one-parameter family of solutions depending on the unfixed value of B (1,2),N 0 . This coefficient is truly arbitrary, and we can set it to zero using the symmetry (A.7). This removes all residual ambiguities from the solution.
• Exclude c A 0 =2,N 0 from the set of variational parameters. This coefficient is not varied freely but instead is computed from the linear system appearing at level n = N 0 , rewritten as (A.11).

B Computation of Y 1,0 and other useful quantities
Let us briefly summarise how Y 1,0 can be computed in terms of the P and Q functions (for the context necessary to understand the following technical details see appendix A in [55]). We start from the definition of Y 1,0 in terms of the T functions where the indexes α and β label the two different wings of the T-diagram of ABJM. Then, using the Hirota equation T functions are affected by gauge ambiguities. There is a convenient gauge choice where the T functions -denoted as T -appearing in (B.1) admit a nice representation in terms of the building blocks of the Pν-system Considering that both the Y -and T -systems are naturally defined on the mirror Riemann section and ν 1 ν 4 is i-periodic on this section, T ± 1,0 can be expressed as Inverting the definition of τ in (2.16): we can write Let us now show how the functions f j i (u) can be expressed in terms of Q functions. To do this we introduce a matrix Ω j i (u) which relates Q a|i to its conjugate Q a|i (u) = Q a|i (u * ) * (see [55]): Q a|i (u) = Q a|j (u) (Ω j i (u)) + , (B.8) where Ω j i (u + i) = Ω j i (u) due to the fact that (2.14) is a real equation. As explained in [55] this matrix can be used to construct a constant gluing matrix so that the product ν 1 ν 4 can be expressed in terms of the Q functions as 1|k . (B.12) In this way, Y 1,0 is computed from the P and Q functions only. Finally, using (B.11), formula (5.17) for the integrand of the integral formula for P(h) can be recast into where the r.h.s. is expressed in terms of the Q functions computed by our numerical algorithm.

C Explicit P functions and ∆(h) at weak coupling
Here below we report the perturbative results used as input data of the numerical algorithm in the case of the sl(2|1) operator with L = 2, K 4 = K4 = 1, i.e. the analytic solutions of ∆ L=2,S=1 (h) and the related P functions up to the fourth non-trivial order at small h: Up to the fifth non-trivial order, ∆ L=4,S=1 (h) and the P-functions for the sl(2|1) operator with L = 4, K 4 = 2, K4 = 0 are given by