Strong coupling results from the numerical solution of the quantum spectral curve

In this paper, we solved numerically the Quantum Spectral Curve (QSC) equations corresponding to some twist-2 single trace operators with even spin from the $sl(2)$ sector of $AdS_5/CFT_4$ correspondence. We describe all technical details of the numerical method which are necessary to implement it in C++ language. In the $S=2,4,6,8$ cases, our numerical results confirm the analytical results, known in the literature for the first 4 coefficients of the strong coupling expansion for the anomalous dimensions of twist-2 operators. In the case of the Konishi operator, due to the high precision of the numerical data we could give numerical predictions to the values of two further coefficients, as well. The strong coupling behaviour of the coefficients $c_{a,n}$ in the power series representation of the ${\bf P}_{\!a}$-functions is also investigated. Based on our numerical data, in the regime, where the index of the coefficients is much smaller than $\lambda^{1/4}$, we conjecture that the coefficients have polynomial index dependence at strong coupling. This allows one to propose a strong coupling series representation for the ${\bf P}$-functions being valid far enough from the real short cut. In the paper the qualitative strong coupling behaviour of the ${\bf P}$-functions at the branch points is also discussed.


Introduction
Maldacena's famous AdS/CFT correspondence [1,2,3] is the best elaborated holographic duality conjecture between gauge and string theories. The discovery of integrability on both sides of the correspondence [4], created a hope to find the exact solution of the theory in the planar limit. The mathematical apparatus offered by integrability, proved to be the most efficient in computing the planar spectrum of anomalous dimensions/string energies. In the large volume limit the spectrum 1 was described by the Asymptotic Bethe Ansatz (ABA) equations [5] which account for all power-like corrections in volume, but neglects the exponentially small wrapping corrections. The wrapping corrections [6] were taken into account by the so-called Lüscher-formulae [7,8,9,10,11] which are now available up to the second order in wrapping [12,13]. The Thermodynamic Bethe Ansatz (TBA) technique was the first method which could sum up all wrapping corrections to the ABA in the form of a set of infinite component nonlinear integral equations [14,15,16,17,18,19]. Though the TBA equations could provide important results, both in the weak [20,21,22] and in the strong [23,24,25,26] coupling regimes 2 , its analytical and numerical treatment proved to be tedious, due to the cumbersome kernels and the infinite number of unknown functions. Later the FiNLIE method [27], which can be considered as an improved finite version of the TBA, allowed one to reach better results in the perturbative regime [28,29], but the structure of the equations was still so complicated that it required reasonable human effort to reach higher and higher orders in the perturbative regime. the perturbative regime. In [32,33] even 10-loop analytical results were obtained for some operators in the sl(2) sector. QSC was powerful to get analytical results also in the near-BPS regimes [30,34]. In [34] analytical next-to leading order results were obtained in the small spin expansion for the anomalous dimensions of twist operators in the sl(2) sector, providing also analytical predictions for the strong coupling expansion coefficients of the anomalous dimensions for some local operators and for the BFKL pomeron intercept. In [35] leading order BFKL equation was derived by performing the S → −1 analytical continuation.
Later, in [36] an efficient numerical algorithm was proposed for solving the Pµsystem and it was used to confirm 2 previously known and to predict several previously unknown coefficients in the weak coupling expansion of the BFKL pomeron intercept.
Recently, analytical expression was obtained for the next-to-next-to leading order of the BFKL pomeron eigenvalue in [37], and the QSC description of cusped Wilsonlines [38] and of the quark-anti-quark potential [39] were worked out.
In this paper we consider twist-2 operators with even positive integer spin. Using the numerical method of [36], we perform the numerical solution of the Pµ-system for the twist-2 states with S = 2, 4, 6, 8 in a wide range of the t'Hooft coupling.
Though analytical strong coupling results are available in the literature for the anomalous dimensions of the states under consideration, they come from small spin results matched with classical and quasi classical string-theory results [34] and not directly from the strong coupling solution of the Pµ-system. This is why the aim of the paper is to gain a deeper insight into the strong coupling behaviour of the solutions of the Pµ-system.
In the S = 2, 4, 6, 8 cases, our accurate numerical results confirmed the analytical predictions of [34] for the first 4 coefficients of the strong coupling expansion for ∆. In the case of the Konishi operator, due to the high precision of the numerical data, we could give numerical predictions to the values of two further coefficients.
Beyond the numerical investigation of the anomalous dimensions, we investigated numerically the strong coupling behaviour of the coefficients c a,n in the power series representation of the P a -functions. Based on our high precision numerical data, in the regime, where the index of the coefficients is much smaller than λ 1/4 , we conjectured that the coefficients have polynomial index dependence at strong coupling. This allowed us to propose a strong coupling series representation for the P a -functions being valid far enough from the real short cut. To get some insight into the behaviour of P a close to the real branch cut, we also investigated the qualitative strong coupling behaviour of the P-functions at the branch points.
The paper is organized as follows: In sections 2. and 3. we recall the Pµand Qω-descriptions of the states under consideration and explain, how the free parameters coming from the symmetries of the QSC are fixed. The next section contains the detailed description of the numerical method together with all necessary technical subtleties which make it possible to implement the numerical code in C++ programming language. The analysis of the numerical data is presented in sections 5. and 6. The paper is closed by the summary of our results. Some technical details of the numerical method and some tables of numerical data are placed into the appendices of the paper.

Preliminaries
In this paper adapting the method of [36], we solve numerically the QSC equations for some twist-2 operators in the sl(2)-sector of the theory. The corresponding operators can be schematically represented as: O = Tr(D S Z L ) + . . . , (2.1) where Z is a complex scalar field of the theory, D denotes the light-cone covariant derivative, L is the twist, and S is the spin of the state. Here we investigate the case when L = 2 and S, the spin of the state, is even. The reason for this choice is to avoid treating null vectors in the internal linear problems of the numerical method (See remark at the end of subsection 4.2). So that we could use the high order perturbative results of [32] as initial values for the numerical iterative algorithm, we parametrized the P-functions and fixed the symmetries of the Pµ-system in the same way as it was done in [32]. Now, we recall the most necessary equations and relations of the QSC framework.
The QSC method [30,31] describes the full planar spectrum of AdS 5 /CF T 4 by the solutions of a set of nonlinear Riemann-Hilbert equations. The fundamental objects of QSC are the eight P-and Q-functions which separately form a basis on the 2 8 element of the Q-system of AdS 5 /CF T 4 . In the sl(2) sector, due to the left-right symmetry of the T-hook, one can describe the whole Q-system by only four P a , a = 1, .., 4 or four Q i , i = 1, ..4-functions, such that the other four (upper indexed) components are simple linear combinations of them: P a = χ ab P b , P a P a = 0, a = 1, ..., 4 (2.2) where χ is a constant matrix: The P a and Q i functions are analytic in the spectral parameter u with branch cuts.
The positions of the branch points depend on the 't Hooft coupling: λ and they may be located at u = ±2g + iZ, where g = √ λ 4π . All branch points are assumed to be of square root type. This means that, the result of two subsequent analytical continuations around a branch point is an identity transformation. The advantage of the choice of P a s or Q i s as basis is their very simple discontinuity structure. On the complex u-plane, P a has a single short cut, while Q i has only a single long cut, such that the discontinuities lie on the real axis.

The Pµ-system and the H-symmetry fixing
Since the states we study lie in the left-right symmetric sl(2) sector of the theory, we specify the presentation of the Riemann-Hilbert equations of the QSC for this sector. For any function f (u), denotef (u) the analytical continuation around the branch point ±2g and for short f [±n] (u) stands for f (u±i n 2 ). Then the Pµ-equations take the form [30]:

