Quantum Spectral Curve and the Numerical Solution of the Spectral Problem in AdS5/CFT4

We developed an efficient numerical algorithm for computing the spectrum of anomalous dimensions of the planar ${\cal N}=4$ Super-Yang--Mills at finite coupling. The method is based on the Quantum Spectral Curve formalism. In contrast to Thermodynamic Bethe Ansatz, worked out only for some very special operators, this method is applicable for generic states/operators and is much faster and precise due to its Q-quadratic convergence rate. To demonstrate the method we evaluate the dimensions $\Delta$ of twist operators in $sl(2)$ sector directly for any value of the spin $S$ including non-integer values. In particular, we compute the BFKL pomeron intercept in a wide range of the 't Hooft coupling constant with up to $20$ significant figures precision, confirming two previously known from the perturbation theory orders and giving prediction for several new coefficients. Furthermore, we explore numerically a rich branch cut structure for complexified spin $S$.


Introduction
Many years of exploring integrable structures in planar N = 4 super Yang-Mills and its AdS 5 × S 5 string dual have led to a remarkably simple system of equations for the exact spectrum of the theory known as Quantum Spectral Curve (QSC) equations. They are expected to capture the conformal dimensions/string state energies at any value of the 't Hooft coupling [1,2]. The QSC equations are formulated as a set of Riemann-Hilbert type equations for a few functions. As a limiting case this system incorporates the renowned asymptotic Bethe ansatz, but also includes all wrapping corrections essential for finite length operators and perhaps constitutes the ultimate solution of the spectral problem. The Quantum Spectral Curve has a transparent algebraic origin related to the underlying psu(2, 2|4) symmetry of the problem. It was derived in [1,2] from the discrete integrability of Hirota dynamics underlying the Y-system enhanced with a specific analyticity condition worked out in [3,4,5,6,7]. In contrast to some explicitly known integral forms of the Y-system [5,8,9,10,11] the QSC is applicable for any local operator or string state.
The striking simplicity of the QSC formulation has already led to a rapidly growing body of exact as well as perturbative results. First, in [1,12] it was pointed out that it can be solved analytically in a near-BPS regime, e.g. for expansion in the coupling λ or in the spin S, in particular in [12] the anomalous dimensions for twist operators in the sl(2) sector were computed at any coupling to order S 2 , giving new analytical predictions at strong coupling for some local operators and the BFKL pomeron intercept. As an extension of the observation of [1] powerful techniques for expansion in the coupling were developed in [13,14] which led, impressively, to 10-loop anomalous dimensions for numerous operators, as well as 6-loop predictions at any S for twist-2 operators. Finally, in [15] the analytic continuation to the BFKL region S ∼ −1 was explored, and the leading order BFKL equation in SYM was derived in this approach. It is also important to mention that the QSC method was also applied to ABJM theory [16] opening a way for various explicit calculations in particular to that of the interpolation function h(λ) [17], the mysterious extra component of the spectral problem in this model.
Even though many explicit analytic results are now available both at strong and weak coupling, one important range of applications of the QSC that has remained unexplored till now is the numerical investigation of the spectrum at finite coupling.
Previous numerical methods based on TBA, even limited to a few operators 1 , low precision and slow convergence rate gave, nevertheless, several highly important results, allowing, in particular, the computation of the anomalous dimension of a nonprotected (Konishi) operator in a planar 4d theory at finite coupling [18]. Numerics also gave a prediction for the strong coupling Konishi anomalous dimension which was later confirmed by several methods [19,20,21,22,12,23,24,25] . The main goal of the present work is to remove the limitations of the previously known methods by developing an algorithm for a numerical solution of the QSC.
The low precision and performance of the TBA-like approach was mainly due to the complicated infinite system of equations and cumbersome integration kernels. The QSC includes only a few unknown functions and thus can be expected to give highly precise numerical results. However, the QSC equations are functional equations supplemented with some analyticity constraints of a novel type which makes it a priori not a trivial task to develop a robust numerical approach.
In this paper we propose an efficient method to solve the QSC numerically and illustrate our method by a few examples. Among the several equivalent formulations of the QSC we identified the equations which are best-suited for numerical solution 2 . We implemented our algorithm in Mathematica and were able to get a massive increase in efficiency compared to the TBA or FiNLIE systems [18,8,24,9]. With one iteration taking about 2 seconds we only need 2 − 3 iterations (depending on the starting points) to reach at least 10 digits of precision. Quite expectedly, the precision gets lost for very large values of the 't Hooft coupling. Nevertheless, without any extra effort we reached λ ∼ 1000 keeping a good precision, which should be more than enough for any practical goal.
Not only does our approach work for any finite length single trace operator and in particular for any value of the spin, it also works with minimal changes even away from integer quantum numbers! We demonstrate this in the particularly interesting case of the sl(2) twist-2 operators. Their anomalous dimension analytically continued to complex 1 only for a few operators the complicated structure of the "driving terms" was deduced explicitly in a closed form. Even for those operators the driving terms may change depending on the value of the coupling.
2 one may call this sub-system of equations as Pω-system, in contrast to previously used Pµ-system or values of the spin S is known to have a very rich structure, in particular the region S −1 is described by BFKL physics. As we show, within the framework of QSC it is not hard to specify any value of the Lorentz spin S as the conserved charges enter the equations through the asymptotics which can in principle take any complex values. Then we can compute the analytically continued scaling dimension ∆ directly for complex S (or even interchange their roles and study S as a function of ∆). The result of this calculation can be seen on Fig. 1. Let us stress that the algorithm is very simple and mainly consists of elementary matrix operations. As such it can be easily implemented on various platforms. In particular, we believe the performance could be increased by a few orders with a lower level, e.g. C++, implementation. In this paper we mostly aim to demonstrate our algorithm, prototyped in Mathematica.
The structure of the paper is as follows. In section 2 we introduce the Quantum Spectral Curve and discuss the key equations for the numerical implementation. In section 3 we describe how our algorithm works in general. We also demonstrate our method on specific examples for twist-2 operators in the sl(2) sector, starting with physical operators (such as the well studied Konishi operator). Then in section 4 we explain how to extend u u Figure 2: P a and Q i have one cut on the real axis in the representations with short and long cuts respectively. The ellipse shows the region of convergence of the series (2.4) the method to non-integer values of the spin. We discuss Fig. 1 in a bit more detail, and also describe our high precision evaluation of the BFKL pomeron intercept 3 . The final section contains our conclusions and speculations.

