Strong coupling results in the AdS5/CF T4 correspondence 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 AdS5/CFT4 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 ca,n in the power series representation of the Pa-functions is also investigated. Based on our numerical data, in the regime, where the index of the coefficients is much smaller than λ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 P-functions being valid far enough from the real short cut. In the paper the qualitative strong coupling behaviour of the P-functions at the branch points is also discussed.


JHEP08(2016)061
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: 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).

JHEP08(2016)061
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 Pand 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 Thook, 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]:

JHEP08(2016)061
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: 16iL(L + 1) , . (2.10) Following the lines of [32] we also introduce the p a functions by a rescaling of the original P a s; p a ≡ (g x) L 2 P a . (2.11) Here x ≡ x s (u/g), where x s (u) = u 2 1 + 1 − 4 u 2 , |x s (u)| > 1, (2.12) 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: • A 1 ≡ g 2 and A 2 ≡ 1, • 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.

JHEP08(2016)061
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 a -functions 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 4,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].

JHEP08(2016)061
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 the Q i functions of the Qω-system and the coefficients are determined from the discontinuity equations of the Qω-system. To do so, we have to recall the Qω-system and its relation to the Pµ-system.
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]: 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:

JHEP08(2016)061
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 Qfunctions 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.

JHEP08(2016)061 4 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 JHEP08(2016)061 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 Chebyshev-polynomials. 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 Chebyshevpolynomial 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: 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: 7 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]. 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 Θ is the unit-step function and where with α a|i = −M a +M i . The source term is the difference of two terms: where q ab n = n l=0 k a,n−l k b l q ab 0 ≡ 1, (4.18) and in the summation limits [. . .] stands for integer part. To avoid any confusion, we note that throughout the paper, in case the letter i stands for an index, then 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 10 ω ij at the positions u A +i 0. From (3.1) and (3.3) the following integral representation can be derived [36]:

JHEP08(2016)061
where ω (0) ij accounts for the discontinuity relations and periodicity, 25) 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, then the coefficients can be computed from well known formulae for the Chebyshev-polynomials. First, we introduce the matrix: Then we compute the expansion coefficients with respect to the Chebyshev-polynomials of the first kind:

JHEP08(2016)061
and finally using the identity (B.6), the coefficients of (4.29) are given by: (4.32) Using the results of appendix C, ω ij and ω reg ij can be expressed in terms of the coefficients a (n) 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, 11 while others are pure imaginary. 12 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 Λ 1 = 2 N 0 + 1. 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 Thus, if all infinite sums are truncated at the N 0 th term, then c is a Λ = 4 N 0 + 1 component vector. For short, we introduce the multi-index I = (i, A), i = 1, . . . , 4, A = 1, . . . , 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 tõ c, 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 derivative: with h being a small number. Thus J Ik (c) is 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:

JHEP08(2016)061
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 unitmatrix 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) ν .

JHEP08(2016)061
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: c • 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 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 JHEP08(2016)061 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) can be found in the DELTAdata.nb ancillary notebook file, while the numerical values of c a,n (g) can be downloaded from the arxiv site of the paper 13 [arXiv:1604.02346].

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 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 . .. 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 13 We did not attach the text files containing all numerical data, because of their very large size. 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 . For n = 4, 5 in the lack of analytical results, δ rel ∆ (n) was computed as the ratio of the estimated error for ∆ (n) fitted and ∆ (n) fitted . 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 description of our Pade-approximation like formula for the anomalous dimension of the Konishi state can be found in appendix G, and its actual form can be found in the ancillary approx.nb notebook file.

The strong coupling behaviour of p a
In this subsection the strong coupling behavior of the p a -functions is studied through the investigation of the strong coupling behavior 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: In order for the readers to get a taste about the n-dependence ofĉ a,n (g), as an example we showĉ 1,n (g) at g g, the picture is structurally very similar. Namely, at fixed g, the n-dependence is given by a "wave-like" function with an enveloping curve which decays at large n. The most important properties ofĉ a,n (g) at fixed g, can be summarized as follows (For more details see appendix F): • The enveloping curve ofĉ a,n (g) has a power like decay at large n:ĉ a,n (g) ∼ n − a(g) .
• 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.
• At fixed g the characteristic wavelength Λ a (g) of the large n periodicity is independent of the value of the index a.
Moreover our numerical data suggests the following strong coupling behavior for the characteristic wavelength Λ a (g) and the exponent a (g): (5) and dots stand for terms negligible at g → ∞.
• At strong coupling the powers a (g) tend to constant values, which lie in the interval 14 One can recognize another interesting property of the coefficients, if one plotsĉ a,n (g) at all available values of g on the same plot. They all have very similar shape, which JHEP08(2016)061  suggests that in the strong coupling limit they can be transformed into a universal gindependent 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 15 in 1/g. Our numerical data was consistent with the series expansions as follows: 15 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.