6)
µ ab = µ [2] ab , where µ ab = −µ ba and (µχ) b a = µ ac χ cb . The equations are valid in the strip 0 < Imu < 1, and elsewhere by their analytical continuations. In this representation µ ab has infinitely many short cuts and as a consequence of (2.5-2.7), it satisfies the Pfaffian-relation: Pf(µ) ≡ µ 12 µ 34 − µ 13 µ 24 + µ 14 µ 23 = 1. (2.8) In the sl(2) sector µ 14 = µ 23 . For twist-L states, the large u behaviour of P a and µ ab is fixed to [30]: where S is the spin of the state and ∆ is its conformal dimension. In addition the prefactors are constrained by the relations: Following the lines of [32] we also introduce the p a functions by a rescaling of the original P a s; is the short cut solution of the equation x + 1 x = u. By the introduction of p a , the sign ambiguity arising in the cases of odd L can be eliminated. In addition to the previously listed equations and properties, analyticity constraints are also imposed on the possible solutions of (2.5-2.7). Namely, in the QSC formulation of the spectral problem of AdS 4 /CF T 5 correspondence, it is postulated [30] that P a and µ ab have no poles on the first sheet and their absolute value is bounded at the branch points. The Pµ-system (2.5-2.7) is invariant under the linear redefinitions (H-symmetry [32]): where H is a constant matrix with detH = 1. In principle H might have 15 components, but if one would like to preserve the prescriptions (2.9) for the large u asymptotics, then only 6 non-zero elements remain to be fixed. These elements can be fixed by fixing the values of A 1 and A 2 and by imposing the value of 4 other coefficients in the large u expansion of p a . In our numerical framework, we used the H-symmetry fixing conditions of [32]. The requirements are as follows: • p 2 has no term proportional to u −1 in its large u expansion, • p 3 has no term proportional to u 0 in its large u expansion, • p 4 has no terms proportional to u 0 and u −1 in its large u expansion.
We used this H-symmetry fixing scheme, so that we could use the high order perturbative results of [32] as initial values for our numerical iterative algorithm. Nevertheless, since we study left-right symmetric states, also parity symmetries can be imposed on the first sheet. For the twist-2 case, we required that on the first sheet: • P 1 is even and real 3 function of u.
• P 2 is odd and real function of u.
• P 3 is even and imaginary 4 function of u.
• P 4 is odd and imaginary function of u.
These conditions allow us to use the following series representations for the p afunctions at L = 2: 14) The coefficients c a,n are functions of the coupling constant g. In our case c 1,n and c 2,n are real, while c 3,n and c 3,n are pure imaginary 5 . In (2.14) the leading terms of the 1/x expansion are fixed by the H-symmetry fixing conditions A 1 = g 2 and A 2 = 1.
In (2.15) A 3 and A 4 are considered as functions of ∆ and g, if we express them by the fixed A 1 = g 2 and A 2 = 1 coefficients through (2.10). These series representations automatically satisfy all the symmetry requirements discussed above and converge on the entire u plane [36]. The radius of convergence in 1/x is R = |x s (2 + i g )|. As a consequenceP a can also be represented by the analytical continuation (x → 1/x) of the series (2.14) and (2.15), but it is not convergent on the entire u plane. Its convergence is restricted to a oval domain lying around the real short cut of p a [36].
In the numerical solution, p a s are represented as truncated versions of (2.14) and (2.15), thus only a finite number of coefficients are to be determined. The concrete numerical solution of QSC [36] is implemented through the Pωsystem. This means that starting from the P a functions, one should determine 3 The Qω-system and its relation to the Pµ-system The nonlinear Riemann-Hilbert equations for the Qω-system are very similar to those of the Pµ-system [31]: where ω ij = −ω ji and (ωχ) j i = ω ik χ kj . The equations are valid in the strip 0 < Imu < 1, and elsewhere by their analytical continuations. In this representation ω ij has infinitely many short cuts and as a consequence of (3.1-3.3), it satisfies the Pfaffian-relation: In the sl (2) sector ω 14 = ω 23 . For large u, ω ij tends to a constant and the large u asymptotics of Q i is governed by the global charges of AdS 5 [31]: ). (3.5) In the sl(2)-sector, the prefactors B i satisfy an equation similar to (2.10): This means that fixing two of the coefficients B i is in our hand. For the sake of brevity, we introduce the vectors [36]: Then the large u asymptotics can be given by the short formulae: The Q-functions can be constructed from the P-functions in the following way. First, one should find 16 upper half plane analytic functions Q a|i as solutions of a set of homogeneous linear difference equations: The index i of Q a|i labels the 4 linearly independent solutions of (3.10). Then the Q-functions are defined by the formula: Since Q a|i is upper half plane analytic, the determination ofQ i is simple: As a consequence, (3.10) can be rephrased as follows: From this equation the leading order large u behaviour of Q a|i can be determined [31]: (3.14)

The brief description of the numerical method
The strategy of the numerical method is as follows [36]. One starts from the series representations (2.14,2.15) of P a and the goal is to compute numerically ∆ and those coefficients of the series, which are left undetermined after fixing the symmetries of QSC. Then from the representations (2.14,2.15),P a can be determined by an x → 1/x transformation. This representation ofP a is convergent in an oval shaped region containing entirely the branch cut on the real axis. The next step is to solve the recursion for Q a|i . This is done in two steps: first it is solved in the large u limit, and then the recurrence relations (3.10) are used to pull back the solution to the real axis. Then Q i andQ i are constructed from (3.11,3.12).
In order to exploit the Qω-equations, one has to determine ω ij , as well. It is computed from Q i andQ i by an integral expression derived from (3.1) and (3.3) (See (4.24) later).
All the quantities computed so far, are considered as functions of ∆ and the unknown coefficients of the series (2.14,2.15). This discrete set of variables is determined by imposing the equations (3.2).
In practice the whole process goes iteratively. One starts from a "good" approximation for the unknown coefficients and ∆, and goes through the steps discussed above. By the solution of (3.2), one gets the new initial values for the unknowns and the procedure is repeated until convergent result is obtained.
In the next section we describe the numerical method in detail, this is why the reader, who is interested in only the numerical results, might skip the next section.

The numerical method
In this section we describe our implementation of the numerical solution of QSC equations. We try to write down all important details and subtleties, in order to give help to those, who would like to solve numerically QSC equations in a fundamental programming language like C++ or Fortran. The technical details, we are going to write down, help to reduce each step of the numerical method to solving linear equations and to summations. The numerical implementation of these two simple mathematical problems is quite straightforward in any fundamental programming language.