Review of the AdS/CFT Quantum Spectral Curve
Let us start by introducing the main equations of the Quantum Spectral Curve framework. For a full review of the QSC we refer the reader to [2], but we will try to make the paper self-contained.
For integrable spin chains Baxter Q-polynomials play a central role. They carry complete information about the state and the spectrum and can also be thought of as wave functions of the model. They are completely determined by a set of functional relations (known as Q-system) and by the polynomiality condition.
In N = 4 SYM the situation is similar. Its spectrum of anomalous dimensions is described by a set of Q-functions satisfying the same Q-system relations. However, in this model the Q-functions are no longer polynomials, but are analytic functions of the spectral parameter u with branch cuts. The positions of the branch points depend on the 't Hooft coupling λ: they are situated at ±2g + in, n ∈ Z, where g = √ λ 4π . In the limit when λ is small the system reduces to the usual psu(2, 2|4) Heisenberg spin chain. For the class of functions with cuts the Q-system alone is not constraining enough and one has to specify the monodromies of the Q-functions in order to close the system. These monodromies were found in [1] and take the form of a Riemann-Hilbert problem as described below. The resulting set of equations found in [1,2] is known under the name Quantum Spectral Curve (QSC). The Q-system of N = 4 SYM is composed of 2 8 Q-functions. However, the algebraic relations between them allow us to choose a much smaller subset, which will be complete in the sense that the rest of Q-functions can be generated from the selected ones algebraically. A convenient choice for such a subset consists of 4 + 4 functions P a (u) and Q i (u) (a, i = 1, . . . , 4). One can say that P a describe the S 5 degrees of freedom whereas Q i correspond to the AdS 5 part. A particularly nice property of P's is that they have only two branch points at ±2g when they are connected by a "short" cut [−2g; 2g] (see Fig. 2). This means that there are more branch points on the next sheet, but for this choice of the cut they do not appear on the first sheet. Very similarly Q's have only two branch point on the main sheet if the cut is taken to go through infinity. In a sense this reflects the non-compactness of the AdS 5 part of the space.
Whereas the coupling determines the position of the branch points, the quantum numbers of the state are specified through the large u asymptotics of Q-functions. P a encode the compact bosonic subgroup SO(6) quantum numbers (J 1 , J 2 , J 3 ), while Q i gives the SO(4, 2) charges (∆, S 1 , S 2 ), which include the conformal dimension of the state ∆. Explicitly Note that for the numerical implementation P a are more handy, in the sense that they can be expressed as a series in the Zhukovsky variable x(u) defined by u = g(x + 1/x), .
This series is convergent everywhere on the upper sheet and also in an elliptic region around the cut on the next sheet (see Fig. 2). A similar parametrization for Q a will not cover even the upper sheet. Fortunately, in the whole set of 2 8 Q-functions there are other 4 functions with one single cut, which are denoted as P a (u), a = 1, . . . , 4. Together with P a (u) they also form a complete set of Q-functions. In particular, one can reconstruct Q i from them. The procedure for this, which will be crucially important in our numerical implementation, is the following: • Find a set of 16 functions Q a|i , satisfying Note that this is a 4-th order finite difference equation, which entangles all Q a|i with fixed i. Different values of i label the 4 linearly independent solutions of this equation. One could also equivalently use Q b|i (u − i 2 ) in place of Q b|i (u + i 2 ) in the r.h.s., due to the constraint [2] P a P a = 0 . (2.6) This constraint also fixes some of the coefficients c a,n .
• The matrix Q a|i can then be used to pass to Q i from P a 's, The equations (2.5) and (2.7) are simply two of the Q-system relations as explained in [2]. We also introduce a matrix Q a|i such that Q a|i Q a|j = −δ i j and use it to define Q's with an upper index: Note that since Q a|i (u) is analytic in the upper-half-plane we can also analytically continue these relations around the branch point at u = 2g to get where the tilde denotes analytic continuation to the next sheet.
Since Q i can now be recovered from P a and P a it is not surprising that actually all information we needed, in particular all the charges (including those in AdS 5 ), are encoded in P's alone, through whereM j andM a are defined in (2.2), (2.3) and there is no summation over a 0 in l.h.s. In particular, one can extract ∆ from the last equation.
The coefficients c a,n and corresponding coefficients c a,n of the expansion of P a (u) need to be found. The constraint (2.6) fixes some of them (for example, we can use it to fix all c 1,n ). The condition (2.11) gives the leading coefficients c a,Ma . The remaining coefficients should be fixed from the analyticity constraints on P's as prescribed by QSC. Let us describe these constraints. The analytic continuation of P a to the second sheet, which we denote byP a , in terms of our ansatz (2.4) becomes simplỹ (2.12) According to [2] we should havẽ where µ ab (u) is an antisymmetric matrix with unit Pfaffian, i-periodic as a function with long cuts, with the discontinuity fixed in terms of P ã (2.14) Knowing the r.h.s. of (2.14) it is straightforward to reconstruct µ ab itself using the spectral representation of µ ab In theory one could reconstruct µ ab for some fixed coefficients c a,n from (2.14) and then impose (2.13) to fix the unknown c a,n . However, in practice that is very hard to do as the power series (2.4) is only convergent inside the region shown on Fig 2. The problem is that this region does not cover the entire cut of µ, which stretches to infinity. In other words, it would be very hard to reconstruct µ ab from some given coefficients c a,n in this direct way. We found that it is much more advantageous to close the analyticity conditions at the level of Q i , which obey very similar equations. The rule is quite simple -one has to interchange short and long cuts. That is, we have to introduce an i-periodic with short cuts function ω ij (u) such that (2. 16) with Pf ω = 1. Now we see that to recover ω ij one only needs to know its discontinuity on the interval [−2g, 2g], which is completely inside the region of convergence on Fig. 2! Thus the discontinuity of ω can be expressed in terms of the coefficients c a,n via (2.7) and (2.9), provided we know how to solve (2.5) for arbitrary P a and P a . In the next section we will describe an algorithm which allows to solve (2.5) very efficiently and then find the coefficients c a,n , which yields the solution of the QSC.