JHEP08(2016)061
n c where the integer leading power n a and the numerical values of c (k) a,n were determined from the fitting process. The best fits yield the following values for the leading powers: 16 (n 1 , n 2 , n 3 , n 4 ) = (1, 0, 3, 2). (5.8) For a = 1 and a = 2 we know from our H-symmetry fixing conditions that c 1,0 ≡ g and c 2,0 ≡ 1 exactly. For a = 1, 2, (5.8) shows that at large g in leading order all coefficients behave in the same way, and this leading order power behaviour is determined by the Hsymmetry 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: A 3 ∼ g 2 /A 2 = g 2 , i.e. g A 3 ∼ g 3 ⇒ n 3 = 3. Similarly: A 4 ∼ g 2 /A 1 = 1, i.e. g 2 A 4 ∼ g 2 ⇒ n 4 = 2.
Next, we can concentrate on the first, leading order coefficients 17 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 nindependent. 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.

JHEP08(2016)061
To guess the exact values of c (0) 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 coefficient. 18 Then one can suspect that the same thing might happen for the cases a = 3, 4. Such an assumption gives analytical predictions for the differences 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) 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 18 In leading order for large g.

JHEP08(2016)061
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.
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 a,n can be perfectly described by polynomials of order 2k. This is why, we make the following conjecture: • The coefficients c (k) a,n are polynomials of order 2k in n.
As a consequence, the polynomials can be given by 2k+1 n-independent parameters, which, for practical purposes, we parametrized as follows: a,k c (m,a) n , n ≥ 1 + δ a,4 , k = 0, 1, 2, . . . , (5.17) where:   Table 3. Numerical values of c We note that in the k = 0 special case, by definition α a,n implies the following series representation for p a (x) at strong coupling: 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) 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 (k) 3 and A (k) 4 of (5.21) are zero for all k ≥ 0.
As a consequence A 3 (g) = A 4 (g) ≡ 0, which implies that besides of the   4 . All values are in the magnitude of the numerical errors.
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,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 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: a,n ∼ 4 2n , α (2n) a,n ∼ 4 2n JHEP08(2016)061  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: (5.25) In the strong coupling limit, (5.25) may fail, if x is close to ±1. In the language of the rapidity 19 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: n c a,n x −2n ∼ g na n K a n √ g x −2 n ∼ g na+ 1 (5.29) 19 Throughout this section, we use the convention, when the branch points are scaled to be located at ±2.

JHEP08(2016)061
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 |x(u 0 + δu)| − √ g ∼ 1 if δu ∼ 1 √ g . 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 20 (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 21 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. 21 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,

JHEP08(2016)061
we made the following proposals for the exact values of the coefficients: We also constructed Pade-approximation like formulae to determine numerically ∆ in the whole range the coupling constant. The description of out Pade-approximation like formulae for the cases S = 4, 6, 8 can be found in appendix G, while their concrete form can be found in the approx.nb notebook file. 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. 22 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. 23

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.
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 Konishi-operator. 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 22 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. 23 At least in case the spin is an even and positive integer number.

JHEP08(2016)061
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 G and ancillary notebookfile approx.nb.) 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 in n. This recognition led us to propose a strong coupling series representation (5.20) for the p a -functions. 24 We argued that (5.20) is an appropriate strong coupling representation of p a (u), if the rapidity u lies outside of an oval domain 25 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, whose strong coupling behavior is discussed in appendix F. 24 The fundamental functions Pa of the QSC method are connected to pa by the simple formula (2.11), this why the results given for pa in the previous sections, can be translated to the language of Pa in a straightforward manner. 25 Here, the rapidity convention is the one, when the branch points are scaled to be at ±2.

JHEP08(2016)061
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.

JHEP08(2016)061
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: The Chebyshev-polynomials can be given by the explicit formula: T n (u) = cos(n arccos u), n = 0, 1, 2, . . .

(B.2)
For practical purposes, we define their slightly modified version: T n (u) = 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, . . .

JHEP08(2016)061
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).

JHEP08(2016)061
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 27 of p a s. In this sense the formulae (4.33) and (4.34) are the periodic analogs of (C.4). Now we show, how to derive (4.33), (4.34) and (4.35) from (4.24), (4.25), (4.26). The derivation of (4.35) goes as follows. One inserts (4.29) into (4.28) and the result into I ij of (4.26). Then evaluating the integrals with the help of the appropriately scaled 28 version (B.5) taken at m = 0, one ends up with (4.35).
To derive (4.33), first one has to rephrase the kernel as an infinite sum: In the sense of (C.2) and (C.3). 28 I.e. u → u 2g substitution in the integral.
To derive (4.34), one should represent ω reg ij by the formula: ω reg ij (u + i 0) = 1 2 (ω ij (u + i 0) + ω ij (u − i 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: We introduce also an integer cutoff Λ X to write the infinite sum as a sum of two terms: 3) The first term in the r.h.s. 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 r.h.s. 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).