Initial values and the discretization
In the previous section we described the set of unknown coefficients to be determined by the numerical method. The H-symmetry of the Pµ-system was partly fixed by fixing the values of A 1 = g 2 and A 2 = 1. Then A 3 and A 4 are given by (2.10) and they depend on A 1 , A 2 and ∆, provided L and S are fixed previously. As we mentioned, this choice of H-symmetry fixing was made to be able to use the perturbative results of [32] as initial values. Thus, for the twist-2 states with even S, in the weak coupling regime, where g 1 4 , we used the six-loop perturbative results of [32] for the unknowns as initial values for the iterations. According to our experience beyond the radius of convergence of the perturbative series (i.e. g = 1/4), the perturbative results were not good initial values for the iterations anymore. For 1 4 g, the numerical method failed to converge if we used the high loop perturbative results of [32] as initial values. For higher values of the coupling constant g, the initial values of the unknowns should be made out of the numerical data belonging to smaller values of g. This means that beyond g ≃ 1/4, one should increase g in small steps, and the initial values should be determined as appropriate compositions of the previously computed data. In our concrete numerical studies, we increased g with ∆g = 0.1, 0.05, 0.02, 0.01 and the initial values were given by a 4, 5, or 6 order Taylor-series composed of the previously computed numerical data. This construction of initial values is given in appendix A.
Since the numerical method uses also the Qω-system, we have further freedom to fix 2 of the coefficients B i . We fixed the values of B 1 and B 2 , then B 3 and B 4 are completely determined by (3.6). For the sake of simplicity, for small g we used the choice: For higher values of g, the choice of these coefficients play important role in the convergence of the numerical algorithm. Our experience suggests decreasing their values as g is increased. For example, in case of the Konishi operator (S = 2) the B 1 = B 2 = 1/g 2 choice was necessary 6 to reach satisfying convergence in the regime g > 2.
So far we explained, how to fix the "free" coefficients and how to construct good initial values for the iterative numerical algorithm. The next step is to choose the discretization points for our functions. The final equation (3.2) is imposed on the short cut of the real axis, this is why we need to give an appropriate discretization of the interval [−2g, 2g]. The discretization should be dense enough to be able to compute the integral expressions for ω ij with high enough numerical precision. Since all functions in the QSC framework have square root-type behaviour at the branch points, it is plausible to choose the discretization points as zeros of the Chebyshevpolynomials. The reason is that on the interval [−1, 1] the Chebyshev-polynomials of the second kind form an orthonormal basis with respect to the square-root type weight function √ 1 − u 2 . A summary on the necessary properties and identities of the Chebyshev-polynomials is given in appendix B.
In order to be able to use the advantages of formulae (B.9) and (B.10), the discretization points are chosen to be the zeros of the appropriately scaled 7 l c th Chebyshev-polynomial of the first kind (T lc ( u 2g )). The integer number l c measures, how dense the discretization is. Then the formula for our discretization points reads as 8 :

The determination of Q a|i
The necessary values: Q a|i (u A + i 2 ), A =, ..., l c are determined by (3.10) in two steps. In the first step, (3.10) is solved analytically for large u in the context of a 1/u expansion. One introduces an integer index cutoff N I , such that the first N I terms of the 1/u series are computed. Then another integer truncation index N u is introduced, such that at the points u ′ A = u A + i (N u + 1 2 ), the series representation of Q a|i truncated at N I , should approximate Q a|i (u ′ A ) within the required numerical accuracy. Then, in the second step, the desired discrete values Q a|i (u A + i 2 ), are computed from Q a|i (u ′ A ) by the successive application of the recurrence relation (3.10).
In the large u regime the following series representations are used: Scaling means only a u → u 2g scaling of the argument, such that the polynomial to be defined on [−2g, 2g] instead of the usual interval of definition [−1, 1]. 8 The same set of discretization points were chosen in [38].
As a consequence of the parity symmetries of P a , only even powers of u appear in the sums. From (2.2) it follows that: A a = χ ab A b and k a n = |χ ab | k b,n . The relation among the coefficients of the 1/u (4.4,4.5) and the 1/x (2.14,2.15) expansions can be computed by the x ↔ u relation:  Formulae (4.6) and (4.7) are valid for non-integer values of n, as well. In the twist-2 case the concrete forms of the k a,n ↔ c a,n relations read as follows:  where with α a|i = −M a +M i . The source term is the difference of two terms: we note that throughout the paper, in case the letter i stands for an index, than it denotes a positive integer number running from 1 to 4. In any other cases it denotes the imaginary unit i.e. i 2 = −1. The solution of (4.13) for m = 1, .., N I ,through (4.3), gives a numerically accurate approximation of Q a|i (u A + i(N u + 1 2 )). Then Q a|i (u A + i 2 ) is computed by the successive application of the recurrence relation (3.10): (4.19) where the 4 × 4 matrix U(u) is given by [36]: With the help of (3.11) and (3.12) it is easy to determine Q i andQ i at the discretization points: The +i 0 prescription is to avoid the evaluation of functions on their branch cuts. When one takes the series representations (2.14,2.15) at u A + i 0, it is better to use the mirror x, the long cut version of x, since it is regular in [−2g, 2g]: We close this subsection with a remark, which explains why we choose even integer values for S in the numerical studies. The reason is that in case of left-right symmetric states: det T a b|i m ∼ S ± 2 m − 1, which 9 means that for odd values of S, one should take care of the zero modes of T a b|i m . This problem is absent in the even S case.

The computation of ω ij
For the numerical algorithm we need to determine ω ij at the positions u A +i 0. From (3.1) and (3.3) the following integral representation can be derived [36]: where ω ij accounts for the discontinuity relations and periodicity, and ω c ij is a constant matrix to fulfill (3.2) close to infinity [36]: In the sl(2)-sector, the antisymmetry of ω c ij is ensured by I 12 = I 21 = I 14 = I 41 = I 23 = I 32 = I 24 = I 42 ≡ 0. In [36], it was explained that for numerical purposes, instead of using (3.2) as a final equation to fix the unknown coefficients, it is better to use a more regular version:Q where ω reg ij (u) = 1 2 (ω ij (u) +ω ij (u)) has no branch cut along the real axis. Our task is to compute ω reg ij (u A ), A = 1, ..., l c from the, so far computed, discrete set of The strategy goes as follows. Since Q i andQ i are bounded at the branch points ±2g, their antisymmetric combination can be represented as: where ρ ij (u) is a smooth bounded function on the real short cut. This allows one to represent ρ ij (u) as a convergent series with respect to some sequence of orthogonal polynomials.
For practical purposes explained in appendices B and C, we choose the Chebyshevpolynomials of the second kind U n ( u 2g ) as basis for this expansion: As a consequence of the convergence of this series, the coefficients quite fast tend to zero. Thus, ρ ij can be computed very accurately from the appropriately truncated version of (4.29). If the first l c terms are left from (4.29) after truncation, than the coefficients can be computed from well known formulae for the Chebyshevpolynomials. First, we introduce the matrix: Then we compute the expansion coefficients with respect to the Chebyshev-polynomials of the first kind: and finally using the identity (B.6), the coefficients of (4.29) are given by:  ij by the formulae: where T n denotes nth Chebyshev-polynomial of the first kind, and the expression of I ij entering ω c ij is also simple in terms of a (n) ij : One can recognize that in (4.34) the multiplier of a (n) ij depend on only g and the discretization points u A . This is why it is useful to compute it at the beginning of the numerical method. The computation of the quantity: involves an infinite sum. The numerical method for computing it within a given numerical accuracy, is described in appendix D.
The coefficients of (2.14,2.15) and ∆ are determined by imposing the equations: Instead of solving numerically (4.37) as an equation, [36] proposed to solve it as an optimization problem. This means that one tries to find the numerical solution of (4.37) by minimizing the quantity: This is performed by the Levenberg-Marquardt algorithm described in detail in the next subsection.