Step 1: Solving the equation for Q a|i
As we explained above the quantity Q a|i is at the heart of our procedure. In this section we will demonstrate how this set of 16 functions can be found for arbitrary P a and P a .
In this procedure the precise ansatz for P is not important. However, as we will see later, we should be able to compute P a (u)P b (u) on the upper sheet for u with large imaginary part. Of course, having the ansatz in the form of a (truncated) series expansion (2.4) we can easily evaluate it everywhere on the upper sheet numerically very fast. The process of finding Q a|i is divided into two parts. Firstly, we find a good approximation for Q a|i at some u with large imaginary part (in the examples we will need Im u ∼ 10 − 100). At the next step we apply to this large u approximation of Q a|i a recursive procedure which produces Q a|i at u ∼ 1.
Large u approximation. For Im u ∼ 10 − 100 we can build the solution of (2.5) as a 1/u expansion. This is done by plugging the (asymptotic) series expansion of Q a|i into (2.5): where N is some cutoff (usually ∼ 10 − 20). This produces a simple linear problem for the coefficients B a|i,n , which can be even solved analytically to a rather high order. The leading order coefficients of Q a|i can be chosen arbitrarily. After that the linear system of equations becomes non-homogenous and gives a unique solution in a generic case. 4 Finite u approximation. Once we have a good approximation at large u we can simply use the equation (2.5) to recursively decrease u. Indeed defining a 4 × 4 matrix Iterating this equation we get, in matrix notation For large enough N we can use the large u approximation (3.1) for Q b|i in the r.h.s. As a result we obtain the functions Q a|i for finite u with high precision.