JHEP08(2016)061
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): is given by (4.7). The final form of the expansion takes the form: where [. . .] stands for integer part. (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), then 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:

JHEP08(2016)061
We close this appendix with a remark on the usage of the ζ-function in C++. During the development of our C++ code, we recognized that neither double nor long double precisions are not enough to get accurate results at strong coupling. These built in precisions were not enough even to reach some kind of convergence. This is why, we used an arbitrary precision package to C++, called CLN (Class Library of Numbers). In the CLN library ζ(z) is a built in function and it could be used to our purposes. If one uses pure C, or C++, it should be recognized that we need ζ(z) at a finite number of integers. Thus one can compute the necessary values e.g. in Mathematica with high precision and then they can be copied into the C-code and stored in a constant array.
E Some remarks on the implementation of the numerical method at strong coupling In section 4 we explained, that at strong coupling the convergence of the numerical method requires initial values being sufficiently close to the exact solution. To get such good initial values, the value of g was increased in small steps and for a given value of g the initial values were constructed from the previously computed numerical data belonging to neighboring values of g. If one is interested in the strong coupling regime of the solutions, this method seems disadvantageous, since one needs to run the numerical code many times in order to reach the strong coupling regime. One might wonder, whether the strong coupling results of section 5 make it possible to skip this lengthy process and to jump directly in the strong coupling regime?
To answer this question, first let us recall the parameters of the numerical method: • "PRECISION" is the number of decimal digits used in the computations.
• N 0 is the cutoff of the series (2.14) and (2.15).
• l c is the number of discretization points (4.2).
• N I is the cutoff parameter in the large u series (4.3), (4.4) and (4.5).
• N u is an integer being large enough to get Q a|i (u A + i(N u + 1 2 )) within the desired precision.
• h is small shift parameter to compute the derivatives: (4.41).
• "SWITCHOFF" is an integer number which determines the number of iterations beyond which the Levenberg-Marquardt damping parameter is switched off (I.e. λ = 0).
• "RESULT PRECISION" is an integer which stops the run of the numerical algorithm, if S of (4.38) becomes smaller than 10 −RESULT PRECISION .

JHEP08(2016)061
Looking at the list of parameters, it is obvious that a lot of parameters must be kept under control. Each parameter has an effect on the speed of convergence and on the runtime of the numerical code, but the first six parameters have the greatest influence on the convergence, since they strongly affect the magnitudes of the numerical errors within the numerical algorithm. It follows, that not only the inappropriate choice of the initial values can lead to the loss of convergence, but also the numerical errors caused by the inappropriate choice of the parameters of the numerical method.
Let us assume, that one can construct appropriate initial values at strong coupling. Then increasing g in large steps induces the following problem: • The parameters of the numerical method must be changed appropriately, as well. If they are set inappropriately, the numerical algorithm would fail to converge. Then due to the high dimension of the parameter space, it would be very hard (practically almost impossible) to find the optimal choice of parameters which would allow the convergence.
In case g is increased in small steps, the parameters are to be changed by only "small" amounts, ensuring a good control over them. Another practical difficulty is that a strong coupling initial value formula will never allow to construct such accurate initial values, which can be obtained from the high order Taylor series methods of appendix A with g increased in small steps. If an initial value configuration lies farther from the exact solution in the space of solutions, then it requires more number of iterations to reach the exact solution. Thus, it might happen 30 that the increment of g in several smaller units can lead to the numerical solution of the problem at a larger value of g even faster than g was increased in one single step.
Here we would also like to mention some technical limitations concerning the strong coupling implementation of the numerical algorithm. If one uses an ordinary PC with an i7 Intel core to run the numerical code, then at g ∼ 7 the running job uses ∼ 12 − 16 GigaByte of memory and the runtime is ∼ 3 − 4 days 31 if one would like to get the solution with ≈ 20 − digits of precision. So, we think, that ordinary PC's we used, are appropriate for the numerical computations only in the regime g 10 − 14. This gives a technical limitation on the highest available value of g.
We just mention, that the large memory requirement and the long runtime are partly the consequence of the necessity of using an arbitrary precision package for the C++ implementation. In our work we used the CLN (Class Library for Numbers) package. 32 Finally, as closure we would like shed light on some difficulties in constructing appropriate strong coupling initial values from the results of section 5. The main difficulty is that one should construct good initial values for c a,n (g) for all values of n. The conjectured polynomial in n behavior allows one to get good initial values for c a,n (g), when n √ g. 30 The authors had such concrete experience. 31 Without exploiting the possibility of parallelization! 32 http://www.ginac.de/CLN.