The Levenberg-Marquardt algorithm
The minimization of S is achieved via the Levenberg-Marquardt algorithm. To describe it, we put all unknowns into a single vector c. In our case certain unknowns are real 10 , while others are pure imaginary 11 . The real unknowns are put into the first Λ 1 components of c, while the other components are the imaginary ones: If we truncate the sums in (2.14) at N 0 th term, then the number of real unknowns is The reason is that the number of coefficients in the truncated versions of (2.14) is 2 N 0 , plus 1, because ∆ is also a real unknown. If the sums in (2.15) are also truncated at the N 0 th term, then the number of imaginary components is Λ − Λ 1 = 2 N 0 . Thus, if all infinite sums are truncated at the N 0 th term, then c is .., l c and denote F I = F i (u A ). In this notation (4.38) takes the form: and our task is to find the vectorc, which minimizes S(c). Assuming that c is close toc, S(c) can be linearized around the minimum and the minimization process consists of subsequent iterative minimizations of the linearized approximations of S(c).
To expand (4.39) around the minimum one needs to compute the derivative matrix: In practice it is done with the help of a second order formula for the first deriva- numerically approximated by the formula: It is worth to introduce its sign modified conjugate: If c is close to the minimumc of S(c), then using a linear approximation: and imposing the minimum condition S(c) ∂c k = 0, one gets a set of linear equations for the components of the minimum vector: In practice, during the iteration, equation (4.44) determines the new values of the unknowns from the old ones. Namely, if c (n) denotes the value of c after the nth iteration, then its value after the n + 1st iteration is given by: The iterational prescription (4.47) works very well, if the initial value of c is very close to the exact solution. Otherwise, it does not define a convergent iteration. In such cases the Levenberg-Marquardt (LM) modification of (4.47) is needed to decrease the difference |c (n+1) − c (n) | at each step of the iteration [36], and so to slow down and stabilize the iteration process. In the Levenberg-method, equation (4.47) is modified by adding a unit-matrix multiplied with an iteration number dependent number to M. In case of the Marquardt-method the unit-matrix is changed to the diagonal part of M: where λ (n) is an iteration number dependent number. The main drawback of the Levenberg-Marquardt modification is that, it defines a quite stable, but very slowly converging algorithm. To find the minimum of S(c) within practically acceptable amount of time, the term proportional to λ should be switched off after a few number of iterations. Here, we have to mention, another important property of the LM-algorithm, namely the larger the value of λ, the slower the convergence is. This is why, it is also desirable to decrease the value of λ at each step of the iteration. Taking into account the facts and experiences above, we used the LM-algorithm in the following way: First, we choose a not too large initial value for λ (0) and a divisor ν > 1. For the states under consideration we took λ (0) = 2.1 and ν = 2.0. At the nth step starting from c (n) , we go through the whole iteration process with λ (n) and get the new vector c (n+1) . If S(c (n+1) ) < S(c (n) ), then we decrease the value of λ by dividing it by ν, i.e. λ (n+1) = λ (n) ν . Otherwise we increase the value of λ by multiplying it by ν: λ (n+1) = λ (n) ν and the new iteration starts from the old initial values i.e. c (n+1) = c (n) . After a certain number of such iterations, when S(c (n) ) becomes small enough (∼ 1), the action of λ is switched off and the further iterations are done with the λ (n) ≡ 0 formula (4.47). We note that in our concrete numerical computations we used the Marquardt-type (4.48) modification of (4.47) and in practice we do not compute the inverse of M, but solve the following set of linear equations for c (n+1) :

The complete algorithm
In this subsection we write down the process of the numerical algorithm.
• Going through the process described in the previous subsections, we compute F I (c (0) ).
• To compute the derivative J jk , one does the same computation another 2 Λ times, but starting from the 1-component shifted initial value vectors: Λ }, where H = ±h or H = ± i h depending on the properties of c k under complex conjugation.
• Then the quantities J jk ,J jk , M jk , v j and S(c (0) ) are computed.
• The corrected values of the unknowns (i.e. c (1) ) are computed by the Marquardtversion of (4.47).
• The initial values of the next iteration are chosen by the rule: ν and the next iteration starts from c (1) . Otherwise λ (1) = λ (0) ν and the new iteration starts from the old initial values i.e. c (1) = c (0) .
• The whole process starts from the beginning...
• After several such iterations λ is set to be zero, and (4.47) determines the new approximations for the unknowns.

Numerical results for the Konishi operator
The Konishi operator is the most studied element of the set of single trace operators in the N = 4 super Yang-Mills (SYM) theory. The set of twist-2 operators with even spin also includes it as the L = S = 2 special case. In this section we summarize our numerical results obtained for the Konishi operator.
We solved the QSC equations in the range g ∈ [0.1, 7.0] and by fitting the numerical data, we determined numerically the first few coefficients of the large g expansion of some important quantities. Previous numerical investigations [23,24,25,26] could determine the first few coefficients of the large g series of the anomalous dimension ∆. Now, beyond the numerical determination of the coefficients of the strong coupling series of ∆, we also determine the large g behaviour of the coefficients of the 1/x series in (2.14,2.15). We also study the strong coupling behaviour of the p a functions around the branch points u = ±2g.
We note that the numerical data for ∆(g) and c a,n (g) are available in the corresponding text file 12 uploaded together with the paper. The pure numerical data can be read in a Mathematica notebook with the help of the DATAIN.nb notebook file 13 , where it is also explained, how to get a required quantity out of the huge array of numerical data.

Numerical results for ∆
We are interested in the coefficients of the strong coupling expansion of ∆: For the twist-L operators in the sl(2) sector, there are analytical predictions for the first four coefficients of (5.1). The coefficients depend on L and S and take the form [34]: The name of the corresponding text file is: L2S2data.txt. 13 It is also uploaded with this paper.
The first two coefficients in (5.2) can be determined either from Basso's slope function [43] or from semi-classical computations in string theory [40,41,42]. The next two coefficients were determined by matching the O(S 2 ) term of the small spin expansion with classical and semi-classical results [34].
To determine numerically the coefficients in (5.1), we computed ∆ numerically in the range g ∈ [0.1, 7.0] range with approximately 20 digits of accuracy and in the range g ∈ [4.6, 7] we fitted the numerical data with a power series of the form of (5.1).
The fitting method went as follows. We fitted a power series of type (5.1) to the numerical data. We increased the order of the truncation of the series until the numerical values of the coefficients stabilized. First, we concentrated on the first coefficient ∆ (0) . We experienced that it is very close to the exact value (5.2). This is why we assumed that its value is equal to the analytical prediction. Then we subtracted ∆ (0) λ 1 4 from the numerical data and fitted the new set of data with a truncated power series of type ∆ (1) ... Again, we increased the order of the truncation of the series until the numerical values of the coefficients stabilized. Then we concentrated on the coefficient ∆ (1) . We experienced that, the fitted value of the coefficient ∆ (1) is very close to the analytical prediction given by (5.2). Again, we assumed that the exact value of ∆ (1) is given by (5.2), and we subtracted also the second term of (5.1) from the numerical data. Then to get ∆ (2) , we fitted the new set of data with a series starting at of order λ − 3 4 etc. Our results for the fitted values of the coefficients of (5.1) are shown in table 1. The numerical data confirms with high precision the analytical predictions for the n = 0, 1, 2, 3 cases. Table 1 contains fitted values for the n = 4, 5 cases as well.
Since so far there are no available analytical predictions for these coefficients, we gave numerical estimations for further two previously unknown coefficients of the strong coupling expansion of the anomalous dimension for the Konishi state .
In table 1. δ rel ∆ (n) denotes the relative error defined by  f itted . Apart from fitting the coefficients of the strong coupling expansion of ∆, we also constructed a Pade-approximation like formula for ∆. According to our estimation, our approximation formula gives the values of ∆ with 14-digits of accuracy in the range of available numerical data i.e g ∈ [0.1, 7.0] and with at least 9-digits of accuracy for g > 7.0. The actual form of the Pade-approximation like formula for the anomalous dimension of the Konishi state can be found in appendix E.