Step 2: Recovering ω ab
Now when we have a good numerical approximation for Q a|i (u) we can compute Q i and Q i from (2.7) and (2.9) and plug them into the discontinuity of ω ij (2.14). After that we can recover ω ij from its discontinuity modulo an analytic function using its spectral representation where the "zero mode" ω 0 ij (u) is the analytic part of ω ij -it has to be periodic, antisymmetric in i, j and should not have cuts. We will fix it below. We note that we only need to know values of Q andQ on the cut. In our implementation we use a finite number of sampling points on the cut given by zeros of Chebyshev polynomials. One can then fit the values ofQ i Q j − Q iQj at those points with a polynomial times the square root u 2 − 4g 2 . After that we can use precomputed integrals of the form (3.5) with high precision by a simple matrix multiplication, which produces the result at the sampling points u A in an instant.
One more point to mention here is that in our implementation we only compute ω reg ij = 1 2 (ω ij −ω ij ) at the sampling points to avoid the problem of dealing with the singularity of the integration kernel. Note that ω reg ij can be used instead of ω ij in the equations like (2.16), because the difference is proportional to Q i Q i which is zero similarly to (2.6), as can be shown by combining (2.6) with (2.7), (2.8).
Finding zero modes. It only remains to fix ω 0 ij (u). First we notice that for all physical operators ω ij should not grow faster than constant at infinity [2]. As the integral part of (3.5) does not grow either and since ω 0 ij (u) is i-periodic it can only be a constant. To fix this constant we use the following observation [2]: the constant matrix α + ij which ω ij approaches at u → +∞ and the constant matrix α − ij which it reaches at u → −∞ are restricted by the quantum numbers [2]. To see this we can pick some point on the real axis far away from the origin and shift it slightly up into the complex plane, then from (2.16) we have Next, notice that since Q j is analytic everywhere except the cut on the real axis, it can be replaced by its asymptotics above the real axis, i.e. Q j (u + i0) ∼ B j u −M j , and also Q j (−u + i0) ∼ B j u −M j e −iπM j , as we find from the previous expression by a rotation by π in the complex plane. As a result we get the asymptotics of Q i at infinities and slightly below the real axis Using that they are related by the analytic continuation in the lower half plane the first equation also gives Combining this with (3.8) we get a relation between the constant asymptotics of ω at the two infinities At the same time from (3.5) we get which implies that We see that the zero modes can be also computed from the values of Q andQ on the cut. Note also that the r.h.s. is not explicitly antisymmetric. Imposing the antisymmetry gives I kl (cot πM l − cot πM k ) = 0, (3.13) so either I kl = 0 or cot πM l = cot πM k . As Pf ω = 1, all I kl can not be equal to zero simultaneously. Having I kl non-zero implies quantization of charges: for example, the choice I 12 = 0 and I 34 = 0, which is consistent with perturbative data, requires cot πM 1 = cot πM 2 and cot πM 3 = cot πM 4 , and so S 1 , S 2 have to be integer or half integer. In section 4 we will see how to relax this condition and do an analytic continuation in the spin S 1 to the whole complex plane.