JHEP08(2016)061
Taking into account, that the characteristic wavelength of c a,n (g) in n is ∼ √ g, the polynomial behavior allows one to construct good initial values only for the "first wave" of c a,n (g). (See figure 1.) Then at large n, the strong coupling fits of appendix F allows one to construct initial values for c a,n (g). Though the precision of these initial values is not as high as that of the n √ g "polynomial" regime. For the middle n ∼ √ g regime we have no formula to construct good initial values. Nevertheless, based on some extrapolations we tried to construct initial values for this regime, as well. Unfortunately, such a construction of initial values did not lead to the convergence of the numerical algorithm for the values of g we tried.
F The large n behavior of c a,n (g) at strong coupling for the Konishi operator As it is demonstrated in figure 1, the large n behavior ofĉ a,n (g) defined in (5.5) looks as if it was a sine-function with some power-like enveloping curve. This is why, at large n we fitted our numerical data by the following large n Ansatz: The fitting process was done by Mathematica's NonlinearModelFit function. The error bars indicated at the figures correspond to the range of 95% confidence intervals. 33 The fit ranges in n and g were [40,80] and [5.1, 7.0] respectively. Figure 5 shows the numerical values of Λ a (g) together with the error bars at the discrete values of g in the range [5.1, 7.0]. The numerical data shown in figure 5, suggests that the value 34 of Λ a (g) is independent of the value of the index a. On the other hand at strong coupling we fitted the data points for Λa(g) √ g with a series in 1 √ g . The first, constant coefficient of the fits proved to be stable under changing the different parameters of the fitting process. Thus we made the following conjecture: with a common coefficient c 0 = 4.35 (5) and the dots stand for terms negligable for g → ∞.
Another important parameter of the large n behavior is the exponent a (g). The numerical data together with the estimated error bars can be seen in figure 6. Looking at figure 6, one can immediately see that the errors are much larger, than those of the characteristic wavelengths. Figure 6 suggests that the exponents a (g) tend to constant values when g → ∞. The 4 constant values are close to each other, but because of the large error bars, one cannot conclude that they would be the same.
The data for the phase factors ϕ a (g) are shown in figure 7. Figure 7 suggests, though with large error bars, that the phases also tend to constant values at g → ∞, but these constants depend on a. 33 We just note, that the actual errors can be even larger due to systematic errors. 34 At least within the numerical precison of Λa(g), which is somewhere between 10 −3 and 10 Here the data are so close to each other, that one cannot make a difference between them. The plotted line corresponds to the ∼ √ g fit of the numerical data. Finally, we close this appendix with some remarks concerning the estimated errors of the parameters of the large n Ansatz (F.1). From figures 5, 6 and 7 it can be seen that the parameters of the large n Ansatz (F.1) cannot be determined with high accuracy from the numerical data of c a,n (g). There are 3 reasons for that: • In the numerical method the series representaions (2.14) and (2.15) are truncated in n. These truncations lead to increasing relative errors for c a,n with large n. Unfortunately, this is the regime, where the Ansatz (F.1) is used.
• Another source of the large errors might be that the intuitively chosen Ansatz (F.1) is not the 100% correct form of the asymptotic large n behavior ofĉ a,n .
• Even if the intuitive Ansatz (F.1) gives the correct leading order large n behavior ofĉ a,n , the numerically accurate determination of its parameters would require the knowledge of correction terms as well. Unfortunately, the analytic form of these terms are unkown.
G Pade-approximation like formulae for the anomalous dimensions   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 35 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]: ∆(g) = (g 2 + g 2 b ) 1/4 a 0 + a 1 h + . . . a n h n 1 + b 1 h + . . . b n h n . Where h = g 2 g 2 +( 1 4 ) 2 and g b is a suitable constant, whose value was chosen to be 2 in the cases of S = 2, 4 and it was set to be 1 in the cases of S = 6, 8.
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 36 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" 37 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 8.
Because of the small magnitude of the deviations, we show separately the residual plot of the data in figure 9. Figure 9 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 10 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. 36 I.e. ∆(0) = L + S, where L = 2 for twist-2 operators and S = 2 for the Konishi-state. 37 Interested readers can gain more information about this function in the help of Mathematica. JHEP08(2016)061  Analogously to (G.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 h-dependent factor in (G.1). The concrete form of the bulky approximation formulae can be found in the approx.nb notebook file attached to the paper.

H Various tables of numerical data
This appendix contains some tables of numerical data which demonstrates that the coefficients c (k) a,n of (5.7) are polynomials of order 2k in n.