The strong coupling behaviour of p a
In this subsection the strong coupling behaviour of the p a -functions is studied through the investigation of the strong coupling behaviour of the coefficients of the series (2.14) and (2.15). First, let us see, how the coefficients c a,n (g), look as functions of n at fixed g. Since the coefficients decay exponentially fast with a rate determined by the radius of convergence R(g) = |x s (2 + i g )| of the problem, for demonstrational purposes it is worth to introduceĉ a,n (g) by the definition: c a,n (g) = c a,n (g) R(g) 2 n+da , d a = δ a,1 + δ a,3 .

(5.5)
In order for the readers to get a taste about the n-dependence ofĉ a,n (g), we shoŵ c 1,n (g) at g = 4.4 in figure 1. In the other a = 2, 3, 4 cases, the picture is structurally very similar. The most important properties ofĉ a,n (g) at fixed g, can be summarized as follows: • The enveloping curve ofĉ a,n (g) has a power like decay with an exponent being close to 1.5. I.e.ĉ a,n (g) ∼ n −ǫa(g) , where ǫ a (g) ∼ 1.5 ± 0.2.
• Ifĉ a,n (g) is considered as a continuous function of n, then it has infinitely many zeros.
• In the large n regime the zeros are located periodically, such that the characteristic wavelength of this periodicity Λ a (g) ∼ a 0 √ g at strong coupling, with a 0 ∼ 4.4.
One can recognize another interesting property of the coefficients, if one plotŝ c a,n (g) at all available values of g on the same plot. They all have very similar shape, which suggests that in the strong coupling limit they can be transformed into a universal g-independent function with some scale transformation. Indeed, figures 2 and 3 show that the transformed coefficients g −naĉ a, √ gν with (n 1 ,n 2 ,n 3 ,n 4 ) = (1, 0, 3, 2) tend to universal g-independent functions K a (ν) at strong coupling. For later purposes, we write it down in a formula as well: where the dots stand for negligible terms for g → ∞.
This fact shows that the in the strong coupling limit the relevant scale of the problem is given by √ g or equivalently λ 1 4 as it is expected from the strong coupling behaviour of the anomalous dimension.

Strong coupling behaviour of c a,n for fixed n
In this subsection we investigate, how the coefficients of the series (2.14) and (2.15) behave at strong coupling, if we fix the value of the index n. We considered the first 12 or 14 coefficients of the series (2.14) and (2.15). I.e. c a,n with a = 1, .., 4 and n = 0, ..., 14. Then in the range g ∈ [4.6, 7.0] we fitted the numerical data with a series 14 in 1/g. Our numerical data was consistent with the series expansions as follows: where the integer leading power n a and the numerical values of c all coefficients behave in the same way, and this leading order power behaviour is determined by the H-symmetry fixing condition. The situation is very similar in the a = 3, 4 cases. There the leading powers are the same as those of A 3 u = A 3 g (x+ 1 x ) and A 4 u 2 = A 3 g 2 (x + 1 x ) 2 with x being fixed. From (2.10) and (5.1) it follows that, at large g: Next, we can concentrate on the first, leading order coefficients 16 c (0) a,n in (5.7). Table 2. shows their fitted values. Looking at the data, one can recognize the remarkable fact that for fixed values of the index a, and for n ≥ 1+δ a,4 the coefficients c (0) a,n seem to be n-independent. The difference between the numerical values of the columns are supposed to be the consequence of numerical errors. Then, it is tempting to guess the exact values of c (0) a,n from the available numerical data of table 2. It is not hard to make good proposals for the cases a = 1, 2: a,n for a = 3, 4 seem to be more difficult, but the following train of thoughts leads to reasonable proposals. One can recognize that based on (5.9), in the case of a = 1, 2, in (2.14) all 1/x powers has the same 14 We tried to fit other types of series in g, like series in 1/ √ g etc., but only the 1/g case gave numerically stable coefficients. 15 We note that n a =n a of (5.6) for a = 1, 2, 3, 4. 16 We just recall that c Then substituting u → g(x + 1 x ) into (2.15) and imposing that the coefficients of each 1/x power are equal, one gets the analytical predictions:  a,n , we can sum up the emerging geometrical series and give analytical formulae for the leading order large g behaviour of the functions p a . The results of the summations take the forms: 14) 17 In leading order for large g.
The above formulae has the common property that they have poles at x = ±1. The positions of these poles are in accordance with the g → ∞ limit of the radius of convergence R. Nevertheless, there are two facts, which indicate that (5.14,5.15,5.16) cannot be good approximations of the functions p a on the entire u-plane at strong coupling. First, in (5.14,5.15,5.16) the neglected terms are O(1/g) with respect to the leading ones, in case the multipliers of 1/g in the correction terms are bounded functions of u with g independent upper and lower bounds. We will see in the next subsection that this is not the case.
Another problem, which indicates the restricted validity of (5.14,5.15,5.16), emerges when one would like to computep a at strong coupling. Naively, it can be done by a simple x → 1/x transformation in (5.14,5.15,5.16). But the result does not account for the thep a (u) ∼ u 4 √ π g+... large u asymptotics expected from (2.9) and (5.1,5.2).
The main reason for these discrepancies is that the coefficients c a,n (g) depend on n and g. This is why the result of the g → ∞ limit depends on the relative magnitude of these two variables. In the expansion (5.7) we considered the limit, when n ∼ 1 and g → ∞.
To be more precise, we will see later that, the n ≪ √ g limit is the one, which corresponds to the expansion (5.7).