Step 3: Reducing to an optimization problem
Having ω ij and Q a|i at hand we can try to impose the remaining equations of the QSC (2.16). We notice that there are two different ways of computingQ i , which should give the same result when we have a true solution: (2.9) and (2.16). Their difference, which at the end should be zero, is 14) The problem is now to find c a,n for which F i (u) is as close as possible to zero. Here we have some freedom in how to measure its deviation from zero, but in our implementation we use the sum of squares of F i at the sampling points u A . Then the problem reduces to the classical optimization problem of the least squares type. In our implementation we found it to be particular efficient to use the Levenberg-Marquardt algorithm (LMA), which we briefly describe in the next section. The LMA is known to have a Q-quadratic convergence rate, which means that the error n decreases with the iteration number n as fast as e −c 2 n . The convergence is indeed so fast that normally it is enough to do 2 or 3 iterations to get the result with 10 digits precision. We give some examples in the next section.
Levenberg-Marquardt algorithm Our problem can be reformulated as follows: given a vector function f i (c a ) of a set of variables c a (which we can always assume to be real) find the configuration which minimizes Assuming we are close to the minimum we can approximate f i by a linear function: which gives the following approximation for S(c a ): The approximate position of the minimum is then at ∂c a S = 0 for which we get from which, in matrix notation, We see that for this method we should also know the derivatives of our F a (u) w.r.t. the parameters c a,n , which in our implementation we find numerically by shifting a bit the corresponding parameter. In some cases, when the starting points are far away from the minimum, the above procedure may start to diverge. In such cases one can switch to a slower, but more stable, gradient descent method for a few steps. The Levenberg-Marquardt algorithm provides a nice way to interpolate between the two algorithms by inserting a positive parameter Λ into the above procedure, c n+1 = c n − (JJ T +JJ T + ΛI) −1 (Jf + Jf ) . (3.20) The point is that for large Λ this is equivalent to the gradient descent method. Thus one can try to increase Λ from its zero value until S(c n+1 ) < S(c n ) and only then accept the new value c n+1 . This helps a lot to ensure stable convergence.
In the next section we demonstrate the performance of our numerical method by applying it to the twist-2 operators in sl(2) sector.

Implementation for the sl(2) Sector and Comparison with Existing Data
The sl(2) sector in the QSC framework Although our method can be used for any state of the N = 4 SYM theory, the examples we provide in this paper are for states in the sl(2) subsector. In this section we will discuss the physical operators which have integer spin, and demonstrate our numerical method in action for the Konishi operator. Then in section 4 we will show how the algorithm works for other states with S no longer an integer.
Let us sketch the neeeded information about the QSC in the sl(2) subsector (more details can be found in [2,1,15]). Operators in this sector have only three non-zero quantum numbers: spin S ≡ S 1 , twist L ≡ J 1 and conformal dimension ∆. Twist-L single-trace operators of this subsector can be schematically represented as where D is a light-cone derivative, Z is a scalar of the theory and the dots stand for permutations. For such states the QSC enjoys several simplifications. First, quantities with upper and lower indices are now related to each other: indices can be raised or lowered using a simple matrix 5 for example, It is also easy to show that in this sector ω ij should satisfy ω 14 = ω 23 in addition to antisymmetry. Second, Q-functions have now definite parity in u, which decreases the number of expansion coefficients in series like (3.1) two times. Finally, one can write down a simplified version of asymptotics (2.1): Coefficients A a and B i are related to the global charges L, S, ∆ (see [1,2]).

Implementation for Konishi
Here we discuss the convergence on a particular example of the Konishi operator which corresponds to S = 2, L = 2. The reason we start from this operator is because it is very well studied both analytically at weak and strong coupling and also numerically. So we will have lots of data to compare with. To start the iteration process described in the previous sections, we need some reasonably good starting points for the coefficients c a,n . For the iterative methods, like, for instance, Newton's method, good starting points are normally very important. Depending on them the procedure may converge very slowly or even diverge. We made a rather radical test of the convergence of our method by setting all coefficients to zero, except the leading ones, which are fixed by the charges. For ∆ we took the initial value 4.1 at the value of 't Hooft coupling g = 0.2. To our great surprise it took only 12 steps to converge from the huge value of S(c a ) ∼ 10 +7 (defined in (3.15)) to S(c a ) ∼ 10 −9 . The whole process took about 20 seconds on a usual laptop (see Fig. 4), producing the value ∆ = 4.4188599, consistent with the best TBA estimates [18,24]. After that we used the obtained coefficients as starting points for other values of the coupling to produce the Table 1. All the values obtained are consistent with the TBA results within the precision of the latter, being considerably more accurate at the same time.
The reason for such an excel-  In the next section we discuss the analytic continuation in S away from its integer values. This is an important calculation which bring us to a highly accurate numerical estimate for the pomeron intercept at finite coupling -a quantity which can be studied exclusively by our methods.