Terms beyond the leading order
From the available numerical data, one can fit further coefficients in (5.7), as well.
We determined numerically the coefficients c implies that c (k) a,n ∼ n 2k at large n. The simplest function, which accounts for this behaviour is a polynomial of order 2k. Indeed, table 3. and the tables of appendix F. show that the numerical values of c (k) a,n can be perfectly described by polynomials of order 2k. This is why, we make the following conjecture: • The coefficients c   of the polynomial Ansatz (5.17). ∆P rel is the relative error measuring, how precise the polynomial description of the various coefficients.
As a consequence, the polynomials can be given by 2k+1 n-independent parameters, which, for practical purposes, we parametrized as follows: We note that in the k = 0 special case, by definition α (m) a,a δ m,1 and that (5.17) can be used only when n ≥ 1 + δ a,4 .
The conjectured (5.17) representation of c (k) a,n implies the following series representation for p a (x) at strong coupling: where A 3 (g) and A 4 (g) admit the strong coupling series representations: The first few values of A 3 (g) and A 4 (g) are given in the table 4. All elements of table 4 are small numbers, lying in the range of numerical errors. This fact suggests us to make the following conjecture: • A   and A (k) 4 . All values are in the magnitude of the numerical errors.
As a consequence A 3 (g) = A 4 (g) ≡ 0, which implies that besides of the 1 (x 2 −1) m type of terms, there are no 1 x or 1 x 2 terms present in the strong coupling series (5.20). The formula (5.20) indicates that in the a = 1 case there is some simplification due to the H-symmetry fixing condition c 1,0 ≡ g. This implies that in the large x expansion of (5.20) the coefficient of 1 x does not get 1 g corrections. As a consequence: α (1) 1,k ≡ 0 for k ≥ 1. This means that in the a = 1 case only 2k parameters describe the conjectured polynomials of order 2k. This fact was built in the polynomial fits as it is demonstrated by table 3.
Reshuffling the series part of (5.20), it can be written as a series in 1 g(x 2 −1) 2 : p a (x) = δ a,2 +δ a,3 gA 3 (g)x + A 3 (g) x +δ a,4 g 2 A 4 (g)(x 2 + 2)+ A 4 (g) Now, we are in the position to discuss the regime of validity of (5.22) in the rapidity plane. Formula (5.22) implies that at strong coupling the variable z = 1 g(x 2 −1) 2 becomes relevant and within the range of convergence, apart from sum trivial factors, p series a (x) can be represented as a sum of functions of z, such that each function is suppressed with an inverse power of g: p series a (x) = g na x δ a,1 +δ a, 3 1 To study the range of validity of (5.22), one has to determine the radius of convergence of the series representations of f odd a,0 (z) and f even a,0 (z). We just recall: The radius of convergence of these series is determined by the large n behaviour of the coefficients. Our numerical data suggests that: for large n. This implies that the radius of convergence of f odd/even a,0 (z) is 1 4 . Thus one can conclude that the validity of the series representation (5.22) is restricted by the inequality: In the strong coupling limit, (5.25) may fail, if x is close to ±1. In the language of the rapidity 18 u, this means that u is close to the branch points ±2. Using the series representation: Thus, naively one might conclude that the series representation (5.22) gives the correct strong coupling approximation of p a in the domain where, the distance of the rapidity u from the branch points is larger than 4 g . Unfortunately the situation is a bit worse. The series (5.22) will be an appropriate strong coupling approximation for p a (u) only outside of an oval region containing the real short cut [−2, 2], such that the horizontal dimension of the oval region is 4 plus a number of order 1 g , and its vertical dimension is of order 1 √ g . See figure 4. The reason is as follows. Rephrasing (5.7) one obtains that: c a,n (g) = g na K a n √ g · (1 + O( 1 √ g )).
(5.28)  The O( 1 √ g ) magnitude of the corrections is a consequence of (5.17). From (5.28) it follows that the n = fixed, g → ∞ limit corresponds to the n √ g → 0 limit. This implies that the strong coupling series representation (5.7) of the coefficients is a good approximation until n ≪ √ g. (5.28) also implies that, at strong coupling a typical sum appearing in p a can be roughly estimated by an integral: (5.29) The strong coupling series (5.22) was obtained by inserting the series (5.7) into (2.14) and (2.15) and evaluating the sums from 1 to infinity. In this representation the strong coupling corrections go as inverse powers of g. Since the validity of (5.7) is restricted to n ≪ √ g, (5.22) can be appropriate representation of p a , if the neglected contributions coming from the √ g n region are exponentially small in g. As (5.29) shows, the exponentially small corrections grow up to power like in the regime, where √ g ln x or equivalently |x| − √ g becomes of order 1. Now we will show that this can happen in an appropriate neighborhood of the real short cut of the u-plane.
At the branch points, x is given by (5.26), therefore √ g ln x ∼ 1, when u lies within a circle of radius ∼ 1 g , whose center is located at the branch points ±2. On the other hand x is a pure phase on the real cut, i.e. |x| = 1. If u 0 ∈ [−2, 2], then ln x(u) can be expanded in a regular Taylor-series around u 0 . This yields that To summarize, the contributions of the √ g n terms are not negligible in (2.14) and (2.15) if u lies in an oval domain containing the real short cut [−2, 2], such that the horizontal dimension of the oval region is 4 plus a number of order 1 g , and its vertical dimension is of order 1 √ g . (See figure 4.) This is the region, where the strong coupling formula (5.22) becomes invalid. To be more precise, the neglected contributions of the √ g n terms are exponentially small outside of this oval domain, and become power-like inside the domain. Now, we have shown that conjecture (5.20) cannot be an appropriate approximation of p a close to the real short cut, this is why we also studied the behaviour of p a close to the branch points in the context of a series expansion in the deviation from the branch points.

Series expansion around the branch points
Now, we study the behaviour of p a at the branch points. Inserting the power series 19 (5.26) into the series representations (2.14) and (2.15), one ends up with the expansions: where we use the convention, when the rapidity is scaled, such that the branch points are at ±2 and v denotes the deviation from them. The coefficients β a,k (g) are certain linear combinations of the momenta 20 of the coefficients c a,n (g). For example the first coefficient is just the sum of the coefficients c a,n (g), i.e. β a,0 (g) = ∞ n=0 c a,n (g).
We fitted the coefficients β a,k (g) by a power series in √ g. The coefficients of the numerical fits proved to be stable with respect to increasing the truncation index of the series, in case the following g dependence was assumed: Its infinite order version. 20 Here, by momentum we mean sums like:    The numerical values of the first few coefficients γ (n) a,k can be found in tables 5. and 6. Concentrating on only the leading order behaviour of (5.30), the following pattern arises: a,k (g v) k/2 + ..., (5.32) where dots mean terms negligible for large g.
As a consequence we can conclude that for large g, close to the branch points p a behaves like a function of gv and the sub-leading corrections are suppressed by positive integer powers of 1 g :

Higher spin results
In this section we publish the numerical results obtained in the S = 4, 6, 8 cases. For these higher spin values, we could not reach as large values of the coupling constant g as it was done in the case of the Konishi operator. The reason for this, is that increasing the spin, the numerical algorithm becomes more and more sensible to the choice of initial values. This fact forced us to increase g in very small ∆g ∼ 0.02 steps. As a consequence, we needed to run 50 jobs subsequently in order to increase g with one single unit. Unfortunately, this process proved to be very time consuming.
By increasing S, also the internal precision of the computations must have been increased, in order to get convergence and reach the required precision for ∆ and c a,n . For example at strong coupling g 2.7, the S = 4, 6, 8 cases required 60-, 80-and 100-digits of precision respectively. The necessity of the application of such high precisions made also the runtime of the jobs very long. Because of these difficulties, in the S = 4, 6, 8 cases, the numerical results we obtained were less accurate than those of the Konishi state. This is why, in the higher spin cases, we restricted our numerical work to 3 types of investigations. Namely, • Numerical determination of the first 4 coefficients ∆ (n) of the strong coupling series of ∆.
• Numerical determination of the coefficients c (0) a,n of (5.7).
• Investigation of the qualitative strong coupling behaviour of the p a -functions at the branch points.
The fitted values of the coefficients in (5.1) at different values of S can be found in tables 7., 8.,and 9. The numerical estimations of the first coefficients beyond the analytical prediction (i.e. ∆ (4) ) are also presented, but only to "give a taste" about their magnitude. Though the precision of the coefficients is not so high as it was in the Konishi case, the first four coefficients can be compared to the analytical predictions (5.2), (5.3) and (5.4). Our numerical data confirms the analytical predictions within the range of numerical errors.
In the higher spin cases, we also computed numerically the first few coefficients from the set of c          Though the numerical values of the coefficients are not as accurate as they were in the case of the Konishi operator, one can see that the same structure shows up.
Namely, for n ≥ 1 + δ a,4 the coefficients seem to be n-independent. Using the same train of thoughts, as it was done in the Konishi case, based on the numerical data of tables 10., 11., and 12., we made the following proposals for the exact values of the coefficients: 2,n = 4 3 , n = 1, 2, ...

(6.2)
We also constructed Pade-approximation like formulae to determine numerically  a,n at S = 8.
∆ in the whole range the coupling constant. The Pade-approximation like formulae for the cases S = 4, 6, 8 can be found in appendix E. Unfortunately, these approximations are not so accurate as that of the Konishi operator. The reason for that is two-fold. First, because we did not reach too large values of g during our numerical work 21 . The second reason is the lower precision of the available numerical data. Nevertheless, according to our estimations, our Pade-approximation like formulae give the numerical values of ∆ with 8-digits of accuracy in the range, where numerical data are available, and with 4-5 digits of accuracy for higher values of g.
The last problem, we studied in the higher spin cases, is the strong coupling behaviour of p a functions at the branch points. Without listing any fitted numerical data, we just note that our numerical results suggest that close to the branch points the qualitative strong coupling behaviour of p a functions is given by (5.33). Thus, it is independent of the concrete value of the spin 22 .

Summary
In this paper, we solved numerically the QSC equations corresponding to some twist-2 single trace operators from the sl(2) sector of AdS 5 /CF T 4 correspondence. Namely, we considered the twist-2 operators with spins S = 2, 4, 6, 8. The primary purpose of the numerical study was to gain some information about the strong coupling behaviour of the solutions of the Pµ-system. 21 The largest values of g reached during the numerical work were 4.1, 3.5 and 2.74 in the cases S = 4, 6, 8 respectively. 22 At least in case the spin is an even and positive integer number.
We applied the numerical method of [36] to solve the QSC equations and we wrote down all technical details, which were necessary to implement the numerical code in C++ language. Roughly speaking, the whole numerical algorithm consist of summations and of numerical solutions of linear sets of equations. Both mathematical problems can be easily programmed in any fundamental programming languages.
The most accurate numerical results were obtained in the case of the Konishioperator. There, λ ∼ 7737 was the highest value of the 't Hooft coupling, which was reached by the numerical computations. From our high precision numerical data, we could numerically confirm the analytical predictions of [34] for the first 4 coefficients of the strong coupling series of ∆. Moreover, due to the high precision of the numerical data, we could give numerical predictions for 2 further coefficients in the strong coupling expansion of ∆. In the cases of S = 4, 6, 8 the numerical data were less precise, nevertheless they proved to be precise enough to confirm the analytical predictions of [34], though with much less precision. We also constructed Pade-approximation like formulas which allow one to compute the anomalous dimensions of the states under consideration within short time and with satisfactory high precision. (See appendix E.) Beyond the numerical determination of ∆, we also focused our attention to determine the strong coupling limit of the p a functions. Since, in the numerical method the coefficients of their series representations (2.14), (2.15) were the basic objects, we tried to determine the strong coupling behaviour of these coefficients. From the numerical data, we found that, at strong coupling, when n ≪ √ g , the coefficients admit the series representations (5.7) with n a given by (5.8). The accurate numerical values obtained for the coefficients of (5.7), inspired us to make analytical proposals for the values of the leading order coefficients (6.1,6.2).
For the Konishi operator, based on the high precision numerical data, we conjectured that the coefficients c (k) a,n in (5.7) are polynomials of order 2k+1 in n. This recognition led us to propose a strong coupling series representation (5.20) for the p a -functions 23 . We argued that (5.20) is an appropriate strong coupling represen- 23 The fundamental functions P a of the QSC method are connected to p a by the simple formula (2.11), this why the results given for p a in the previous sections, can be translated to the language of P a in a straightforward manner. tation of p a (u), if the rapidity u lies outside of an oval domain 24 containing the short real cut, such that its horizontal dimension is equal to 4 plus a number of order 1 g and its vertical dimension is ∼ 1 √ g . (See figure 4.) Furthermore, outside of this domain (5.20) accounts for all power like contributions in g, but neglects the exponentially small ones, which come from the index range √ g n.
Because of this restricted validity of (5.20), we also studied the behaviour of the solutions close to the branch points. The result of this investigation can be summarized by the scaling behaviour given by (5.33).
The strong coupling investigation of the numerical data suggested the strong coupling scaling behaviour (5.28) for the coefficients. This indicates that √ g is the relevant scale of the problem at strong coupling and it tells us that there are 3 important regimes of n in the strong coupling limit. These are the n ≪ √ g, n ∼ √ g and n ≫ √ g regimes. In the 3 different regimes the coefficients have different strong coupling behaviours.
We also discussed some general properties of the coefficients at fixed values of the coupling constant. If c a,n is considered as a continuous function of n, the numerical data implied that • that c a,n has infinitely many zeros located periodically at large n, and • that c a,n decays as ∼ n −ǫa(g) R −2 n at large n, where R = |x s (2+ i g )| is the radius of convergence of the series (2.14), (2.15) and ǫ a (g) is a numerical constant with approximate value ǫ a (g) ∼ 1.5 ± 0.2.
Our numerical work contributes to the deeper understanding of the strong coupling behaviour of the solutions of the QSC-equations and hopefully it will help in finding the a method for the systematic analytical solution of the Pµ-system in the strong coupling limit.
of Sciences and OTKA K109312. The authors also would like to thank the support of an MTA-Lendület Grant.

A Construction of initial values at strong coupling
For small values of the coupling constant g, the numerical iterations can start from the perturbative solution of the problem [32]. This strategy works for g In this appendix we describe, how to construct good initial values for the numerical iterations, provided we have the numerical solution of the problem for several smaller values of g. To construct good initial values, one should increase the value of g in small steps. We increased the value of g uniformly at each step by ∆g = 0.1, 0.05, or 0.02. If we assume that every unknown coefficient is a smooth function of g, then a good initial value of the numerical problem can be given by a numerical Taylor-series, constructed from the numerical data belonging to previous values of g. Here, let f a function of g. f should be considered here as the analog of any unknown coefficient of the numerical problem. E.g. ∆(g) is one of them.
In case ∆g is small enough, a good initial value can be constructed as a second order Taylor-series: For the numerical implementation of (A.1), one needs to compute the appropriately accurate numerical formulae for the derivatives: Inserting (A.2) and (A.3) into (A.1), and making the g → g − ∆g substitution, one gets the formula: 25 Good initial value means that the numerical algorithm converges if the process starts from it.
where for later convenience we introduced the notation: g n = g−n ∆g. By increasing the order of the Taylor-series method and using the same procedure, higher order formulae can be derived. Here, we list them up to the sixth order. The forms of the 4-, 5-and 6-order formulae take the form: Finally, we mention that, in case we had numerical data at least for six consecutive values of g, then we used the 6-point rule (A.7) to construct the initial values of the numerical algorithm for the next value of g.

B Chebyshev-polynomials
In this appendix we summarize some useful properties and integral formulae of the Chebyshev-polynomials. The Chebyshev-polynomials of the first kind T n (u) form a sequence of orthogonal polynomials on [−1, 1] with respect to the weight function: The orthogonality relation is given by the integral formula: (1 + δ n,0 ), n, m = 0, 1, 2, ...