Extension to Non-Integer Lorentz Spin
In this section we explain which modifications are needed in order to extend our method to non-integer values of spin S, and give two specific examples of calculations for such S.

Modification of the Algorithm for Non-Integer Spin
First we need to discuss how the procedure of fixing zero modes of ω's described in section 3.2 is modified for non-integer S. The main difference stems from the fact that analytic continuation to non-integer S changes the asymptotic behavior of ω ij at large u, as described in [12,15,27]. While for integer S asymptotics of ω are constant, for non-integer S some components of ω have to grow exponentially. Without this modification the system has no solution: indeed, in section 3.3 we assumed constant asymptotics of all ω's and derived quantization condition for global charges. A minimal modification would be to allow exponential asymptotics in one of the components of ω. In order to understand which of the components can it be, let us recall the Pfaffian constraint satisfied by ω ij Pf ω = ω 12 ω 34 − ω 13 ω 24 + ω 2 14 = 1.
First, it is clear ω 14 alone can not have exponential asymptotics. Second, in the case of integer S both ω 12 of ω 34 are non-zero constants at infinity [15,12]; then shifting S infinitesimally away from an integer we see that it would be impossible to satisfy the condition (4.1) if we allow one of them to have exponential asymptotics at infinity: this exponent will multiply the constant in the other one. So the only two possibilities left are ω 13 and ω 24 , which are both zeros at infinity for integer S. From perturbative data we know that it is ω 24 which should have exponential asymptotics. Thus we formulate the "minimal" prescription for analytic continuation of Q-system to non-integer S: e 2π|u| asymptotic has to be allowed in ω 24 . This prescription was tested thoroughly on a variety of examples [28,29,15,12,27], but it would be interesting to derive it rigorously and generalize it to states outside of the sl(2) sector. Of course, one can also consider adding exponents to more than one component of ω ij : in this case the solution will not be unique. A complete classification of solutions of Q-system according to exponents in their asymptotics might be interesting. For example it is known that allowing for an exponent in some other components corresponds to the generalized cusp anomalous dimension [12,30]. Because of the exponential asymptotics of ω 24 , the argument in section 3.2, which fixes the zero modes of ω, has to be modified. First, formula (3.12) still holds true for i = 1 or i = 3, as ω 24 does not enter anywhere in the derivation. Thus Consequently, one can use (3.12) for both ω 13 and ω 31 , and reproduce the quantization condition (3.13) for global charges, which in this case implies that either ∆ = 0 or ω 13 = 0. Equation (3.12) can also be used for ω 14 and ω 23 (which are equal) and imposes that either ∆ = 0 or ω 14 = 0. It remains to fix the zero mode in ω 0 24 , for which we use an ansatz The constants a 1 , a 2 , a 3 can be found from the Pfaffian constraint (4.1). To this end we expand the hyperbolic cotangent in (3.5) to get where the terms of the expansion are integrals similar to I ij with additional factors of e 2πu or e 4πu in the integrand 6 . Analogous expansion can be obtained at u → −∞. Then plugging these expansions into (4.1) we get formulas for the coefficients a 1 , a 2 , a 3 . For example, With these modifications we can reconstruct all ω ij including the zero modes and then proceed with our algorithm as in the case of integer S.

Exploring Complex Spin
In this section we briefly describe the results of our numerical exploration of ∆(S) as an analytic function of a complexified spin S. As explained in the previous section the generalization of our numerical method to arbitrary values of spin requires minimal modifications of our main code. Thus we are able to generate numerous values of the anomalous dimension for any S with high precision in seconds. In fact both S and ∆ enter the QSC formalism on almost equal footing and we can also switch quite easily to finding S for given ∆. This is what is adopted in the vast literature on the subject and what we are going to consider below. This viewpoint is particular convenient due to the symmetry ∆ → −∆ which makes the pictures particularly nice.
Starting from S = 2 (Konishi operator) we decreased the value of S or ∆ in small steps using the solution at the previous step as a starting point for the next value. In this way we built the upper two curves on Fig. 6. Let us point out one curious technical problem -one can see for instance from (4.5) that the lines S = ±∆ + Z are potentially dangerous due to the divergence. In fact we found that near these dangerous points on the line the factor I 12 I 34 also vanishes canceling the potential divergence. This however affect the convergence "radius" of our iterative procedure and we found it quite complicated to cross those lines, even though in very small steps we were able to reach close to them. The way out is to go around these lines in the complex plane ∆. To make sure there is no true singularity or branch point we also explored a big patch of the complex plane ∆, indeed finding some branch points, but deep inside the complex plane, having nothing to do with these lines. For example when g = 0.2 we found 4 closest branch points at roughly ±1 ± i, see Fig. 1. By making an analytic continuation (described above) through those cuts we found another sheet of the Riemann surface S(∆). On this sheet we have found four cuts: two are connecting it to the first sheet and two other ones, located symmetrically on the imaginary axis, lead to further sheets. We expect an infinite set of sheets hidden below and also more cuts on both sheets outside of the area that we have explored.
It is instructive to see how this Riemann surface behaves as g → 0. First, the real parts of branch points on the physical sheet are very close to ±1, but the imaginary part goes to zero. Thus at infinitely small g the cuts collide, isolating the region |Re ∆| < 1  Figure 7: The BFKL intercept j as a function of coupling λ. The red solid line with tiny red dots is obtained by our numerical procedure. It interpolates perfectly between the known perturbative predictions (the blue dashed lines) at weak [31,32] and strong coupling [33,34,35,12].
from the rest of the complex plane. These two separated regions become then the areas of applicability of two different approximations: for |Re ∆| > 1 one can apply the usual perturbation theory and Beisert-Eden-Staudacher Asymptotic Bethe Ansatz, whereas the region |Re ∆| < 1 is described by BFKL approximation and so-called Asymptotic BFKL Ansatz [27]. The presence of the cut can be to some extent deduced from perturbative perspective in each region: in the regime of usual perturbation theory ∆(S) = 2 + S − 8g 2 H S + O(g 4 ), (4.6) where H S is the harmonic number. It has poles for all negative integer values of Sthese poles are weak-coupling remnants of the cuts we see at finite coupling. In the BFKL regime one should instead look at the leading order BFKL equation [36,37,31] To make sense of this equation one has to take the limit g → 0, S → −1 so that the l.h.s stays finite. Then the ψ-functions in the r.h.s generate poles at odd values of ∆, which, again, are cuts degenerated at weak coupling. Fig 6 represents a section of the Riemann surface by the plane Im u = 0, i.e. dependence of S on ∆ for real ∆, which, of course, consists of two curves, originating from the two sheets we explored. At weak coupling the upper curve becomes piecewise linear, approaching different parts of the dotted line: for |∆| > 1 it coincides with S = ±∆ − 2 and for |∆| < 1 it becomes S = −1. One could expect a similar piecewise linear behavior for the lower curve: it approaches S = ±∆ − 2 for |∆| < 1, approaches S = 0 in some region outside of |∆| < 1 and becomes a certain linear function even further away from ∆ = 0. It would be interesting to explore the complete analytic structure of this Riemann surface, and understand what describes its asymptotics when g → 0. It should produce a hierarchy of "Asymptotic Bethe Ansätze" each responsible for its own linear part of the limiting surface.