(B.2)
For practical purposes, we define their slightly modified version: In the QSC method, close to the branch points, the relevant functions behave like 4 g 2 − u 2 . This is why, in our numerical studies the Chebyshev-polynomials of the second kind U n (u) become important, since they form an orthonormal basis on [−1, 1] with respect to the weight function √ 1 − u 2 . They can be given by the explicit formula: U n (u) = sin((n + 1) arccos u) sin(arccos u) , n = 0, 1, 2, ... (B.4) and the orthogonality relations they satisfy, read as: The two kinds of Chebyshev-polynomials are related by a simple recurrence relation: where U n (u) for n < 0 is zero by definition. According to the theory of orthogonal polynomials, on [−1, 1] any smooth function f can be represented as a convergent series in either T n or U n : As a consequence of (B.6), the coefficients are related by: where the discretization points are chosen to be zeros of T lc : Here, it is assumed that l c is so large that the coefficients with higher index are so small that they are irrelevant up to the numerical precision required. Thus the series is truncated at the index l c . In our actual numerical computations, the following integral formulae for U n are important: where x s is given in (2.12) and (B.12) contains a principal value integration.
C The derivation of formulae (4.33) and (4.34) In this appendix we show, how to use the Chebyshev-expansions to the derivation of the formulae (4.33) and (4.34) for ω ij and ω reg ij . First, we start with some remarks concerning the coefficients of the series representations (2.14,2.15).
Let f (u) be a function on C with the properties as follows: • It has no poles, • It has a single branch cut at [−2g, 2g] with square root type discontinuity.
• The discontinuity on the branch cut is given by i ρ(u).
• The discontinuity becomes zero at the branch points, which means that it behaves like ∼ 4 g 2 − u 2 at ±2g.
• f decays at least as fast as 1 u at infinity.
Then f (u) can be expressed by its discontinuity by the formula: Moreover, since ρ(±2g) = 0, it can be represented as: where ρ 0 (u) is a smooth regular function on [−2g, 2g]. This is why it can be expanded in a convergent series with respect to U n s: As a consequence of (C.1), (C.2), (C.3) and (B.11) f (u) admits the convergent series representation as follows: Consequently, we can conclude that the coefficients in the expansions (2.14) and (2.15) are nothing else, but the coefficients of the Chebyshev-series of the discontinuity functions 26 of p a s. In this sense the formulae (4.33) and (4.34) are the periodic analogs of (C.4).
To derive (4.33), first one has to rephrase the kernel as an infinite sum: Then inserting (C.5), (4.28) and (4.29) into (4.25) and evaluating the integrals with the help of (B.11) one ends up with (4.33).
To derive (4.34), one should represent ω reg ij by the formula: 0)). The derivation of (4.34) is very similar to that of (4.33). The only difference comes from the ∼ 1 u−v term of (C.5). Now, the ±i 0 prescriptions become important. If they are treated by the Sokhotski-Plemelj formula, only the principal value part remains. This principal value integral can be evaluated with the help of (B.12), which gives the term T n+1 (u) in (4.34).

D A method to compute (4.36) numerically
In the implementation of the numerical method for solving QSC equations, only such simple mathematical operations appear, like summations and finding the solutions of some linear equations. Both methods can be easily implemented in C++ language. There is only one subtle quantity Ω A,n (g) defined in (4.36), which requires the accurate computation of an infinite sum. In this appendix, we describe, how to reduce the computation of this quantity to finite summations, provided one needs the result with a given numerical accuracy. Here, we recall the definition of Ω A,n (g), where u A ∈ [−2 g, 2 g] are the discretization points. For the sake of simplicity, in the sequel we will omit the index A from u A . First, we sketch the idea of the numerical computation and the deeper technical details will be given in the subsequent paragraphs. For practical purposes, we introduce a short notation for the summand: (D.2) We introduce also an integer cutoff Λ X to write the infinite sum as a sum of two terms: 3) The first term in the rhs. of (D.3) is a finite sum, so it can be evaluated numerically by a computer. Since Λ X is chosen to be large, in the second term on the rhs. we can use the large k expansion of the summand. It defines a series in 1/k, and the explicit sums of the 1/k powers can be expressed by the Riemann-zeta function. To reach a given accuracy, only a finite number of terms of the 1/k series needed to be taken into account. If 1 k Nx is the last term, which is summed in the large k series, then the magnitude of the numerical error is ∼ 1 Unfortunately, this naive estimation needs to be corrected, when one takes a deeper look at the structure of the summand (D.2). This is why, in the next paragraphs, we write down in more detail the numerical computation of (D.1).
The first ingredient is the large k expansion of the summand I (n) X (k, u). It can be obtained by inserting the following two series expansions into (D.2): where κ (α) s is given by (4.7). The final form of the expansion takes the form: (D.6) allows us to make the appropriate choice for the cutoff parameters Λ X and N x . For the sake of simplicity concentrate on the power like terms in (D.6). A typical such term looks like ∼ g n−q u q k n . In the numerical algorithm, we need to compute (D.3) at the discretization points, which lie in the interval [−2g, 2g]. This is why we can give an upper estimation for this typical power-like term: This inequality tells us that, not the powers of 1/k determine the magnitudes of the terms in the 1/k series, but the powers of 2 g k . This means that, if 1 k Nx is the last term, we sum from Λ X to infinity in (D.3), than the numerical error can be estimated by 2g Λ X Nx instead of the naively expected value 1 Λ X Nx . Now, we are in the position to make a choice for the values of Λ X and N x . We require N c digits of accuracy for (D.3). This means that the estimated error term should be ∼ 10 −Nc . In accordance with the content of the previous paragraph, this requirement imposes an inequality among the parameters Λ X , N x and N c . 2g Λ X Nx 10 −Nc . (D.8) The value of Λ X is chosen to "maximize" the inequality: Certainly, (D.9) does not allow to determine both Λ X and N x . One of them is free to choose and the other one is given by (D.9). In our actual numerical computations, we made the choices:   Apart from the numerical values we listed in the tables, the interested readers can find all the numerical data we obtained, in the Mathematica notebook and text files attached to the text file of the paper.
Apart from fitting the strong coupling series coefficients of the the anomalous dimensions, we also used the numerical data to construct Pade-approximation like formulas in order to describe the anomalous dimensions of the operators under consideration at all values of the coupling constant with satisfying numerical precision.
Instead of the computation of an interpolating function composed of rational polynomials, we performed a nonlinear model fit to the data points. This approach gave smooth approximants for real values of the coupling constant, and could inform us about the validity of the approximation as well.
We found that fitting a naive rational polynomial approximation for ∆(g) does not give stable 29 values for the coefficients of the rational polynomial. This is not surprising, if one observes that in the perturbative expansion around g = 0 only even powers are present, while in the strong coupling regime the leading term is ∼ √ g and the corrections go as inverse powers of g.
To have an optimal form for the approximation, we basically followed the Ansatz used in [23]: 2 and g b is a suitable constant, whose value is chosen to be 2 in the case of the S = 2, 4, 6, 8 twist-2 operators.
In principle some analytical information can be built into the Ansatz from the perturbative results [32], by fixing some relations between coefficients. For practical calculations however, we exploited only the known value 30 of ∆(0) and the leading order strong coupling asymptotics of ∆(g) given in (5.2). These data fixed a 0 and the ratio of a n and b n .
Because of the high precision of the numerical data, an unusually high number of coefficients could be fitted. For the Konishi operator, we stopped at n = 15, where the coefficients seem to be still stable with respect to changing the value of n. We performed the fits by Mathematica's build in NonlinearModelFit function, which provides "prediction bands" 31 allowing one to infer to the accuracy of the Pade-approximation like formula, as well.
The measured points and the fitted curve are shown in figure 5.
Because of the small magnitude of the deviations, we show separately the residual plot of the data in figure 6. Figure 6 shows that the data points are so close to the fitted curve that the data points are approximated with the Pade-approximation like formula with 14 digits of accuracy.
To predict the accuracy of the fitted curve beyond the measured interval, we used Mathematica's build in "MeanPredictionBands" function and we set the confidence level to 99%. Figure 7 shows that even outside of the range of available numerical data, the fitted Pade-approximation like formula can be taken seriously up to 9 digits of accuracy.
Analogously to (E.1), Pade-approximation like formulas were constructed for the S = 4, 6, 8 cases, as well. The structure of the approximation formulae are the same as that of the Konishi operator, the only difference is the actual form of the rational           of the polynomial Ansatz (5.17). ∆P rel is the relative error measuring, how precise the polynomial description of the various coefficients.