BFKL Pomeron Intercept
The pomeron intercept j(λ) is a quantity which relates spectrum of single-trace operators and scaling of high energy scattering amplitudes in the Regge regime [38]. This regime is particularly interesting, since it establishes a connection between results in N = 4 SYM and multicolor QCD: the non-trivial leading order of so-called BFKL eigenvalue is the same in two theories, and in the higher orders N = 4 SYM is expected to reproduce at least the most complicated part of the QCD result.
Our goal is to demonstrate the universal power of our approach by giving an extremely precise numerical estimate for this important quantity at finite coupling in a wide range of couplings.
One defines the intercept as j = S(∆ = 0) + 2, where S is the spin of the twist-2 operator such that ∆(S) = 0. Having formulated the problem like this, we can in principle apply the algorithm described in section 3 to find the correct value of S, while keeping ∆ at zero. However, one may already suspect that the point ∆ = 0 is special. Indeed, we know that for any solution of QSC there is always another one related by ∆ → −∆ symmetry. At the level of Q i functions this allows simultaneously interchanging Q 1 ↔ Q 3 and Q 2 ↔ Q 4 as one can see from the asymptotics (3.25). From this we see that at small ∆ two different solutions of QSC (related by the symmetry) approach each other, making the convergence slow, exactly like Newton's method becomes inefficient for degenerate zeros. In other words, in the limit ∆ → 0 the Q's related by the symmetry become linearly dependent in the leading order. Furthermore, since the matrix Q a|i should stay invertible, the leading coefficients B i of asymptotic expansion of Q i diverge at ∆ → 0.
The way out is to perform a linear transformation of Q's preserving the equations: it will replace two of them by linear combinations Q 3 − γQ 1 and Q 4 + γQ 2 with some coefficient γ, so that the divergent leading order cancels and the four functions Q i become linearly independent.
For the gauge choice in which B 1 = B 2 = 1 the transformation acts on i-indices of Q-functions with a matrix 7 .
One can check that rotation by this matrix will render Q a|i finite and linearly independent, and moreover, preserve relations (3.23). After this one can apply the standard procedure from section 3 with the only modification that the large u expansion of Q a|i will contain log u/u n terms in addition to the usual 1/u n .
Having done this, we can readily generate lots of numerical results. In particular we built numerically the function j(λ) which interpolates perfectly between the weak and strong coupling predictions. We have found j(λ) with high precision (up to 20 digits) for a wide range of 't Hooft coupling (going up to λ 1000). The results are also summarized in the Table 2.
1.000 000 000 000 000 000 0 0.  Table 2: Numerical data for the pomeron intercept for various values of the 't Hooft coupling. All digits are expected to be significant but some additional tests are in progress, and will be reported in second version of this preprint. Table 2 represents a small portion of all data we generated, which is available by request. In particular we generated ∼ 100 points with small g in the range 0.1 . . . 0.017 each with more than 20 digits precision. Fitting this data with powers of g 2 we found j = 1 + 11.09035488895912g 2 − 84.0785668075g 4 − 2543.0481652g 6 + 156244.8086g 8 where the first 3 terms are known analytically from Feynman diagram perturbation theory calculations [31,32] and their numerical values coincide in all digits with our prediction above. The last two terms give our numerical prediction for the numerical values of the NNLO and NNNLO BFKL pomeron intercept. Our fit also gives predictions for the higher corrections but with a smaller precision.

Conclusions and Future Directions
In this paper we have demonstrated that in addition to their analytic power, the QSC equations can give highly precise numerical results at finite coupling. We develop a numerical procedure which applies to generic single trace operators and as such it is unique in its kind. Furthermore, the algorithm converges at a remarkably high rate which gives us access to high numerical precision results -up to 20 digits or even more in a few iterations.
The efficiency of our method is demonstrated on the example of sl(2) sector operators. We also formulated how to extend our procedure to non-integer quantum numbers. We studied the twist-2 operators for complex values of the spin discovering a fascinating Riemann surface (see Fig.1). In addition we reformulated our equations to be directly applicable to the BFKL pomeron intercept and evaluated the intercept j with high precision of up to 20 significant figures. By fitting our data we also gave a prediction for the perturbation theory expansion j(λ) = 1 + 0.07023049277268284 λ − 0.00337167607361 λ 2 (5.1) − 0.00064579607573 λ 3 + 0.0002512619258 λ 4 + . . .
reproducing correctly the first two nontrivial orders [31,32] and giving a prediction for higher orders. The range of possible applications of our method is vast. First, it is not limited solely to the sl(2) sector of N = 4 SYM, but is directly applicable to any single trace operators of the theory. It would be interesting to do an explicit example of a numerical calculation with our algorithm outside of the sl(2) sector. For example, the wider class of sl(2, C) operators (identified in [27]), also exhibiting a BFKL regime, could be a good candidate to begin with. Second, we expect our method to be applicable for such non-local operators as the generalized cusp anomalous dimension and quark-anti-quark potential, DD-brane and other boundary problems [39,40,11,41,30]. Third, it may be interesting to generalize our method to ABJM theory as well as to various integrable deformations of N = 4 SYM theory.
The numerical results could also be helpful for the analytical exploration of the spectrum -for instance, in such regimes as BFKL and at strong coupling, which remains almost unexplored, and various limiting cases of the generalized cusp. Furthermore, studying numerical results and the behavior of the generated Q-functions in various limits can reveal new analytically solvable regimes.
To facilitate further applications and development of our method we are planning to release to the public a user-friendly version of our code as a Mathematica package. It may be also useful to convert our Mathematica code into a more decent and fast C++. It should not be difficult as our algorithm is quite simple and only uses some basic matrix operations.